Exploring the Compositional Space of High-Entropy Alloys for Cost-Effective High-Temperature Applications

High-entropy alloys (HEAs) are nearly equimolar multi-principal element alloys, exhibiting exceptional thermal and mechanical properties at extreme conditions such as high-temperatures and stresses. Since the first discovery and early conceptualization of conventional HEAs nearly two decades ago, HEAs with far-from-equimolar compositions have attracted substantial interest to provide a broader range of material properties and to adjust price fluctuations and availability of commodities. Here, we present a first-principles investigation of non-equimolar chromium-manganese-iron-cobalt-nickel (CrMnFeCoNi) HEAs and effects of molybdenum (Mo) and niobium (Nb) substitutions on cost, phase stability and solubility, and mechanical and thermal performance up to 1000 K operational temperature. Virtual-crystal approximation is used to expediently approximate random solid solutions at the disordered mean-field limit. Using multi-objective metaheuristics built on a first-principles database, golden compositions are predicted for thermally well-insulated components and effective heat sinks. Replacing Co with Fe lowers commodity costs without hindering phase stability and solubility. Lower Ni concentration leads to lower thermal conductivity, indicating better thermal insulation, while reducing Mn concentration significantly increases the thermal conductivity, indicating better performing heat sinks. Moving away from equimolar ratios commonly increases the thermal expansion coefficient, which could generate higher thermal stresses. Nb and Mo substitution always lead to substantially higher commodity cost and density but with an increment in the mechanical performance due to solid-solution hardening. However, alloying with Mo and Nb is the only compositional space that reduces the thermal conductivity and thermal expansion coefficient.


INTRODUCTION
High-entropy alloys (HEAs) are single-phase multi-principal element alloys exhibiting solid solution and are conventionally mixed with equimolar fractions . Among conventional HEAs, the face-centered cubic (FCC) chromium-manganese-iron-cobalt-nickel (CrMnFeCoNi), also called Cantor alloy, (Cantor et al., 2004) has been extensively studied due to its remarkable mechanical properties such as high yield strength (Kim et al., 2020), tensile strength (Gali and George, 2013;Otto et al., 2013), hardness (Liu et al., 2013), and fracture toughness (Gludovatz et al., 2014). The early concept behind the formation of HEAs has been mostly attributed to high configurational entropy (S conf ) due to configurational disorderliness (Yeh et al., 2004). By definition, S conf is the highest for equimolar mixture, for which it substantially reduces mixing Gibbs free energy. However, there has been a rapidly growing literature on the development and manufacturing strategies of non-equimolar HEAs (Miracle and Senkov, 2017), in particular, Cantor-like HEAs (Cantor, 2021).
Non-equimolar HEAs provide a broad compositional space for multi-objective materials design for diverse engineering applications for extreme conditions such as aerospace (Dada et al., 2021a), nuclear reactors (Tong et al., 2005), and biomedical (Varalakshmi et al., 2010) applications. In such applications, conventional alloys commonly have severe limitations to achieve more cost-effective, long-living, and more green designs (Dada et al., 2021b). In particular, low Co non-equimolar HEAs are an appealing system for diverse applications and the replacement of high-strength steel components .
Ad-hoc design of engineering materials is a multi-faceted task requiring consideration of material properties, manufacturability, operational conditions, compatibility to the current design, and environmental impact (Mouritz, 2012;Edwards, 2014). Multiobjective metaheuristics (Liu et al., 2020) serve as powerful tools for compositional optimization during materials selection (Ashby, 2000). There have been previous works such as composition selection based on commodity availability and cost (Fu et al., 2017) and phase-stability prediction based on thermodynamic multi-objective criteria (Gheribi et al., 2018). However, such works are severely limited to a single aspect of materials selection due to a lack of available materials data. Firstprinciples simulations provide tractable tools to calculate material properties highly suitable for multi-objective optimization (Liu et al., 2020).
This work presents a systematic first-principles investigation of composition-and temperature-dependent material properties of non-equimolar Cantor-like HEAs. We also investigate the effects of Nb and Mo substitutions. First-principle materials database and available commodity cost are used to perform multi-objective optimizations within fundamental materials objectives and constraints for high-temperature applications (T > 800 K), particularly thermally well-insulated components and effective heat sinks.

MATERIALS AND METHODS
The equimolar CrMnFeCoNi, labeled by X 5 , serves as the reference system. Two compositional spaces are derived starting from X 5 . The first compositional space was spanned by non-equimolar five-principal-elements (5-PEs) HEAs. It was sampled using a discrete mesh of compositions, obtained by individually varying molar fractions of each PE in between x 0.05 − 0.3, where x is the atomic molar fraction while keeping the remaining PEs equimolar among themselves. For instance, the non-equimolar 5-PE HEAs are Cr x (MnFeCoNi) 1−x , labeled by Cr x X 4 1−x , where x = 0.05, 0.1, 0.15, 0.25, 0.3 when varying Cr concentration. The second compositional space was spanned by 6-and 7-PE HEAs. It was sampled using a discrete mesh of compositions, obtained one-by-one by substituting half each PE of X 5 with Nb and Mo. For instance, the 6-and 7-PE HEAs are Cr 0.1 X 4 0.8 Nb y Mo 0.1−y where y 0.1, 0.075, 0.05, 0.025 when the molar fraction of Cr is halved. Relevant symbols and molar fractions of available compositions are listed in Supplementary Tables S1, S2.

Objectives and Constraints for Optimization
Inspired by the most common ambient conditions for high-stress and high-temperature applications, four main objective categories were determined for multi-objective optimization, schematically illustrated in Figure 1. The main objectives and constraints (whenever present) are also summarized in Figure 1. We focus on fundamental material properties rather than designspecific performance metrics as they are more suitable for general assessments. Moreover, they are also used to determine the most common specific performance metrics (Ashby, 2000).
The first objective category is economics, aiming to minimize commodity price and solid density (ρ). Economics, associated with materials design, is far more complex due to manufacturing, repair, replacement, and recycling costs. While these are highly dependent on the particular design, commodity cost is a fundamental factor. It is particularly crucial when considering large fluctuations in unit prices, shown in Supplementary Figure  S1. Solid density is an indirect contributing factor, particularly for transportation technologies, determining fuel consumption. For instance, estimations indicate that a 1 kg reduction on an averagesized airplane may reduce fuel consumption by ∼ 800 liters per year (Mouritz, 2012). Beyond reducing operation costs, lower density also reduces carbon footprint and associated costs.
The second objective category is related to the phase-stability and solubility of a multi-principal-elements alloy. Relative phase stability of non-equimolar compositions with respect to X 5 was assessed by Gibbs free energy (GFE) difference (ΔG free ). In addition to the GFE difference, likelihood of solid-solution formation is assessed by valence electron concentration (Δn v ), Mulliken electronegativity (Δχ), and atomic-size mismatch parameter (Δ(δr)) differences (Qin et al., 2019;Rojas et al., 2022). A limit of 5% increase G free compared to that of X 5 is chosen as the upper-bound constraint to ensure that comparable phase stability is attained. Rojas et al. (2022) have shown that good solubility of Mo and Nb solute atoms was achieved when Δn v and Δχ were below ±2%. These conditions are imposed as constraints during the optimizations. Lastly, minimizing Δ(δr) serves as an optimization objective when Mo and Nb are considered.
The third objective category is mechanical properties for which Young's modulus (E), yield strength (σ y ), and Pugh ratio (ratio between bulk (K) and shear (G) moduli) serve as objectives. E and σ y commonly serve as figures of merits (FOMs) for mechanical strength (Ashby, 2000), while K/G serves as a figure of merit for ductility, thus leading to better workability (Boucetta, 2014). For most engineering applications, these FOMs are desired to be maximized.
The last objective category is related to the thermal behavior, expediently assessed through thermal conductivity (κ) and thermal-expansion coefficient (β). A low κ is preferred for thermal insulation such as in jet-engine components (Mouritz, 2012), while a higher κ is required for effective heat sinks to dissipate the excess of heat rapidly. In general, a low β is highly preferable to avoid thermal dilatation in integrated components.

Theoretical Methodology for First-Principles Material Properties
The Kohn-Sham density-functional theory (KS-DFT) (Hohenberg and Kohn, 1964;Kohn and Sham, 1965) with semi-local electronic exchange-correlation functionals (Kohn and Sham, 1965;Perdew et al., 1992) is an almost-universal approach to study ground-state properties of solids. Its perturbative extension, namely, density-functional perturbation theory (DFPT) (Gonze, 1997;Gonze and Lee, 1997), also provides a parameter-free approach to simulate lattice dynamics. These theories offer parameter-free approaches to calculate elastic and anharmonic thermodynamic properties (Corso, 2016;Palumbo and Corso, 2017).
Despite rapidly growing computational resources and state-ofthe-art algorithms in recent decades, practical implementations of the approximate KS-DFT are still severely limited by system size and spatial complexity, particularly for metallic systems (Ponga et al., 2016;Ponga et al., 2020). It becomes more infeasible for HEAs as they require many well-converged super-cell simulations to achieve an accurate statistical representation of configurational disorderliness. Virtual-crystal approximation (VCA) (Bellaiche and Vanderbilt, 2000) offers a less robust yet expedient method to simulate such systems. A linear mixing scheme is applied to atomic pseudo-potentials (Heine, 1970) to achieve a virtual representation at the disordered mean-field limit (Bellaiche and Vanderbilt, 2000). Although VCA has been commonly deemed over-simplified, it is entirely accurate when calculating elastic constant, simplephase transformation, and thermodynamic functions of HEAs (Tian, 2017).
Random solid solutions such as HEAs do not exhibit longrange ordering due to spatial disorder; however, they still may have short-range ordering (SRO), which can play crucial roles in the local chemical environment and lattice distortions (Miracle and Senkov, 2017). Such local effects may significantly contribute to superior thermal stability, sluggish diffusion, and exceptional mechanical strength Miracle and Senkov, 2017). SRO may also be necessary for capturing fine details of electronic band structures (Maurizio et al., 2003). By definition, VCA does not include any SRO. There have been considerable efforts to include SRO within VCA by using tight-binding modellike approaches to explicitly include interactions of unlike first and second neighbors (Baldereschi and Maschke, 1975;Porod and Ferry, 1983). Despite their relative simplicity, such corrective approaches are unsuitable for high-throughput materials research due to additional computational costs.

Figures of Merits for Phase Stability and Solubility
GFE serves as a universal FOM when comparing the relative phase stability of solids. For non-magnetic HEAs, it is given by the sum of the enthalpy of mixing (H mix ), the electronic Helmholtz free energy (F el ), the vibrational Helmholtz energy (F vib ), and the configurational-entropy contribution ( − TS conf ). H mix is commonly approximated using the enthalpy of mixing of binary bulk metallic glasses (Youssef et al., 2015), while F el and F vib can be approximately calculated within the KS-DFT-based methods (Wang et al., 2004;Shang et al., 2010), and S conf is given by FIGURE 1 | Schematic illustration of multi-objective optimization for thermal insulation and heat sinks, operating at high-stress and high-temperature. Abbreviation list: ρ: Solid density, G free : Gibbs free energy, n v : Valence electron concentration, χ: Mulliken electronegativity, δr: Atomic-size mismatch parameter, E: Young's modulus, σ y : Yield strength, K: Bulk modulus, G: Shear modulus, κ: Thermal conductivity, β: Thermal expansion coefficient.
Frontiers in Materials | www.frontiersin.org January 2022 | Volume 8 | Article 816610 where R is the gas constant, M is the number of PEs, and c i is the molar fraction (Santodonato et al., 2015). By definition, Eq. 1 reaches its maximum for equimolar mixture of M PEs. Antiferromagnetic ordering of Cr atoms against significant local magnetic moments of Fe and Mn atoms commonly leads to weak paramagnetism in X 5 and its FCC derivatives (Tian et al., 2013;Schneeweiss et al., 2017;Šebesta et al., 2019). Hence, the magnetic contribution to GFE due to the total magnetic moment is negligible. For instance, it is −3.23 × 10 −7 eV per unit cell for a magnetic field of 500 Oe which induces a total magnetic moment of ∼ 4 × 10 −3 emu on a sample of X 5 with a dimension of 10, ×, 2 × 0.02 mm at 300 K.
We define n v as the electron number at outermost shell per unit volume to include temperature effects due to thermal expansion. It is also more convenient when comparing different phases of an alloy with the same ratio of PEs. Within KS-DFT, χ is approximately equal to the negative of chemical potential (Parr et al., 1978), which is a temperature-dependent quantity. Finally, δr is simply given by where r i and r avr are the atomic radius of PEs and their arithmetic averaging, respectively.

Approximate Mechanical Properties
It is a challenging task to calculate technologically relevant mechanical properties such as Vickers hardness, yield strength, or ductility/brittleness, for which even experimental measurements are highly dependent on morphology, choice of indenter, the direction of loading forces relative to the microstructure orientation, and impurity of samples (Chen et al., 2011). Despite that, there are some expedient macroscopic semi-empirical connections between such quantities and the HEA's second-order elastic tensor, which is well-defined in the linear elastic regime (Corso, 2016). Pugh ratio, given by the ratio between bulk and shear moduli (K/G), is commonly used as a FOM for ductility and workability. A Pugh ratio above 1.75 indicates good workability (Boucetta, 2014). Vickers hardness is commonly used as an industry-standard measure for wear resistance. It is approximated by H Teter V 0.151G (Teter, 1998), or H Tian V 0.92(G/K) 1.137 G 0.708 (Tian et al., 2012). The former one tends to underestimate, whereas the latter expression tends to overestimate the material hardness. Thus, we used their averaging to approximate Vickers hardness, given by H V 0.5(H Teter V + H Tian V ). Following this approximation, the semi-empirical relation between H V and yield strength (σ y ) is given by H V ≈ 3σ y (in MPa units) (Tiryakioğlu, 2015). Within this definition and in the context of our calculations, σ y can be interpreted as the lattice friction as it does not account for dislocation density, grain-size, solid-solution hardening (SSH), and precipitation hardening, which can synergetically increase σ y . In 5-PE CrMnFeCoNi HEAs, it has been shown that the lattice friction value is in the neighborhood of ∼ 164 MPa (Zhao and Nieh, 2017), which is considerably higher than in other single-phase metals such as Ni.
Temperature dependence of macroscopic elastic properties was carried out through temperature dependence of K, E, and G. Temperature-dependent K was obtained from the Murnaghan equation of state at each temperature (Palumbo and Corso, 2017), while temperature-dependent E and G were carried out using temperature-dependent longitudinal and transverse sound velocities, respectively, and solid density (Farraro and Mclellan, 1977).

Approximate Solid-Solution Hardening due to Nb and Mo Substitution
Introducing additional PEs, in particular, with different atomic radii such as Nb and Mo to X 5 , commonly improve mechanical strength through SSH (Fleischer, 1963;Labusch, 1970). SSH can be expediently estimated by the sum of individual SSH of each solute atom, given by (Toda-Caraballo and del Castillo, 2015), where c i is the molar fraction of i th solute atom, α is the theorydependent parameter which is commonly 3 < α < 16 and α > 16 for screw and edge dislocations, respectively.
G is the composition-dependent shear modulus, Z is the fitting constant, and a is the composition-dependent lattice parameter.

Approximate Thermal Properties
Thermal-expansion coefficient is a key FOM to assess thermal dilatation of components exposed to high temperatures. It can be calculated within the quasi-harmonic approximation starting from well-converged electronic structure at different volumes (Palumbo and Corso, 2017). Thermal conductivity was used to assess the ability of dissipating excess heat in solids. For metallic systems, it is predominantly mediated by free electrons and phonons, given by their respective contribution, symbolically κ κ el + κ lat . Electronic thermal conductivity can be approximated by where k B is the Boltzmann constant, N(E F ) is the density-of-state at the Fermi level, V is the unit-cell volume, v F is the Fermi velocity, λ is the electron-phonon mass enhancement parameter, and τ p is plasmon lifetime (Orhan and Ponga, 2021). Lattice thermal conductivity (in W·m −1 ·K −1 unit) is given by the semiempirical formula: Frontiers in Materials | www.frontiersin.org January 2022 | Volume 8 | Article 816610 where m avr (in kg unit) is the average mass, V 0 (in m 3 unit) is the equilibrium volume, N is the number of atoms in unit cell, Θ D is the Debye temperature, and γ is the Grüneisen parameter (Morelli and Slack, 2006). We refer the reader to the recent work of Orhan and Ponga (2021) for details on obtaining κ within fully first principles.

Computational Details
Ground-state and lattice-dynamics simulations were performed using the Quantum ESPRESSO (QE) software (Giannozzi et al., 2009(Giannozzi et al., , 2017. Virtual-atomic pseudo-potentials within VCA were generated using the SG15 optimized norm-conserving Vanderbilt scalar-relativistic pseudo-potentials of PEs (Hamann, 2013;Schlipf and Gygi, 2015). A generic FCC crystal was used for initial crystallographic information. Variable-cell optimizations were performed with a 150 Ry kinetic-energy cutoff for wavefunctions (E cut ), and a choice of smearing parameter for the Marzari-Vanderbilt cold smearing (Marzari et al., 1999) equivalent to 300 K at a 12 × 12 × 12 Monkhorst-Pack-equivalent (Monkhorst and Pack, 1976) uniform Brillouin zone with a 10 −7 Ry total energy, 10 −6 Ry/a 0 total force convergence, and 10 −10 self-consistency thresholds. Equilibrium crystal structures were used to perform three sets of simulations. For subsequent simulations, E cut was set to a common converged value of 100 Ry. Phonon dispersion was calculated on a uniform 4 × 4 × 4 mesh. In the first set, electronic structures, elastic properties, and anharmonic thermodynamic properties were calculated on a shifted 12 × 12 × 12 uniform Brillouin zone on a temperature grid in between 0 − 1000 K with a step size of 5 K. In the second set, a series of sequential simulations were performed to obtain electron-phonon coupling on an electronic temperature grid in between 100 − 1000 K with a step size of 100 K. Ground-state simulations on a highly dense Brillouin zone sampling of 36 × 36 × 36 were performed to be used for interpolating electron-phonon coupling coefficients in the following simulations. Subsequent simulations in the second set were performed with a 12 × 12 × 12 uniform Brillouin zone. High-symmetry points on the first Brillouin zone of primitive FCC symmetry were used for a weighted sum of electron-phonon coupling coefficients to obtain an electron-phonon mass enhancement parameter (λ). In the last set, electronic band structures were calculated on a 36 × 36 × 36 uniform Brillouin zone to be used as input database for calculating Fermi velocities using an in-house code available online (Orhan, 2020).

Multi-Objective Optimization
In this work, our composition spaces represent design spaces where each available combination is a particle. We use the particle-swarm optimization (PSO) (Kennedy and Eberhart, 1995), a stochastic multi-objective metaheuristics method that has been extensively used in diverse optimization problems such as nano-structure design (Jiao et al., 2021) and data processing (Song et al., 2021). A multi-objective function with goals shown in Figure 1 was defined using the constant-weighting method (Lian, 2010). Thus, the multi-objective problem was solved as a linear combination of weights and single-objective function by satisfying the imposed phase stability and solubility constraints through the feasibility method (He and Wang, 2007). Compositions were promoted to explore regions where feasible solutions exist. For our particular problem, feasible-solution regions satisfied phase-stability and solubility conditions.

Interpolating Compositional Spaces Using the Higher-Order Maximum-Entropy Scheme
Conventionally, PSO is performed on continuous design spaces while available compositions are obtained on a discrete mesh of molar fractions at different temperatures. Thus, a fitting scheme is required at the end of each iterative solution to ensure the multiobjective functions represent the original discrete composition space, spanned by the available first-principles database. Two interpolation schemes were tested.
The higher-order maximum entropy scheme (HOLMES) is an approximation method designed to operate over unstructured multidimensional data (Bompadre et al. 2012). Within HOLMES, a shape function provides a balance between maximizing information entropy for least-biased statistical inference and minimizing shape-function width for computational efficiency. Moreover, HOLMES allows higher orders of polynomial consistency, enabling high-order interpolation of sparse data. HOLMES was explicitly designed to operate over unstructured multidimensional data. It is robust for situations where data points are spaced in widely varying manners and do not require complicated tinkering of its shape parameter to achieve consistent results. Additionally, because the shape functions achieve maximum entropy for their width, the resulting approximation is more resilient to uncertainty in the data than other methods. This property makes HOLMES a good choice in scenarios where clumping of data points may occur, such as exploring a design space for optimization. However, being an approximate rather than an exact interpolant means the HOLMES scheme will not necessarily duplicate data points values.
In practice, HOLMES is a set of shape functions, w a (x), centered around a corresponding set of data points or nodes, x a . Interpolation is performed as per any standard interpolation scheme. For some function u(x), the HOLMES approximant is where N a (x) is the number of nodes within a neighborhood of the query point x. Within the context of this work, the functions u(x a ) are the composition-dependent response functions of material properties of HEAs, listed in Figure 1, at each temperature point. Smooth response functions were obtained a high-dimensional space as required by the optimization process.
To compare the optimization procedure results, we also implemented the radial basis function (RBF) interpolation. RBF interpolation is a well-known alternative interpolation technique operating over unstructured multidimensional data (Buhmann, 2003). RBF uses kernel functions centered at each data point to create a system of equations with the shape functions as unknowns. Polynomial consistency of the desired dimension and order can be achieved by expanding the system of equations to include polynomial coefficients and suitable orthogonality constraints Wendland (2005). The reader is referred to Fasshauer (2007) for practical implementation details. RBF interpolation in its standard form is an exact interpolation, meaning it duplicates known values at data points. For some choices of kernel and data set, RBF has even achieved spectral rates of convergence (Madych and Nelson, 1992). However, in scenarios with extreme clumping of data points, the linear system may become ill-conditioned and ruin interpolation accuracy. Furthermore, performance is highly dependent on the choice of kernel and kernel parameters. The compact Wendland C2 kernels from Wendland (2005) produce sparse matrices and are less parameter-sensitive than most other kernels. Given the small size of the data set and its regular spacing, the kernel support radius is allowed to encompass all data points.

Computational Details
Multi-objective optimizations were performed using the numerical computing software MATLAB(R2021a) following (Isiet and Gadala, 2019). Global solutions were obtained by 100 independent runs with each 1000 iterations to ensure the robustness of the stochastic process. 100 randomly generated compositions were used for the discrete representation of composition space. Multi-objective functions were evaluated for each temperature.
Despite evaluating the same amount of multi-objective function, the computational demands of RBF and LME differ due to the nature of the methods. For a single processor, the ratio of time taken for a single evaluation of a 5-PEs composition between LME:RBF was 0.013. However, when evaluating a large number of compositions, as is in the case here, the computational demands for LME increase sharply. For example, for 100 evaluations, the ratio LME:RBF ratio is 0.6181, and for 1000 evaluations, the ratio increases to 5.0426. Therefore, it could be said that LME is more appropriate for problems where limited sampling is required, while RBF performs better when dealing with a large number of evaluations across a compositional space. For multi-dimensional problems with low-order nonlinearity, like here, RBF performs better as it is computationally more efficient, but in problems with high-order nonlinearity, LME has been shown to interpolate data with a higher accuracy (Fan et al., 2018).

RESULTS AND DISCUSSION
We start with relative phase stability and solubility of nonequimolar HEAs with respect to the reference system (X 5 ). In Figure 2, three general trends arise in ΔG free of 5-PE nonequimolar HEAs. The first one is that moving further away in the equimolar case (x 0.2) leads to larger |ΔG free |. The second one is that this behavior is not symmetric, and ΔG free increases more rapidly for compositions that increase the GFE. The last one is that increasing (decreasing) molar fractions of lighter (heavier) PEs decrease (increase) ΔG free . As Fe is at the middle in terms of atomic number and mass, its ΔG free is relatively small, almost within the practical KS-DFT accuracy threshold of ∼ 1 meV.
In Supplementary Figure S2, the four terms contributing to ΔG free are shown. ΔH mix and ΔF el are in the orders ≤ 10 −1 meV which are below practical KS-DFT accuracy threshold and are negligible. Thus, ΔF vib and ΔS conf determine relative phase stability. By definition, S conf is lower for non-equimolar 5-PE cases compared to X 5 ; thus, − (TΔS conf ) is always positive. This leaves ΔF vib to determine the sign of ΔG free . In Supplementary Figure S2, ΔF vib takes positive values for reduced molar-fractions of lighter elements (i.e., Cr and Mn), while it takes negative values for heavier elements (i.e., Co and Ni). It becomes insignificant in the case of Fe x X 4 1−x . When ΔF vib is positive, it works with ΔS conf in the same direction. Thus, |ΔG free | is asymmetrical with respect to deviation from the equimolar case as shown in Figure 2.
The trend in ΔF vib can be traced back to phonon density of states (PDOS), as shown in Supplementary Figure S3. Increasing (decreasing) Cr or Mn (Co or Ni) shifts PDOS to lower (higher) energies and leads to higher (lower) F vib and lower (higher) ΔF vib . Tracing back further, ΔF vib follows an almost identical trend seen in the average atomic mass difference (Δm avr ), shown in Supplementary Figure S4, except Fe x X 1−x for which changes are negligible. This indicates that interatomic force constants (IFC) (hence, the second derivative of energies with respect to ionic displacements) do not change substantially while shifts in PDOS are mostly due to Δm avr as its square-root is used to normalized IFC before calculating phonon frequencies (Baroni et al., 2001).
Supplementary Figure S5 shows the changes in χ, Δn v , and Δ(δr) compared to the X 5 HEAs. The relative small changes, all below 2%, suggest high-solubility of these compositions. Moreover, moving further away from the equimolar case leads to larger solubility FOMs which provide a large compositional space for optimization for 5-PE cases.
In the case of additions of Nb and Mo, despite ΔG free being always negative in Figure 3, the underlying mechanisms are more complex. The main factor leading to negative ΔG free is higher S conf due to additional element(s), resulting in large negative − (TΔS conf ), as shown in Supplementary Figure S6. Furthermore, ΔF vib is also negative except in the cases of Cr 0.1 X 4 0.8 Mo 0.1 and Cr 0.1 X 4 0.8 Nb 0.025 Mo 0.075 for which − TΔS conf is large enough to compensate. Unlike previous cases, ΔF vib is not simply determined by Δm avr . A more complex behavior is present in PDOS, shown in Supplementary Figure S7, such as larger energies shift to lower energies when substituting a heavier element with Nb and Mo despite lower Δm avr as shown in Supplementary Figure S8. For instance, despite overall highest Δm avr , there is almost no shift in the case of Cr 0.1 X 4 0.8 Nb y Mo 0.1−y . On the other hand, for a given substituted PE, higher Δm avr leads to smaller shifts and lower energies. It indicates that effects of Δm avr are over-compensated by higher IFC.
In Supplementary Figure S9, three FOMs for solubility are shown for T 1000 K when No and Mo are used as solute atoms. δχ and Δn v have weak temperature dependence while Δ(δr) is independent of temperature. These FOMs are desired to be kept as low as possible to reduce affinity for phase segregation . For additions of Nb and Mo, the solubility FOMs indicate that some phase segregation in the form of precipitations may not be avoidable .

Temperature Dependence of Material Properties
Temperature dependence of mechanical and thermal properties has to be taken into account when optimizing compositions for a range of operational temperatures. In Supplementary Figure S10 and Supplementary Figure S11, material properties in between the room temperature (300 K) and 1000 K are presented for 5-, 6-, and 7-PE HEAs. Prior to an in-depth analysis, there are some general trends arising throughout compositions. In general, ρ and E do not have strong temperature dependence; Pugh ratios, σ y , and β exhibit moderate temperature dependence, while κ el and κ lat vary substantially with temperature. Since σ y ∼ ∝ G 2 /K while the Pugh ratio is K/G, they always have opposite trends with respect to temperature. Also, β mostly determines temperature dependence of E and G through longitudinal and transverse sound velocities, respectively, as ρ has a weak temperature dependence. Temperature dependence of relevant quantities in Eq. 6 and Eq. 7 are presented in Supplementary Figure S12 and Supplementary Figure S13 to investigate strong temperature dependence in κ el and κ lat . Strong temperature dependence of κ el is mostly due to τ p in Eq. 6, such as however, λ is in the order of 1. In the case of κ lat , it is almost entirely due to κ lat ∝ T −1 , as for which f[γ]/γ is mostly constant.

Composition Dependence of Material Properties at High Temperature
We shift our focus on the high-temperature case where T 1000 K ( ∼ 725°C). In Figure 4, materials objective for nonequimolar 5-PE HEAs alongside X 5 are shown with re-scaled by-values of X 5 (black solid lines). Increasing Fe content or decreasing Co content is the most reliable option to reduce commodity cost. Remarkably, changing Fe content does not significantly change material properties. In the case of Co reduction, κ decreases while β becomes larger, hinting better thermal insulation at the expense of higher thermal expansion. Tuning Cr, Mn, or Ni content gives wider range control over σ y , the Pugh ratio, κ, and β; however, it does not lead to any significant cost reduction. In the case of κ, in particular for Ni x X 4 1−x , it varies over a wide range, in between ∼ 0.25 − 1.25 times compared to that of X 5 .
Results for additions of Nb and Mo in the composition are shown in Figure 5. Substituting any other PE by Nb drives the cost up as it is the most expensive metal (shown in Supplementary Figure S1). Replacing Co or Mn by Mo keeps the cost almost fixed. Substituting Cr with Nb or Mo does not have any significant benefit except slight reduction of thermal conductivity and thermal expansion coefficient. Ni 0.1 X 4 0.8 Nb y Mo 0.1−y is also not quite favorable as they have lower strength at a higher cost. Despite higher mechanical strength, Fe 0.1 X 4 0.8 Nb y Mo 0.1−y is not a good choice as well due to its higher thermal expansion and lower ductility. Among remaining options, Co 0.1 X 4 0.8 Nb y Mo 0.1−y is the most promising as it provides a compositional range for simultaneously improving mechanical strength and ductility while keeping thermal expansion almost constant.
In Supplementary Figure S14 and Supplementary Figure  S15, composition-dependent quantities, used to calculate κ el and κ lat in Eq. 6 and Eq. 7, are shown to further analyze this trend. For the 5-PE cases, κ lat is one order higher than κ el . Thus, κ shows the same trend with κ lat in Supplementary Figure S14 regardless of how κ el behaves. The most visible clue of this behavior is in Ni x X 4 1−x in Supplementary Figure S14, where the term f[γ]/γ in Eq. 7 determines this trend. This term monotonically decreases for increasing γ between 0.5 − 4.5, representing the approximate lower and upper bounds for HEAs in this work. Thus, γ predominantly determines trends in κ of 5-PE HEAs, and higher γ leads to lower κ. For Nb and Mo substitutions, the general trend in κ is subtle changes for Cr and Mn while there are substantial reductions for Fe, Co, and Ni. In the case of substituting Cr, Mn, and some cases of Fe, κ lat does not change significantly, and lower κ el leads to quite subtle changes in total thermal conductivity.
In other cases, κ lat and κ el are substantially lowered by Mo and Nb substitution for which κ lat sets the trends in κ. Similar to the previous case, changes in κ lat (or lack of) are driven by γ as other terms in Eq. 7 are not strongly composition-dependent. Since τ p ∼ ∝ λ −1 , κ el ∼ ∝ λ −3 while κ el ∝ N(E F )/V and κ el ∝ v 2 F . As shown in Figure 5, N(E F )/V is almost constant throughout the compositional space and v F varies between higher and lower energies, while higher λ compared to that of X 5 leads substantially lower κ el . One step further, this trend can be traced back trend in PDOS (Supplementary Figure S7) for which more shifts to lower energies lead to lower κ el due to higher λ ∼ ∝ PDOS(E)/E 2 .

Solid-Solution Hardening due to Nb and Mo Substitution
It is now interesting to analyze the effect of strengthening of solutes (i.e., Mo and Nb) into the HEAs. Remarkably, their large radii difference will pin dislocations near the solutes, and this can be quantified using the approach presented by (Toda-Caraballo and del Castillo, 2015). We also point out that this effect is added to the one computed through the elastic interactions -also known as lattice frictionof the matrix presented in the previous sections. Figure 6 shows the effect of SSH for several compositions of Nb and Mo at 1000 K Figures 6A,C and the temperature dependence of SSH when the same molar fraction of Mo and Nb is used (y 0.05) Figures 6B,D. SSH for screw (top panel) and FIGURE 5 | Composition dependence of material properties of non-equimolar six-and seven-principal-elements high-entropy alloys, shown by their ratio to those of CrMnFeCoNi (X 5 ), at T 1000 K.
Frontiers in Materials | www.frontiersin.org January 2022 | Volume 8 | Article 816610 edge (bottom panel) dislocations were obtained by setting the theory-dependent parameter α in Eq. 4 to 4 and 16, respectively, and Z 1/700 (see Supplementary Figure S16 for the full composition-temperature space). SSHs were calculated by adding up all the individual contributions due to Nb and Mo contents, shown in Supplementary Figure S16. In Figures 6A,C, SSHs generally reach a maximum between y 0.025 or y 0.075, hinting that 1 : 3 or 3 : 1 ratio of Nb:Mo is most favorable to have maximum additional hardening. This composition also stabilizes the changes in electronegativity with respect to the base alloy, as shown by Rojas et al., (2022). In Figures 6B,D, SSHs have subtle temperature dependence which is a general trend throughout compositional space as seen in Supplementary Figure S16. Temperature and composition dependence is mostly determined by shear modulus rather than lattice parameters (shown in Supplementary Figure S17), which lead to a SSH optimal value when y 0.025 − 0.075 (see Supplementary Figure S18).

Optimized Candidate Compositions
Compositional space (design space) for the 5-PE HEAs was discretely sampled by eleven points for each PE, whereas it was halved for the original PEs with Nb and Mo at four different sampling points in the case of 6-and 7-PEs. Thus, we constrained compositional spaces to the discrete compositions during optimizations for which the 5-PEs case had a large design space compared to the 6-and 7-PE cases. With that in mind, two general cases were studied within two interpolation schemes. Comparing Supplementary Figures S19, S20, the RBF interpolation scheme performs better than HOLMES for the 5-PEs compositional space. This is not surprising as HOLMES maximizes information entropy and thus smooth data out. Thus, they result in different optimal compositions, listed in Supplementary Tables S3, S4 for thermally well-isolated  components, and Supplementary Tables S5, S6 for effective heat sinks. In the case of 7-PE HEAs, the two interpolation schemes perform quite similarly, as shown in Supplementary  Figures S21, S22, resulting in the same optimal compositions, as listed in Supplementary Tables S7, S8 for thermally well-isolated components, and Supplementary Tables S9, S10 for effective heat sinks. Thus, the RBF interpolations scheme is used to determine optimal compositions.

Thermally Well-Insulated Components
The first case is where a component is operating at a high temperature and required to be thermally well-insulated to protect other components. In this scenario, the key quantity is the total thermal thermal conductivity which is required to be as low as possible. In Figure 7, the relative gains/losses on materials objectives of the two optimal compositions are shown for T 1000 K.
For the 5-PEs compositional space, three predominant objectives are the commodity cost, κ, and β during optimization. This is due to fact that these objectives deviate most from the reference system within the feasible compositional space. In Figure 4, changing the Fe concentration has almost no effects on material properties except the average commodity cost, while lower Co and Ni concentrations reduce κ at the expense of increment in β. Thus, the optimum candidate composition is predicted as Cr 0.2 Mn 0.2 Fe 0.3 Co 0.15 Ni 0.15 for which the relative gains/losses in objectives are shown in 7 (left). Since Co and Ni are the two most expansive elements among the original PEs, partially replacing them with Fe, which is the cheapest element, leads to a 15.1% reduction in the average commodity. It also reduces ρ by only 1% which leads to ∼ 90 kg mass reduction per 1 m 3 , implying further cost reduction due to lower density. Lower Co and Ni concentration reduces κ by 18.5% at the expense of a 6.7% increase in β, and sight worsened mechanical strength.
Overall gain may compensate these relatively small losses as long as they are still within safety ranges of the particular design. When Nb and Mo were introduced in the design, it was not feasible to reduce the commodity cost and ρ, as indicated by fittings in Supplementary Figures S21, S22. Furthermore, the cost and κ have generally inverse composition dependence in Figure 5, and it is challenging to simultaneously optimize these two predominant objectives within the feasible-solution region, determined by the phase stability and solubility constraints. Thus, the gain κ is relatively small to keep the cost lower. Nevertheless, Nb and Mo substitutions are quite promising to improve the mechanical and thermal properties. As shown in Figure 7 right panel, the composition Cr 0.1 Mn 0.2 Fe 0.2 Co 0.2 Ni 0.2 Nb 0.075 Mo 0.025 is promising as it provides significant improvements in material properties They simultaneously reduce κ and β, hinting better thermal insulation and less thermal expansion. There is also improvement of the lattice friction value σ y and K/G by 4.4% and 1.2%, respectively. Moreover, the SSH hardening of using simultaneous additions of Nb and Mo could lead to significant increments in the yield strength, as shown in Section 3.3. The 3 : 1 ratio of Nb:Mo substitution leads to 16.7 MPa SSH in the case of screw dislocations (α 4), and 80.6 MPa SSH in the case of edge dislocations (α 16). Supplementary Figure S23 shows the homologous optimized solutions obtained with the HOLMES interpolation scheme.

Effective Heat Sinks
For effective heat sinks, the main objective is to reach the highest κ for an optimal dissipation of excess heat. Two relative gains/ losses in two optimal compositions are shown in Figure 8 for T 1000 K. In this particular case, a lower Mn fraction is the key to achieve the main objective as it increases κ in Figure 4.
Thus, the optimal composition is predicted as Cr 0.25 Mn 0.05 Fe 0.3 Co 0.1 Ni 0.3 . It has a κ 33.4% higher compared to that of X 5 . Similar to the previous case, reducing Co while increasing Fe results in a considerable reduction ( − 33.5%) in the average commodity cost. Unlike the previous case, there is also considerable improvement in mechanical strength such as 21.6% higher σ y and 4.1% higher E at the expense of considerably lower Pugh ratio ( − 24.6%), indicating less ductility.
The 7-PEs optimal composition is predicted as Cr 0.2 Mn 0.1 Fe 0.2 Co 0.2 Ni 0.2 Nb 0.075 Mo 0.025 . It also relies on reducing the Mn concentration to increase κ for which a subtle 2.7% gains is achieved. Similar to the previous case, Nb and Mo substitution substantially increases the commodity cost and density. It also improves mechanical strength; however, not as much as in the case of Cr 0.25 Mn 0.05 Fe 0.3 Co 0.1 Ni 0.3 . On the other FIGURE 7 | Optimized candidate compositions for fiveand seven-principal-elements high-entropy alloys for thermally well-insulated components, operating at high temperatures. The Wendland radial basis functions are used for interpolating the compositional spaces.  Despite some improvements in the mechanical and thermal properties, it is application-dependent to understand whether they are sufficient enough to compensate the increased cost. Supplementary Figure  S24 shows the homologous optimized solution obtained with the HOLMES interpolation scheme.

DISCUSSION
This work explored the non-equimolar compositional space of five multi-principal elements (Co, Cr, Fe, Mn, and Ni) with and without additions of Mo and Nb using first-principles simulations. We found that the solubility FOMs and mechanical properties have commonly weak temperature dependence while thermal properties exhibited significant temperature dependence. This temperature dependence of G free is predominantly due to S conf and ΔF vib . Thus, while the configurational entropy is maximized when equimolar fractions are used, the vibrational contributions play a critical role in stabilizing HEAs. This observation is crucial as it opens up the possibility of stabilizing diverse compositions by contributing to the GFE.
Another important aspect of our first-principles data is that the changes in compositions go beyond affecting the average mass (m avr ), change in electronegativity (Δχ) and charge density (Δn v ), and misfit parameter (Δ(δr) which are widely used and computed using classical models. For instance, the phonon DOS critically changed compared to the equimolar case except for changes in Fe. These results, which were achieved due to the use of DFT calculations, point out the possibility to control phonondominated properties such as thermal conductivity, wave velocity, and phase transformation by tuning the HEAs composition . Indeed, the analyzed compositions, whose values are provided in the SI for the maximum reproducibility, indicated that thermal conductivity and the coefficient of thermal expansion could be drastically changed among all properties.
Of particular interest, we also explored non-equimolar compositions with additions of Mo and Nb. The trend of incorporating Mo and Nb is starting to gain momentum in the HEA community, with some experimental studies showing the synergistic effects of these two solutes (Fan et al., 2022;Rojas et al., 2022). Our results showed that small additions of both Mo and Nb could significantly improve the yield strength thanks to SSH. Indeed, we showed that a ratio of Nb:Mo of 1 : 3 is optimal to retain the electronegativity of the base alloy while providing significant SSH (i.e., ∼ 80 MPa). Remarkably, this ratio of Nb and Mo also maximizes the SSH value, as shown in Section 3.3 and Figure 6. An essential outcome of the SSH calculations is that the increment of the yield strength was weakly dependent on the temperature, and thus, can play a significant role in contributing to mechanical performance at high temperatures.
While the first-principles simulations and optimized compositions provide valuable insights, there are several limitations in our approach. For instance, the use of VCA is limited to HEAs that are in solid solution, and it is not possible to assess the thermodynamic stability of intermetallic phases that could appear in it. This limitation is more critical in HEAs combining Mo and Nb, where experimental works have shown that the volume fraction of σ, μ, and Laves phases increase with Mo and Nb content. For instance, Liu et al. (2016) and Shun et al. (2012) investigated CoCrFeNiMo x with x 0.1 − 0.8. Others have investigated the effect of Nb in non-equimolar HEAs, including Liu et al. (2015), and later He et al. (2016) and Jiang et al. (2017) studied the effect of Nb in the CoCrFeNiNb x system, with x 0.1 − 1.2. More recently, Sunkari et al. (2020) investigated a non-equimolar CoCrFeNi 2.1 Nb 0.2 HEA, and found excellent properties when it was cryo-rolled. These works revealed that these two solute atoms could be solved in the FCC matrix, although precipitation of intermetallic phases could appear for large molar fractions, that is, x ≥ 0.1. Thus, these observations suggest that controlled proportions of Mo and Nb should be employed. Remarkably, our study suggests that small proportions, that is, x 0,−, 0.1, could be employed to maximize SSH. Additionally, the proportion of Mo:Nb 3: 1 seems to be optimal in preserving electronegativity and maximization of SSH while being an economically viable option. This result is in agreement with the recent experimental work presented by Rojas et al. (2022).
In summary, we have systematically investigated nonequimolar CoCrFeMnNi with and without Mo and Nb additions. We have used thermalized first-principles simulations to construct a reasonably large data set. Based on this set, we have created a smooth high-dimensional interpolated space with LME and RDF to predict the highly non-linear response function of the HAEs properties as a function of composition and temperature. Combining the first-principles set with the interpolation scheme allowed us to find optimal compositions for two applications of interest related to optimal heat insulation and heat conduction using multi-objective metaheuristics optimization. The 5-PEs compositional space provides cost and density reductions with good phase stability and solubility thresholds and relatively small performance gains/losses in material properties. Thus, the nonequimolar 5-PE HEAs offer significant advantages when the cost is the most crucial factor when designing a component.

DATA AVAILABILITY STATEMENT
The first-principles materials data is available online in the repository Orhan, 2021. Further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
OO contributed to the original research topic, performed firstprinciples simulations, generated and analyzed materials database, and prepared the original manuscript. MI and LC performed multiobjective optimizations and contributed to the original manuscript. MP developed the original research topic, supervised, managed the funding acquisition, and revised the original manuscript.

FUNDING
We acknowledge the support from the Natural Sciences and Engineering Research Council of Canada (NSERC) through the