First-Principles Calculations of Elastic and Thermodynamic Properties for Multi-component Co-based Superalloys

First-principles calculations were performed to investigate the elastic and thermodynamic properties for multi-component Co-based superalloy systems and explored the effect of alloying on stabilizing the γ′ phase. First, the comparisons were carried out for the γ′ phase in Co3(Al,TM) (TM being transition metals) and Ni3Al systems between the present computational results using the EMTO-CPA method and other available DFT calculations as well as experimental data. The lattice parameters, elastic constants, and Debye temperatures are consistent with experimental results and other calculations. The predicted thermodynamic properties, e.g., the Gibbs free energy, excess entropy, and linear thermal expansion coefficient, agree well with CALPHAD results, experimental results, and other available first-principles calculations. A combination of EMTO-CPA method and Debye–Grüneisen model is utilized in this work to ensure that the alloying effect on the stability of the γ′ phase in a multi-component Co-based system is captured efficiently. This could open the path for designing novel multi-component Co-based alloys based on first-principles calculation. To demonstrate this, predictions for the properties of multicomponent systems were undertaken. Our results show that Ni aids in the stabilization of the (CoNi)3(Al, Mo, Nb) phase.


I. INTRODUCTION
CO-BASED superalloys have garnered a lot of attention since the possibility to form an ordered L1 2 -Co 3 (Al,W) (A 3 B type) phase was demonstrated [1] due to their high-temperature strength, environmental resistance, and cobalt's slightly higher melting point compared to nickel. However, the insufficient stability of the c¢ phase and their high density are still issues that need to be addressed. To achieve this, c¢ hardened Co-based superalloys with increasingly complex multi-component composition and lower density are developed. [2][3][4][5] Recently, Liu et al. [6] designed a multi-component Co-based superalloy based on machine learning results, namely, Co-36Ni-12Al-2Ti-4Ta-1W-2Cr with high c¢ solvus temperature (1266°C), low density (8.68 g/cm 3 ), and good high-temperature oxidation resistance. However, few studies of the thermodynamic and other properties of multi-component Co-based superalloys using first-principles calculations exist. Examples are works of Yao et al., [7,8] but the majority of such studies are experimental ones. Makineni et al. [9] optimized a new W-free Co-30Ni-10Al-5Mo-2Nb alloy with high specific strength (94.3 MPa g À1 cm À3 ) and lower density (8.38 g/cm 3 ) than Co-Al-W-based alloys by experiments. Due to its attractive mechanical properties, it may be a good candidate for high-temperature applications. It was found that Ta and Ti will additionally stabilize the c¢ phase. [10] However, such research if restricted to experimental methods will not be able to systematically scan the large amount of possible alloy compositions with reasonable effort and time consumption. First-principles calculations to accompany such alloy developments using modeling could accelerate such projects but are often hampered by the inability of such methods to tackle multi-component systems efficiently. [7,8,10] Fortunately, ab initio calculation schemes which are better suited to predict the properties of multicomponent systems are available now, such as the exact muffin-tin orbitals (EMTO) method [11][12][13] combined with the coherent potential approximation (CPA). [14,15] It was, for example, successfully applied to study the elastic parameters of partially ordered phases. [16] In this work, the purpose is to demonstrate the potential of the EMTO-CPA method by applying it to the development of multicomponent Co-based alloys. For the partly disordered L1 2 phases, we calculated the thermodynamic properties in order to estimate their thermodynamic stability at temperatures between 0 and 1400 K. To validate the approach, comparisons with available experimental [1,17,18] and other theoretical results [7,8,[18][19][20][21][22][23][24] using VASP and Thermo-Calc were made. The comparison with experimental and Thermo-Calc results can also be used to verify the applicability of the Debye-Gru¨neisen model as used in the present study to determine thermodynamic properties of the c¢ phases. The impact of the individual terms on the Gibbs free energy is briefly discussed. Based on accurate equations of state for the pure metals and L1 2 phases, current theoretical results provide a solid basis to predict the thermodynamic properties of multi-component Co-based systems with reasonable accuracy. The remaining sections of the paper are organized as follows. The computational details and methodologies for calculating thermo-physical properties are presented in Section II. The following Section III presents the calculated lattice parameters, enthalpy of formation, elastic properties, as well as the thermodynamic properties and compares them for selected cases with results yielded by other calculation methods or experiments. Finally, the main findings and conclusions are summarized in Section IV.

A. Theoretical Methods
DFT calculations were done with the exact muffin-tin orbitals method [11][12][13] (EMTO) combined with the coherent potential approximation (CPA) [14,15] for treating the chemical and magnetic disorder in solid solutions. The total energy is computed via the full charge-density technique. [25] Comparisons of the results obtained from EMTO-CPA calculations, with those from Vienna Ab Initio Simulation Package (VASP) calculations were performed. Here, the so-called special quasi-random structures (SQS) were used to model chemical disorder in VASP, which are generated by using the mcsqs code of the Alloy Theoretic Automated Toolkit (ATAT). [26] Finite-temperature properties were determined in the frame of the Debye-Gru¨neisen method in combination with the contribution of electronic free energy. The Gibbs free energy, taking excess entropy contributions (e.g., lattice vibration, thermal electronic excitation, and configurational entropy) into consideration, can be expressed as a function of Wigner-Seitz radius (r) and temperature (T) using the following approximations, see Eqs. [1] through [3], omitting the PV term being negligible at ambient pressure for condensed phases: [27] G r; Here E int is the internal energy, F el is the thermal electronic free energy, and the vibrational energy E D r; T ð Þ as well as vibrational entropy S D r; T ð Þ are calculated using Debye-Gru¨neisen model, which considers the equilibrium volume dependence of the Debye temperature H D with the Gru¨neisen constant. S conf: is the configurational entropy. Here, it also should be noted that the equilibrium volumes can be determined from fits to the curves of Gibbs free energy volume at each temperature. Equations [4] through [9] describe the relationship between elastic constants and Debye temperatures, where V is the equilibrium volume, and h and k B are the Planck and Boltzmann constants, respectively. v m is the average sound velocity with v L being the longitudinal and v T the transversal velocity. q is the density, B and G being polycrystalline bulk and shear modulus, respectively, obtained using Hill averaging. More details can be found in References 7,19,27. Therefore, thermal expansion and other thermal properties are treated within the quasi-harmonic approximation. All elastic parameters including the Debye temperature were calculated as a function of volume. In addition to the theoretical volumes determined by DFT calculation for selected cases for comparison also, available experimental room-temperature volumes were used. Such values are marked as expt. and vol. in the diagrams. As an alternative to the thermal expansion determined from fits to the Gibbs free energy, the thermal expansion coefficient of the whole compound was estimated as a weighted linear combination of the experimental thermal expansion coefficients of the contained pure metals. [28] By this, it is possible to calculate the volume at each temperature and the Gibbs free energy for the respective temperature can be calculated. Those values are marked as theo.vol. in the diagrams.
Besides, a semiempirical method, e.g., CALculation of PHAse Diagram (CALPHAD), is often applied to predict the phase diagram and quantify the thermodynamic properties, such as Gibbs free energy and excess entropy. In this work, all CALPHAD calculations were performed by using Thermo-Calc (version 2019a) software with the database of TCNI8.
For cubic crystals, there are three independent single-crystal elastic constants, C 11 , C 12 , and C 44 . C 11 and C 12 are derived from the bulk modulus B = (C 11 + 2 C 12 )/3 and the tetragonal shear modulus C¢ = ðC 11 À C 12 Þ=2. The theoretical equilibrium lattice parameter (a) and bulk modulus (B) were obtained from a Morse equation of state [19] fitting. C¢ and C 44 were determined by applying selected strains to the deformed unit cell inducing orthorhombic and monoclinic distortions, as described in Reference 19. Herein, six strains d = 0 to 0.05 pct (step = 0.01) were used. From these single-crystal elastic constants, the polycrystalline shear modulus (G) is obtained according to the Hill averaging, which is the arithmetic mean of the Voigt and Reuss average.

B. Computational Details
For EMTO calculations, the Green function was calculated for 16 complex energy points. The s, p, d, and f orbitals were included (l max = 3), and l cutoff l h max ¼ 8 was used for the full charge density. The one-electron Kohn-Sham equations were solved within the scalar-relativistic approximation and the soft-core scheme. The electrostatic correction to the single-site approximation was described using the screened impurity model with screening parameters of 0.6 and 0.9. [22] A 29 Â 29 Â 29 k-mesh was used for equation of state calculations, whereas 29 Â 29 Â 29 and 29 Â 29 Â 37 k-meshes were employed for the orthorhombic and monoclinic distortion calculations of the elastic constants, respectively. The self-consistent calculations were performed with the generalized gradient approximation (GGA) via the Perdew, Burke, and Ernzerholf (PBE) scheme. [29] The formation enthalpies of Co 3 (Al,TM) were also calculated using the Vienna Ab initio Simulation Package (VASP) [30,31] for comparing with EMTO calculations. The exchange-correlation functional was treated within the GGA formulated by PBE. [32] The semi-core p states were considered as valence shells in the pseudopotential for Nb, Mo, and Ta atoms (i.e., 4p 6 5s 1 4d 4 , 4p 6 5s 1 4d 5 , 5p 6 5d 3 6s 2 ), and the standard PAW potentials for Al, Ti, V, Cr, Co, Ni, and W elements (i.e., 3s 2 3p 1 , 3d 2 4s 2 , 3d 3 4s 2 , 4s 1 3d 5 , 3d 7 4s 2 , 3d 8 4s 2 and 5d 4 6s 2 ) were adopted, respectively. To consider chemical disorder, all of the supercells employed in VASP were generated via SQS method with 32 atoms, and we also carefully examined the convergence of energy differences in supercells with 32 and 64 atoms. All the atoms are relaxed including all cell-internal atomic positions and unit cell volume by using a conjugate gradient algorithm. The calculations were done with an energy cutoff of 600 eV and C-centered k-point meshes were constructed so that the k-point distance in all three directions was 2p Â 0:03 A À1 (i.e., k-meshes are 5 Â 5 Â 5 for 32 atoms SQS supercell). Structural relaxation was done with the Methfessel-Paxton method with a smearing width of 0.1 eV, and then the tetrahedron method with Blo¨chl corrections was used for accurate total energy during a final static calculation. The total energy tolerance and ionic force tolerance to reach convergence were set as 10 -5 eV/atom and 10 -3 eV/Å , respectively.

III. RESULTS AND DISCUSSION
A. Lattice Parameter Figure 1 presents the calculated lattice parameters for pure metals and L1 2 -Co 3 (Al,TM) (TM being transition metals). Available DFT [19,20] and experimental [33] results are also included for comparison. For pure metals, our EMTO results show an excellent agreement with previous EMTO and VASP values, as well as experimental data. Before considering the Co 3 (Al,TM) alloy system, we computed and compared the lattice parameter for the Ni 3 Al system, as shown in Figure 1(b). Its lattice parameter agrees well with the experimental one. [17] Furthermore, we also estimated the Curie temperature (T C ) for Ni 3 Al using the mean-field approximation, viz., where c is the concentration of nonmagnetic elements, k B is the Boltzmann constant, and E PM À E FM is the difference of total energy at the paramagnetic (PM) and ferromagnetic (FM) states. The calculated Curie temperature for Ni 3 Al is 29 K and is in good agreement with the experimental value (41 K), [34] which indicates that the Fig. 1-Theoretical and experimental lattice parameters for pure metals and the partial-disordered L1 2 phase alloys (Co 3 (Al,TM) and Ni 3 Al). (a) Comparison between the present EMTO calculations and DFT results from the literature (EMTO [19] and VASP [20] ) as well as available experimental data at room temperature (RT). [33] (b) The calculated lattice parameters for Co 3 (Al,TM) and Ni 3 Al alloys by using EMTO-CPA with two screening parameters (0.6 and 0.9). Two existing experimental values for Co 3 (Al,W) [1] and Ni 3 Al [17] are included for comparison. present EMTO method is able to capture this thermodynamic property of a L1 2 alloy system reliably.
With respect to the partial-disordered L1 2 phase, the present EMTO calculations for Co 3 (Al,TM) alloys were performed within two screening parameters (ALPCPA = 0.6 and ALPCPA = 0.9). The screening parameter represents the effect of the charge misfit on the spherical potential using the screened impurity model (SIM) [22,35] and can be used to match the charge transfer obtained within the CPA and the real charge transfer, controlling the effectiveness of screening around the impurity. [19] As presented in Figure 1(b), we noticed that the choice of the screening parameter can make a significant difference for some of the transition metals alloyed to the Co 3 (Al,TM) system. The screening parameter has no effect on the lattice parameter of L1 2 compounds with 3d elements (Ti, V, and Cr), while showing a significant impact on L1 2 compounds containing 4d elements (Nb, Mo) and even more pronounced if 5d elements (Ta and W) are involved. Although the deviation from the unit slope for Co 3 (Al,TM) alloys is slightly larger than that for the pure metals, still the EMTO lattice parameters for L1 2 phases are in good agreement with those calculated by VASP. Overall, the lattice parameters calculated by EMTO are slightly higher than those calculated by VASP. This does of course not imply that the lattice parameters determined by EMTO are less precise than those determined by VASP but just that the two do not fully coincide. It is noteworthy that the EMTO results are closer to the available experimental values than VASP results as, for example, for Co 3 (Al,W). For the Co 3 (Al,W) phase also, a better agreement between EMTO and experimental results is found for a screening parameter of 0.6 instead of 0.9.

B. Enthalpy of Formation
To evaluate the chemical stability of the L1 2 phases and the difference between the EMTO-CPA and VASP-SQS methods, enthalpies of formation were determined with both methods. The enthalpy of formation and the corresponding relative enthalpy of formation for the L1 2 phases are plotted in Figure 2. The VASP results with and without relaxation are also plotted here for comparison: the corresponding labels are VASP-SQSr and VASP-SQSu, respectively. Due to the effect of the screening parameter on lattice parameters, two screening parameters within the EMTO-CPA scheme were adopted to compute their formation enthalpies for L1 2 -Co 3 (Al,TM) and Ni 3 Al alloys, which are shown in Figure 2(a). The red color represents the enthalpy of formation with screening parameter 0.9, and the green color represents the enthalpy of formation with screening parameter 0.6. The enthalpies of formation of Ni 3 Al calculated by EMTO with both screening parameters and VASP-SQSr fit well with each other, as well as the calculation reported in Reference 21. For the Co 3 (Al,TM) system, the screening parameters have a significant impact on the enthalpy of formation for the partial-disordered L1 2 phases, and the trend is similar to the effect on lattice parameters. As previously noted, when the third metal belongs to the 4d and 5d group, especially the results using a screening parameter of 0.6 show some scatter and a significant difference between EMTO and VASP results, as presented in Figure 2(a). The red data with screening parameter 0.9 are located along the unit slope, also for third elements belonging to the 4d and 5d group, but are systematically higher than the VASP results. Due to those results, the following calculations are all performed with screening parameter 0.9. The fact that the EMTO results are slightly higher than VASP results could stem from the fact that both represent two different methods of DFT calculations. Another possible factor is the local relaxation owing to the fact of large misfits between transition elements alloyed into the Co 3 (Al,TM) lattice. To validate the latter assumption, the formation enthalpies were calculated by VASP with relaxation (VASP-SQSr) and without relaxation (VASP-SQSu). For the sake of comparability, here we present the values of Co 3 (Al,TM) relative to the formation enthalpy of Co 3 (Al,W) in Figure 2(b). A set of EMTO-CPA results using screening parameter of 0.9 is included for comparison. All three methods exhibit the same trend of relative enthalpy of formation. It is noteworthy that the results of EMTO-CPA differ from both sets of VASP results. However, VASP results without (VASP-SQSu) and with relaxation (VASP-SQSr) are in relatively good agreement with each other. This implies that the energy change caused by localized lattice relaxation is fairly small and that the effect on the enthalpy of formation can be ignored. Accordingly, the differences between EMTO-CPA and VASP-SQS seem to stem from the employment of two different DFT methods. Nevertheless, the results in Figure 2(b) demonstrate that both methods depict the same relative trends in energy and enthalpy differences.

C. Elastic Properties
Starting from unit cells with theoretical equilibrium lattice parameters, the single and polycrystalline elastic moduli were calculated at 0 K using the EMTO-CPA method for pure metals (Table I), Ni 3 Al, Co 3 (Al,W) (Table II), and a multicomponent L1 2 phase of (Co-xNi) 3 (Al,Mo,Nb) type (Table III). For comparison, available numerical data from VASP-SQS, [21,24,36,37] WIEN2K, [24] and experimental data [17,18,23] are also Ref. [18]. c Ref. [36]. d Ref. [23]. e Ref. [24]. f Ref. [21]. g Ref. [37]. h Ref. [17]. listed in the appropriate tables. The polycrystalline elastic moduli were derived using Hill's averages, with the exception of a few values taken directly from the references. Figure 3(a) presents the single-crystal elastic constants for Ni 3 Al calculated using EMTO, compared with VASP calculation results [36] as well as the available experimental data. [18] The results derived with EMTO are in good agreement with the experimental data, despite for C 44 and to a lesser degree C 12 , which slightly differ from the experimental values. [18] While the VASP results [36] resemble C 12 , C 44 , and B well, the C 11 determined by VASP is away from the experimental values. This means that DFT calculations at 0 K irrespective of the used program package and method cannot completely capture the experimental data at 300 K. When comparing the polycrystalline elastic constants between EMTO, VASP, and experimental results, as can be done based on Table II, the EMTO   results show an excellent agreement with the experimental results for the calculated B/G ratio, Poisson ratio (m), and Cauchy pressure (C 12 À C 44 ), as well as the Debye temperature. Figure 3(b) presents the elastic constants for Co 3 (Al,W) including ones calculated by EMTO and other methods [21,24] and corresponding experimental results. [23] Here, the EMTO calculations were also performed using two screening parameters (0.6 and 0.9). As shown in Figure 3(b), it can be concluded that different screening parameters have almost no effect on the calculated elastic constants. The single-crystal elastic constants C 11 and C 44 , with screening parameter 0.9 slightly differ from the experimental values, the positive deviations being 11.3 and 23.4 pct, respectively. When using the experimental lattice parameter to calculate the single-crystal elastic constants with screening parameter 0.9, the deviations decrease to 4.08 and 16.95 pct (Table II). In addition to the effect of the cell volume or  [36] and experimental data [18] ; (b) Comparison of the present elastic constants for Co 3 (Al,W) with the previous VASP results, [21] WIEN2K, [24] and experimental data [23] .
lattice parameter, not considering relaxation of geometry and ionic positions could also be responsible for the deviation between EMTO and experimental results. Based on the calculated elastic properties in Table II, the present EMTO results in general show good agreements with the experimental data, except for the Cauchy pressure. The Cauchy pressure indicates the angular character of the bonding orbital [40][41][42] which seems to be not properly depicted in the present case by the DFT calculations. This is evidenced by the fact that not only EMTO but also other calculations (VASP [21] and WIEN2K [23] ), presented in Table II, predict a negative Cauchy pressure. It should be noted that the relaxation of ion positions can be ruled out as a reason for the present negative Cauchy pressure, as for the VASP and WIEN2K results shown, the former one performed with relaxation while the latter one performed without relaxation, both result in negative Cauchy pressure. A possible reason for the difference could be the type of exchange-correlation function used. In Reference 21, this was the generalized gradient approximation (GGA) with the Perdew and Wang parameterization (PW91), while in Reference 43, the GGA-PBE was employed. Both of them used the same special quasi-random structure (SQS), but the latter one predicted positive Cauchy pressure (C 12 À C 44 ). In addition, the sign of the Cauchy pressure seems not to be influenced by the adopted lattice constant in the calculation because using the experimental lattice parameter yields a negative Cauchy pressure as well. When using the experimental lattice constant in EMTO calculations, the deviations for C 11 , C 44 , B/G ratio, Poisson ratio (m), and Debye temperature are 4.8, 16.98, À 13.7, À 10, and 7.8 pct, respectively, compared to the results obtained using the equilibrium lattice volume. With exception of the sign of the Cauchy pressure with the present DFT calculations of elastic constants of Co 3 (Al,W), as evidenced by VASP [21] and WIEN2K [24] results, reasonable results for elastic constants are acquired with EMTO that agreed also with the experimental values.
As described above, the EMTO-CPA method can describe the elastic properties for pure metals and different L1 2 compounds with reasonable accuracy. As the number of multi-component L1 2 -strengthened Co-based superalloys currently developed grows, the L1 2 phase becomes more complex. [2,3,5,9,44] A huge advantage of the EMTO-CPA method is a possibility to calculate properties of multi-component compounds with relatively low computational effort. To demonstrate this the effect of alloying elements on elastic moduli for a multi-component Co-based system as a function of Ni concentration was calculated. The L1 2 phase is based on the composition of the c¢ phase of the alloy Co-xNi-10Al-5Mo-2Nb as reported in Reference 9: Co-(0-37.4)Ni-10.8Al-8.1Mo-3.7Nb with Co and Ni on A site and Al, Mo and Nb on the B site. For sake of brevity, it will be called L1 2 -(Co-xNi) 3 (Al,Mo,Nb) in the following.
The calculated elastic properties are presented in Figure 4, and all numerical data are listed in Table III. In Figure 4(a), it can be seen that all elastic constants (B, C 11, C 12, and C 44 ) first increase with increasing Ni concentration, until x(Ni) reaches 15 at. pct. Subsequently, the bulk modulus (B) and C 12 and C 44 slightly decrease. On the contrary, C 11 slightly increases with increasing Ni concentration. The polycrystalline elastic properties for c¢: Shear Modulus (G), Young's modulus (E), Pugh ratio (B/G), and Cauchy pressure (C 12 À C 44 ), are presented in Figure 4(b). G, E, and Debye temperature (H D ) slightly increase with increasing Ni concentration which are also observable in Table III. The Pugh ratio B/G follows the same trend as B, increasing first with increasing Ni concentration and then decreasing slightly. According to previous studies on the relation between the brittle/ductile behavior of materials and the Pugh ratio B/G, [42] Niu et al. [45] reported that if B/G is in the range of 1.67 to 1.754, the nature of materials is   5-Calculated entropies, Gibbs free energies, and thermal expansion coefficients as a function of temperature for Ni 3 Al (a) to (c) and Co 3 (Al,W) (e) to (f). Available VASP [8,21] and experimental data [48] are included for comparison.
brittle, whereas all values of B/G > 1.754, as found for the (Co-xNi) 3 (Al,Mo,Nb) compounds in our calculations (Figure 4(b) and Table III), indicate a ductile behavior. With increasing Ni concentration also,

B/G increases and in parallel should do the ductility.
Here, it should be noted that the Pugh ratio might provide a general indication of the behavior for materials with the same or similar crystal structure, but the validity is limited and predictions of the behavior of novel materials should only be made with caution. [46] In addition, it is beyond the capabilities of current DFT calculations to allow direct conclusions about the resistance of a material against plastic deformation. Nevertheless, comparing the effects of increased Ni concentration to improve the mechanical properties of the alloys Co-xNi-10Al-5Mo-2Nb as reported in Reference 9, with calculation results provides some insights. Experimentally the best mechanical properties were found with a 30 pct addition of Ni in Co-30Ni-10Al-5Mo-2Nb alloy. This corresponds to the values calculated for x(Ni) of 50 pct in Table III (Co-37.4Ni-10.8Al-8.1Mo-3.7Nb is composition of the c¢-phase measured by experiment), which is also the maximum Ni concentration experimentally tested in Reference 9. It is noteworthy that the values of Cauchy pressure (C 12 À C 44 ) stay positive and follow the same trend as Pugh ratio B/G with increasing Ni content. This would indicate that the bonding character of (Co-xNi) 3 (Al,-Mo,Nb) is more metallic bonding like with increasing Ni concentration, resulting in a better ductility as evidenced by the B/G ratio.

D. Thermodynamic Properties
Here, three terms of entropy, e.g., vibrational (S vib. ), electronic (S el. ), and configurational (S conf. ) contributions, are considered for the L1 2 -type Co-based alloy system. The vibrational entropy contributions are calculated using the Debye-Gru¨neisen approach, with the volume change as a function of temperature [27] introduced via the quasi-harmonic model. [47] The thermal electronic and configurational contributions are included for determining the effect of each one of them on the thermodynamic properties. The calculated thermodynamic properties, i.e., Gibbs free energy, entropy, and thermal expansion coefficient, as a function of temperature for Ni 3 Al and Co 3 (Al,W) alloys are presented in Figure 5. For comparison, available VASP [7,8,21] and experimental data [48] as well as the CALPHAD results derived from Thermo-Calc software, are also plotted.
For the Ni 3 Al phase, the temperature-dependent total entropies ( Figure 5(a)) and Gibbs free energy ( Figure 5(b)), are in excellent agreement with the Thermo-Calc results, especially, when taking the electronic contributions into consideration (e.g., see green symbols). The calculated linear thermal expansion coefficients with and without the electronic contributions for Ni 3 Al are plotted in Figure 5(c). The available experimental data [48] and the calculations by VASP [8] are also plotted for comparison. A conclusion can be drawn that current calculations are in good agreement with the experimental data. The only exception is the linear thermal expansion coefficient which deviates from the experimental results above the Debye temperature (~480 K) as presented in Figure 5(c) and Table II. This means that the present method is no longer available once the service temperature exceeds the Debye temperature. In addition, the major contribution to the different thermodynamic properties stems from the vibrational energy and entropy. The agreement, especially with Thermo-Calc results, is nevertheless improved if the electronic excess contribution is included.
The calculated total entropies, Gibbs free energies, and the linear thermal expansion coefficients at different temperature for Co 3 (Al,W) are presented in Figure 5(d) through (f), respectively. The agreement with the Thermo-Calc results and the results calculated by EMTO-CPA (pink symbols) and VASP (cyan line) using Debye-Gru¨neisen method, is still good especially taking into account the excess contributions (S el. and S conf. ). It is possible that the remaining difference between the EMTO and Thermo-Calc results stems from the magnetic entropy contribution which was not considered in the calculations of thermodynamic properties based on the EMTO-CPA results. With respect to the significant deviation between the present results by EMTO-CPA as well as Thermo-Calc and the VASP results [21] shown in Figures 5(d) and (e), the reason is not fully clear. The phonon frequencies and thermodynamic properties are calculated from the force constants using the PHONOPY package in Reference 21 but this should, due to depicting the real phonon spectrum and not only an assumption of it as in the Debye-Gru¨neisen method, yield better but not worse results.
Based on the predictions in Ni 3 Al and Co 3 (Al,W) above, the Debye-Gru¨neisen approach was also applied to the multi-component Co-based system, and the composition of the c¢ phase are Co-10.6Al-8.7-Mo-4.8Nb and Co-37.4Ni-10.8Al-8.1Mo-3.7Nb, respectively. For sake of brevity, the above two representative c¢ phases are annotated as (Co-xNi) 3 (Al, Mo,Nb) with x = 0 and 50. The calculated total entropies, Gibbs free energies, and thermal expansion coefficients are shown in Figure 6. Note that the nominal composition for Co 3 (Al,Mo,Nb) is slightly different for the calculations done with VASP and EMTO. Since CPA method is able to exactly describe the experimental composition in EMTO scheme, while the VASP method sometimes is limited in its ability to model a complex composition within a reasonably small supercell. Here, Co-9.375Al-9.375Mo-6.25Nb is the composition of the compound used for VASP calculations, while the one for EMTO being exactly the same as in the experimental alloy. [9] The entropy results are in good agreement between the EMTO-CPA and the available VASP result, [7] as can be seen in Figure 6(a). The entropy calculated by Thermo-Calc is systematically higher and accordingly the Gibbs free energy systematically lower compared to the VASP as well as the EMTO-CPA results. The same trend is found for (Co-50Ni) 3 (Al,Mo,Nb) ( Figure 6(d)), when comparing EMTO-CPA and Thermo-Calc results perhaps even to a slightly stronger extent as shown in Figure 6(d). That no similar VASP results for (Co-50Ni) 3 (Al,Mo,Nb) are available in the literature, which probably stems from the necessary large supercell size to depict increasingly complex chemical compositions. It is most probable that the differences between EMTO-CPA and Thermo-Calc for the (Co-xNi) 3 (Al,Mo,Nb) stem from the non-consideration of different excess contributions (for instance, magnetic contribution, total anharmonic contribution, and longitudinal spin fluctuation contribution [49] ).
Nevertheless, by considering the contribution of configurational entropy the match between Thermo-Calc and EMTO-CPA results for the Gibbs free energies can be significantly improved. Especially for the complex multi-component (Co-xNi) 3 (Al,Mo,Nb) alloy, where a higher degree of mixing on the sub-lattices existed. The linear thermal expansion coefficients are presented in Figures 6(c) and (f). With increasing Ni addition, the linear thermal expansion coefficient slightly decreases. Due to the slight difference of nominal compositions for Co 3 (Al,Mo,Nb) between the available VASP and EMTO-CPA results, the difference in linear thermal expansion coefficients between the two methods is explicable (Figure 6(c)). A conclusion can be drawn that the current EMTO-CPA calculations are reasonable, accurate, and as exact as VASP results (if the system is treatable at all with a supercell approach) to predict thermodynamic properties of multi-component Co-based systems. Thus, they are useable with state-of-the-art accuracy to design novel multi-component Co-based alloys based on first principles. It is striking that an increasing c¢ solvus temperature with increasing Ni content is reported in Reference 9 while the difference in the Gibbs free energies calculated by EMTO-CPA but also by Thermo-Calc does not support this finding. This will be discussed more deeply in the following section.
E. Stability of c¢ Phases Figure 7 depicts the formation energies for different L1 2 phases and the influence of Ni addition on the lattice constant and formation energy of the multi-component (Co-xNi) 3 (Al,Mo,Nb) phase. As shown in Figure 7(a), the formation energy for (Co-xNi) 3 (Al, Mo,Nb) decreases with increasing Ni concentration. This is evidence that the higher c¢ solvus temperature found for increasing Ni concentration in alloys with this phase [9] does not stem from different thermodynamic properties but can be attributed to the lower energy of formation of the L1 2 compounds richer in Ni. At the same time, the lattice parameters slightly increase with increasing Ni concentration. The lattice parameters are in good agreement with results from other calculations as well as available experimental results [8,50] also shown in Figure 7(a). The difference of lattice parameter between the EMTO-CPA calculations and experimental lattice parameter is about 0.4 pct. Furthermore, compared to other Co-based L1 2 phases, the formation energy of (Co,Ni) 3 (Al,Mo,Nb) is the lowest, which indicates that the L1 2 -(Co,Ni) 3 (Al,Mo,Nb) is more stable than the other Co-based L1 2 phases, as shown in Figure 7(b). The tendency of formation energies (FE) comparing the four phases agree with similar results in References 7, 8. It is noteworthy that the addition of Ni has a significant effect on stabilizing the multicomponent L1 2 phase. The change in formation energy between the two multicomponent L1 2 phases when adding Ni is larger than the difference between L1 2 -Co 3 (Al,Mo,Nb) and L1 2 -Co 3 (Al,W).

IV. CONCLUSIONS
In this work, the exact muffin-tin orbitals method (EMTO), combined with the coherent potential approximation (CPA) based on density functional theory (DFT), was used to investigate the elastic and thermodynamic properties, as well as the effect of alloying on stabilizing the c¢ phase in multi-component Co-based systems. The EMTO-CPA calculations were compared and validated against other available DFT, CALPHAD, and experimental results. The reasonable agreement in most cases with CALPHAD results, experimental results, and other available calculations demonstrates that the aforementioned properties can be reliably predicted for L1 2 -hardened Co-based alloys using the EMTO-CPA approach with at least the same accuracy as with a VASP supercell approach. Based on the results, the main conclusions can be summarized as follows: (1) The lattice parameter, elastic constants, and Debye temperature are consistent with experimental results and other calculations for the pure elements and Co 3 (Al, TM) (TM being transition metals) and Ni 3 Al systems.