Modeling of vibrational and configurational degrees of freedom in hexagonal and cubic tungsten carbide at high temperatures

Transition metal carbide is a class of materials characterized by high hardness, high melting points, and low chemical reactivity. It is widely used in industrial applications involving exposure to elevated temperatures, aggressive media, and heavy loads, and is thus of technological and industrial importance. In this paper the hightemperature thermodynamic properties of tungsten carbide, WC, is studied. At most temperatures below melting, WC assumes a hexagonal structure with essentially no vacancies (δ-WC). Only at very high temperatures (around 3000 K), WC crystallizes in the cubic rocksalt structure (γ -WC), which is more common for the transition metal carbides and in the case for WC can contain up to 40% carbon vacancies. At lower temperatures, γ -WC can, however, form as thin interfacial structures or nanoparticles. Hence, the thermodynamic properties of both γ -WC and δ-WC are of relevance. Here, we conduct a first-principles density-functional theory based computational investigation of the γ -WC and δ-WC phases, which requires modeling of high carbon vacancy concentrations at high temperatures. The configurational degrees of freedom are modeled with an alloy cluster expansion and sampled through Monte Carlo simulations. To account for the dynamic instability of the cubic γ -WC phase at low temperatures, the vibrational degrees of freedom are treated using effective harmonic models constructed from ab initio molecular dynamics simulations. Finally, we obtain a part of the W-C phase diagram in reasonably quantitative agreement with experimental data.


I. INTRODUCTION
Transition metal carbides are characterized by high hardness, high melting points, and low chemical reactivity. These carbides are, therefore, widely used in manufacturing as structural and tool materials capable of operating at elevated temperatures, in aggressive media, and under heavy loads. Hence, transition metal carbides are of great technological and industrial importance [1][2][3][4][5][6][7][8][9][10][11][12][13]. The group IV (Ti, Zr, Hf) and group V (V, Nb, Ta) carbides crystallize in the cubic rocksalt structure with a large homogeneity range, maintained by the possibility of having a large amount of carbon vacancies (up to 50%) in the material [1][2][3]. Group VI (Cr,Mo,W) carbides, on the other hand, form different cubic and hexagonal structures.
At room temperature, tungsten carbide (WC) only assumes a hexagonal structure, referred to as δ-WC, with essentially no vacancies present, i.e., it has an insignificant homogeneity region [1,[14][15][16]. Above 1523 K also a semicarbide W 2 C phase becomes stable, having several structural modifications, all with rather broad homogeneity ranges [16,17]. At high temperatures, above 2789 K, a cubic tungsten carbide phase γ -WC with rocksalt structure becomes stable in the composition range from WC 0.6 to stoichiometric WC and can contain a lot of vacancies [16]. This is the same phase as formed in the group IV (Ti, Zr, Hf) and group V (V, Nb, Ta) carbides. The cubic tungsten carbide γ -WC starts to melt at T = 3028 K [16].
Due to its favorable mechanical and thermal properties tungsten carbide (WC) is extensively used in wear-resistance hard alloys. If WC powder is sintered together with a metal Co powder, a composite material, cemented carbide, is formed that exhibits an excellent combination of hardness, toughness, and wear resistance [18]. WC-Co is therefore of high technological importance and is extensively used for various metal machining operations such as turning, milling, drilling, threading, grooving, etc. [19,20]. By doping the material with various transition metal (M) atoms one has shown that the growth of the WC grains in the Co-liquid phase sintering process can be retarded, leading to a fine-grained material with improved hardness [21,22]. The growth inhibition is believed to be related to the formation of ultrathin MC films, comprising about two atomic layers, that are formed at phase boundaries between WC grains and the Co binder phase [23][24][25][26][27][28][29][30][31].
Also in undoped WC-Co, ultrathin layers of cubic WC (γ -WC) have been observed at low temperatures [31][32][33][34]. This is interesting since the bulk γ -WC phase is stable only at very high temperatures, above 2789 K. The γ -WC phase has also been identified as the primary WC phase in ultrafine WC particles [35]. These observations indicate that the γ -WC phase is most likely stabilized by interfacial effects. In the case of interfacial layers these structures are referred to as a complexions [36,37].
It is believed that the cubic interfacial layers are formed at the sintering temperature of about 1700 K, both in the undoped and the doped cemented carbides. They are considered to be thermodynamically stable structures and to predict their occurrence, composition, and structure, as function of temperature and doping condition, is a challenge. Due to their technological relevance and the lack of experimental information, it is highly desirable to develop an accurate computational approach to derive phase diagrams describing structure and composition of the interfacial layers as function of external conditions, so-called interfacial phase diagrams [38].
This paper is a first step in obtaining accurate interfacial phase diagrams for cemented carbides. The aim with this paper is to derive a part of the bulk phase diagram for the binary W-C system describing the stability region of the hexagonal δ-WC and the cubic γ -WC phases using first-principles techniques. For the bulk case an experimental phase diagram is available and we can directly make a comparison and obtain a good quality assessment of our method. This study will be followed by studies of interfacial phase diagrams for both undoped [39] and doped [40] cemented carbides.
Generating the W-C phase diagram requires the calculation of the full temperature-and composition-dependent free-energy landscapes, including effects from vibrations and various configurations, for the δ-WC and γ -WC phases at high temperatures. Modeling vibrations and configurations for a wide composition range and at high temperatures is crucial to obtain the correct free energies [41]. Further, the stoichiometric γ -WC phase is dynamically unstable at 0 K, rendering standard vibrational free-energy techniques, such as quasiharmonic approximation (QHA), invalid. It is therefore important to couple vibrational and configurational degrees of freedoms (DOFs), a challenging task [42].
Our modeling is based on density functional theory (DFT) for the electronic-structure with exchange-correlation contributions described by the Perdew-Burke-Ernzerhof (PBE) functional [43] as implemented in VASP [44,45]. More details concerning the electronic structure calculations can be found in Appendix B. The vibrational contributions to the free energy are treated using effective harmonic models (EHMs), where the force constants are fitted to snapshots from ab initio molecular dynamics (AIMD) simulations using the HIPHIVE package [46]. The free energies are then evaluated using the PHONOPY [47] program and are further improved using free-energy perturbation (FEP) [48] theory. The configurational part, which arises from carbon vacancies, is treated by constructing an alloy cluster expansion (CE) using the ICET package [49] after which the CE is used to perform thermodynamic integration based on results from Monte Carlo (MC) simulations. We also incorporate electronic free energy and thermal expansion in the modeling. Finally, we derive a part of the binary W-C phase diagram from first principles, obtaining the stability region for the δ-WC (space group P6m2) and γ -WC (space group F m3m) phases.

A. Phase diagrams
A phase diagram shows the stable phases α given the temperature T , pressure P, and overall composition c = N C /(N C + N W ), where N C and N W are the total number of C and W atoms, respectively. The stable phases are found by minimizing the free energy of the system given these macroscopic thermodynamic parameters, i.e., minimizing the weighted sum of the individual free energies G α (c α , T, P), with the constraint of preserved overall composition Here, w α denotes the fraction of phase α in the system and c α denotes the composition of phase α.
In this study we consider only the zero-pressure case. The Gibbs free energy G then equals the Helmholtz free energy F . To construct the phase diagram, one first has to construct a full temperature-and composition-dependent free-energy landscape for all relevant phases, i.e., derive models and compute results for G α (c α , T ) = G α (c α , T, P = 0), and then minimize the weighted sum according to Eq. (1).

B. Free energy
The δ-WC and γ -WC phases consist each of two sublattices of size N mix , containing either only C or W lattice points. In this study we only consider vacancies in the C sublattice, hence, N W = N mix is fixed while N C N mix . The composition y of the system is uniquely defined by the C concentration relative to the C sublattice and may be expressed as the C to W ratio, i.e., y = N C /N W = N C /N mix .
We start by considering the Helmholtz free energy per atom of phase α, i.e., the total free energy divided by (N C + N W ), and denote it as F α . It depends on the composition y, specific volume v, and temperature T , i.e., F α = F α (y, v, T ). We assume that it can be decomposed into four parts where E 0 α (y, v) is the relaxed total energy from an electronic structure calculation and F el α (y, v, T ), F vib α (y, v, T ), and F conf α (y, T ) are the temperature-dependent free-energy contributions from excitations of electrons, motion of atoms, and atomic configurations, respectively. Furthermore, the free energy relevant for the phase diagram is at zero pressure. We, therefore, minimize F α (y, v, T ) with respect to v, thus incorporating thermal expansion, to get Gibbs free energy per atom at zero pressure according to When the volume is varied, the cell shape of the structures (simulation cells) is kept fixed. This means that the c/a ratio for the hexagonal δ-WC structure is kept fixed for all temperatures. 033804-2

C. Electronic free energy
The electronic contribution to the free energy F el α (y, v, T ) arises from the excitation of the electrons from their ground states. The electronic entropy S el (T ) can be computed from the electronic density of states (EDOS) as [50] Here, n el (ε) is the EDOS and f (ε) the Fermi distribution where β = 1/k B T and μ el (T ) is the electron chemical potential at temperature T . μ el (T ) is computed from the condition of conserved number of valence electrons N el val , i.e., Finally, by using the thermodynamic relation

Sommerfeld approximation
At not too high temperatures and assuming the freeelectron gas model, the free energy can be approximated by the Sommerfeld expression [51] F el (T ) = − π 2 6 n el (ε F )k 2 B T 2 , where ε F = μ el (T = 0) is the Fermi energy, i.e., the electron chemical potential at T = 0 K.

D. Vibrational free energy
The vibrational contribution to the free energy F vib α (y, v, T ) arises from the motion of the atoms and, since the considered phases are solid at the temperatures of interest, a harmonic approximation is a suitable starting point for analyzing the vibrational free energy.

Harmonic approximation
In the harmonic approximation (HA) the potential energy is written as where i j are the harmonic force constants and u i the displacement of atom number i. The vibrational density of states (VDOS) n ph (w) can be calculated from i j and allows the harmonic free energy to be calculated as [52,53] where is the free-energy contribution of a mode with frequency ω at temperature T .

Effective harmonic models
At 0 K the γ -WC phase is dynamically unstable, rendering the harmonic approximation invalid. Further, since the interesting temperatures for the δ-WC to γ -WC phase transition are close to the melting point, it is likely that anharmonicity plays an important role in both phases.
By using EHMs, both dynamic instabilities and anharmonicity in a system can be approximately taken into account [54,55]. In this work, we use temperature-dependent EHMs which are trained using the HIPHIVE package [46] to forces from AIMD simulations. These effective harmonic force constants effectively incorporate anharmonic effects and introduce a temperature dependency = (T ). It follows that also the frequencies have a temperature dependency, i.e., ω i = ω i (T ). The VDOS from the EHMs can be used in Eq. (11) in order to obtain the free energy F EHM (T ).
The EHM is only accurate for temperatures close to the temperature of the AIMD simulation. Hence, multiple EHMs need to be constructed at several different temperatures and the resulting free energies, which are calculated using Eq. (11), need to be interpolated to get a good description across a larger temperature range.

Free-energy perturbation
To improve the accuracy of the free energies obtained using EHMs and Eq. (11), free-energy perturbation (FEP) [48] is used, according to which the free-energy difference between EHM and DFT is given by where U denotes the internal energy and . . . denotes a thermal average. This correction to the EHM free energy yields a significant improvement in accuracy [56,57], and the average can readily be evaluated using DFT since AIMD simulations are used to construct the EHMs in the first place. This correction is often significant for unstable systems as the relaxed internal energy, i.e., E 0 , will not coincide with an extrapolation of temperature-dependent internal energies, i.e., E (T ), down to 0 K [57]. Instead, the extrapolated value will lie lower compared with E 0 , which is picked up in FEP. The total vibrational free-energy contribution is given by

E. Configurational free energy
For solids the configurational contribution to the free energy F conf α (y, T ) arises from changes in the occupation of the lattice sites. The γ -WC phase is known to have up to 40% carbon vacancies [16]. Here, we model the ordering energetics of the carbon atoms in γ -WC by constructing an alloy CE for 033804-3 the total electronic structure energy [58]. The CE model can be viewed as a generalized Ising model and the energy of a configuration is expressed as where σ is the occupation vector where each value represents the occupation of a lattice site (here, −1 and +1 for sites that are occupied by a C atom or are vacant, respectively). The summation in Eq. (14) runs over all orbits a (also referred to as symmetry-inequivalent clusters). Their respective multiplicities are denoted by m a and a (σ) a represents the average occupation of orbit a, obtained by averaging over all clusters within that orbit (i.e., its star).
J a are the effective cluster interactions (ECIs) of the model. They are obtained by training against DFT computed relaxed structures using the ICET package [49]. Atomic relaxation and strain contributions to the energy are thereby effectively incorporated in the ECIs, even though the CE itself is restricted to a rigid lattice [59]. The effect of thermal expansion is, however, not included. The configurational contribution to the free energy will therefore, in our modeling, effectively be volume independent, i.e., F conf . With a CE one can compute the energy of a configuration in a very efficient way, enabling the calculation of thermodynamic averages in MC simulations.
The grand canonical (GC) ensemble as implemented in the MCHAMMER module of ICET [49] is used to compute the configurational part of the free energy. In the GC ensemble the temperature T and carbon chemical potential μ C are fixed and the probability of being in a configuration σ is given by where E (σ) is the energy for configuration σ with N C carbon atoms. The standard Metropolis MC technique [60] is used, where carbon atoms are flipped to vacancies and vice versa, in order to generate configurations according to P(σ ). This allows averages such as E (σ) and N C , and hence also y = N C /N mix to be calculated. For a fixed T we can thus obtain μ C as a function of y, i.e., μ C = μ C (y, T ). The chemical potential μ C is given by the partial derivative of the free energy with respect to the number of carbon atoms, N C , i.e., μ C = ∂G/∂N C . By integrating we thus obtain In terms of the free energy per atoms this can be rewritten as and, hence, where we have used that, by definition, F conf (0, T ) ≡ 0.
FIG. 1. Free energy per atom of the stoichiometric (y = 1) γ -WC and δ-WC phases. The dotted line represents the free energy from electronic excitations. The dashed line represents the vibrational free energy (EHM + FEP) excluding thermal expansion. The solid line represents the total free energy G α (y = 1, T ), which includes electronic and vibrational contributions as well as thermal expansion. For δ-WC the total free energy is shifted by the volumerelaxed total energy of δ-WC. For the cubic phase, all free energies are shifted with the relaxed total-energy difference between the two WC phases, i.e., E 0 γ − E 0 δ .

Dilute limit
In a system of noninteracting vacancies, i.e., the dilute limit, the configurational entropy can be obtained from the ideal mixing entropy The configurational free energy per atom in the dilute limit can thus be written as where the factor 1/(1 + y) ensures that the free energy is per number of atoms.

A. δ-WC and γ-WC at y = 1
To construct the W-C phase diagram the free energy has to be derived for the δ-WC and γ -WC phases as functions of temperature and composition. We start with the case y = 1, i.e., N C = N W = N mix and no vacancies are present. In this case, the configurational contribution to the free energy vanishes and we only need to consider the electronic and vibrational contributions.
In Fig. 1 the free energy per atom of the γ -WC and δ-WC phases is shown as a function of temperature. At 0 K the energy difference between the two phases is E δ→γ = 0.454 eV/atom and the hexagonal δ-WC phase is the stable phase.
If we include the electronic contribution the free energy slightly decreases as a function of temperature, mainly for the γ -WC phase. The vibrational contribution is substantially larger and moreover the γ -WC phase shows a stronger temperature dependence, such that at high temperatures the γ -WC phase becomes the stable phase. If we only include the vibrational contribution without thermal expansion this transition occurs at T δ→γ = 3740 K. If we add the electronic contribution the transition temperature is shifted down to T δ→γ = 3340 K. While the electronic contribution is small, it thus still has a sizable influence on the transition temperature. For the total free energy, which also includes thermal expansion, the transition temperature is reduced even further to T δ→γ = 2780 K.

Electronic contributions
From Fig. 1, we see that the electronic contribution to the free energies is rather small. At 3000 K the electronic free energy of γ -WC and δ-WC is reduced by 57 and 18 meV/atom, respectively.
This can be understood from the EDOS which is shown in Fig. 2 at the respective 0-K volumes. The electronic free energy is approximately proportional to the EDOS at the Fermi energy n el (ε F ). It is evident from Fig. 2 that the γ -WC phase has a higher EDOS at the Fermi level than δ-WC, which exhibits a pronounced dip near the Fermi energy. This reflects the fact that γ -WC is more metallic, while δ-WC is more covalent. The EDOS per atom at the Fermi energy n el (ε F ) for γ -WC and δ-WC is 0.572 and 0.160 1/eV, respectively, which, according to Eq. (9), corresponds to the electronic freeenergy contributions −62.9 and −17.6 meV/atom at 3000 K. This is quite close to the more accurate values obtained from the full EDOS together with Eqs. (5) and (8).

Vibrational contribution
The vibrational contribution to the free energy is considerably larger compared to the electronic contribution. First, let us focus on the vibrational and thermal characteristics of the δ-WC phase. In Fig. 3(a) the phonon dispersion relation for 0 K (i.e., HA) and 3000 K (i.e., EHM) is shown, where the phonon dispersion at 3000 K includes thermal expansion. Two distinct bands are obtained arising from the large mass difference between C and W atoms. The low-frequency band corresponds essentially to W atom vibrations, while the high-frequency band primarily arises from C atom vibrations. There is a clear softening of the modes when the temperature is increased and this frequency renormalization is mainly due to thermal expansion.
Consider now the γ -WC phase. The phonon dispersion relation for γ -WC at 0, 1000, and 3000 K is shown in Fig. 3(c). As in the case of the δ-WC phase, we obtain two distinct bands arising from the large difference in mass between C and W atoms. Further, at 0 K, i.e., using HA, there is clearly an unstable phonon mode at the X point, a mode with imaginary frequency. However, at 1000 and 3000 K the imaginary frequencies are no longer present, meaning the γ -WC phase has been dynamically stabilized.
The potential energy surface (PES) along the unstable mode at the X point is explored within the HA at 0 K, EHMs at 1000 and 3000 K, and DFT ( Fig. 4) for the relaxed volume at 0 K. It reveals a double-well landscape which leads to an imaginary mode in the HA approximation, which is stabilized using EHMs. The PESs constructed with EHMs deviates significantly from the DFT-PES, reflecting the strong anharmonicity. Moreover, for this mode the EHMs for 1000 and 3000 K are quite different, which illustrates the need for temperature-dependent EHMs when determining the free energy over a large temperature range. At 3000 K the displacements sampled are larger than at 1000 K and a stiffer EHM is obtained. Consequently, the EHM for 3000 K has a higher frequency compared to the EHM for 1000 K. The increase in frequency at the X point with temperature is visible in the low-frequency (W-dominated) band in Fig. 3(c). While Fig. 3(c) includes the effect of thermal expansion, which generally leads to a lowering of frequencies with temperature, the effect from the instability at the X point is dominating in this case. According to our calculations, the lowest mode at the X point is stabilized somewhere between 500 and 1000 K. Temperatures lower than 750 K are, therefore, omitted from Fig. 1.
From the VDOS in Fig. 3(b), it is clear that, on average, the phonon frequencies are higher in δ-WC than in γ -WC. The free energy is given by the product of the VDOS and f vib (ω, T ) and the decrease in free energy is, therefore, larger for low frequencies than for high frequencies. Consequently, the vibrational free energy of γ -WC will decrease more strongly with temperature than the one of δ-WC, which is also seen in Fig. 1.
For more information regarding the fitting of the EHMs, see Supplemental Material [61].

Formation free energy
Based on the result for the free energies G α (T ) = G α (y = 1, T ) in Fig. 1, we can now determine the free energy of formation for the δ-WC and γ -WC phases according to where G bcc-W (T ) and G gr (T ) are the free energy per atom, or chemical potential, for body-centered cubic tungsten and graphite, respectively. Data for these pure phases are presented in Appendix A. In Fig. 5 the result as a function of temperature is shown and we compare with experimental data from Ref. [62].
The formation free energy for the δ-WC phase is rather independent of temperature, indicating that the temperature dependence of the vibrational free energy contribution in δ-WC is similar to body-centered cubic (bcc) W and graphite. While δ-WC is rather stiff with correspondingly high frequencies, γ -WC is considerably softer. This explains the much stronger temperature dependence of the formation free energy of the latter and the fact that it is decreasing as function of temperature.
The computed formation free energy for δ-WC agrees very well with the experimental data in Ref. [62]. DFT together with the PBE approximation for the exchange-correlation functional thus gives a good description of the δ-WC phase. In Table I, we compare our data for the formation enthalpy H form α and entropy S form with the experimental data at T = 1500 K from Ref. [62]. We also compare the computed transition temperature T δ→γ with the corresponding experimental result from Ref. [16]. In view of the sensitivity of T δ→γ to the slope of the free energies (see Fig. 1), we also view this agreement to be fully satisfactory.
In Table I, we also include results calculated with the same methodology but using the vdW-DF-cx exchange-correlation TABLE I. Formation enthalpies and entropies for stoichiometric γ -WC and δ-WC at 1500 K, given in eV/atom and k B /atom, respectively, as well as the transition temperature between the δ-WC and γ -WC phases in K. The experimental data are from Refs. [62] and [16]. functional [63]. In this case, the agreement for the formation enthalpy is somewhat less satisfying while the agreement for formation entropy and transition temperature is improved.

B. γ-WC at y 1
For y < 1, the configurational contribution to the free energy has to be considered, in addition to the electronic and vibrational contributions. The γ -WC phase can contain a large amount of carbon vacancies and to obtain the configurational part of the free energy an alloy CE was constructed [49]. To this end, about 520 structures with optimized shape, volume, and positions were selected for the cubic γ -WC phase and the clusters included in the final CE used cutoffs of 8, 8, 6.5, and 6 Å for pairs, triplets, quadruplets, and quintuplets, respectively, amounting to a total of 83 orbits (symmetryinequivalent clusters). For more information regarding the training structures for the CE, see Supplemental Material [61]. The final CE was trained using the automatic relevance detection regression (ARDR) algorithm resulting in 43 nonzero ECIs and a cross-validation (CV) score of 13 meV/N mix . From the ECIs (Fig. 6), we note that all the pair interactions are positive meaning that carbon-carbon and vacancy-vacancy bonds are energetically unfavorable compared to carbon-vacancy bonds. The strongest ECI is the one corresponding to the second-nearest-neighbor pair, as in the case of cubic TiC [64].

Low-temperature structures
The CE was used to determine the configurational contribution to the free energy at high temperatures. We also investigated possible structures predicted by the CE at lower temperatures. This was done for a range of carbon concentrations using simulated annealing in the canonical ensemble, meaning the MC simulations were evolving with decreasing temperature at fixed concentration. For more information regarding the MC simulations, see Supplemental Material [61]. These simulations reveal that at a certain temperature the system transitions from a disordered state to an ordered state. In order to illustrate these order-disorder transitions, the configurational contribution to the heat capacity was computed as a function of both temperature and concentration (Fig. 7). Peaks in the heat capacity reveal disorder-order transitions, which are particularly noticeable for y = 7 8 , 3 4 , and 1 2 . For W 8 C 7 and W 2 C the most pronounced transition occurs around 1800 and 1200 K, respectively. These types of ordered structures are common in nonstoichiometric group IV (Ti, Zr, Hf) [3,64,65] and group V (V, Nb, Ta) [3,65,66] cubic carbides. In WC a lower tungsten carbide phase W 2 C with hexagonal structure is instead formed [16]. However, ordered cubic structures in WC can become important in thin-film structures, which are present at considerably lower temperatures [31][32][33][34], and in ultrafine WC nanoparticles [35].

Free energies
Next, we introduce the mixing free energy per atom according to where data for G bcc-W (T ) and G gr (T ) can be found in Appendix A. For y = 1 the mixing free energy is equal to the formation free energy, introduced in Eq. (25). In Fig. 8(a), we present the result for the mixing free energy of the γ -WC phase, G mix γ (y, T ), as a function of composition and at different temperatures, with and without thermal expansion. A workflow diagram for calculating the free energy of γ -WC can be found in the Supplemental Material [61]. The mixing energy, corresponding to the mixing free energy at 0 K shown in Fig. 8(a), exhibits a pronounced dependence on the vacancy concentration y, decreasing by 0.25 eV/atom as the concentration changes from y = 1 to 0.5, i.e., as the relative amount of W increases. We conclude that at T = 0 K 033804-7 a more energetically stable configuration is formed at the reduced carbon concentration of y = 0.5 compared with y = 1. This is consistent with the repulsive carbon-carbon interaction that is evident from the ECIs of the CE. At higher temperatures the mixing free energy G mix γ (y, T ) is reduced, shows a weaker y dependence, and its minimum is shifted slightly to higher carbon concentrations, as seen in Fig. 8(a).
Let us now focus on the different temperature-dependent contributions to the free energy, which are shown in Fig. 8(b) at 2000 K, excluding the effect from thermal expansion. The electronic contribution is calculated, as described in section II C, for various concentrations y and the result is shown as a blue line in Fig. 8(b). We find that the electronic excitations give the smallest contribution to the mixing free energy. There is a small increase of the electronic contribution at high carbon concentrations, indicating a drop in the EDOS near the Fermi level as the carbon concentration increases to y = 1.
The configurational contribution is obtained using the CE (see section II E). As seen in Fig. 8(b), this contribution exhibits more variation as a function of composition and is generally larger than the electronic contribution, but still quite small. The difference between our model and the ideal mixing model according to Eq. (18) is primarily due to the energetic contributions from the CE model. The entropic contribution at this high temperature is similar in the two different models but the CE model exhibits more structure with variations in the composition. We also notice a small increase at the concentrations y = 1 2 , 3 4 , and 7 8 , which coincide with the ordering tendencies shown in Fig. 7.
The largest temperature-dependent contribution to the free energy comes from the vibrations. The green circles in Fig. 8(b) show the result for nine selected high-symmetry ground states at 2000 K, obtained using EHMs in combination with FEP. At y = 0.5, where an energetically stable configuration is formed, the frequencies are high and comparable to the frequencies for the pure elements. This implies a rather small contribution to the mixing free energy from the vibrations at y = 0.5. However, for y = 1, where a less energetically stable configuration is formed, the frequencies become softer as seen in Fig. 9(a), where a clear shift to lower frequencies is obtained both for the lower W-dominated and higher C-dominated groups of bands. Lower frequencies imply a larger vibrational free-energy contribution to the mixing free energy and hence the difference between y = 1 and 0.5 is large (0.18 eV/atom).
The above results from the vibrational DOFs [green circles in Fig. 8(b)] were obtained using high-symmetry groundstate structures. However, the vibrational properties for a ground-state configuration can be quite different compared to a disordered configuration [67]. To investigate this effect, calculations in the HA were carried out for the ground state and a few disordered (high-temperature) structures at y = 0.5 and 0.75. The disordered structures were chosen from MC 033804-8 simulations at 3000 K. We find that for y = 0.5 the HA free energy differs by roughly 35 meV/atom between the ground state [green circles in Fig. 8(b)] and disordered structures [red circles in Fig. 8(b)] with a spread of about ±15 meV/atom at 2000 K. At y = 0.75, the difference is smaller and for y = 1 there is, by definition, no difference.
To understand the differences between ground-state and disordered structures, we show in Fig. 9(b) the VDOS for the y = 0.5 ground state and one representative disordered structure. The disordered structure has an increased VDOS for low frequencies in the W-dominated (low-frequency) bands compared to the ground state. This explains the lower vibrational free energy for the disordered configurations compared to the ground states.
In order to include the above-described effect over the entire concentration range, we fitted a second-order polynomial to the vibrational free-energy data at y = 0.5, 0.75, and 1, as shown by the red line in Fig. 8(b). We can now add the electronic (blue line), configurational (orange line), and vibrational (red line) contributions to obtain the total temperature-dependent contribution to the mixing free energy. This is added to the 0-K contribution and we obtain the total contribution at different temperatures, shown by dashed lines in Fig. 8(a).
The final step is to minimize the total free energy with respect to volume to obtain the free energy at zero pressure G γ (y, T ) (see section II B), which is done in order to take into account thermal expansion. This is done by first deriving the shift of the free energy for three different compositions y = 0.5, 0.75, and 1, at different temperatures. This shift is added and interpolated across the entire concentration range. The final result is presented in Fig. 8(a) at two different temperatures. We find that the mixing free energy at T = 2000 K has a minimum at y = 0.52 and the mixing free energy is negative over the entire concentration range. At T = 2800 K, the minimum is shifted to y = 0.58, at which point the mixing free energy equals −0.262 eV/atom. These final results for the mixing free energy of the γ -WC phase enable the construction of the binary W-C phase diagram in Fig. 13.

Stability of the γ-WC phase
In section III A, we found that the stoichiometric (y = 1) γ -WC phase is dynamically unstable at low temperatures while it is dynamically stabilized above T ≈ 750 K. The phase can also be stabilized by introducing vacancies, i.e., for y < 1. To demonstrate this effect, we determined the VDOS in the HA for the nine high-symmetry ground states selected above. The results for the structures at y = 0.5, 0.75, and 1 are shown in Fig. 10, which shows that the γ -WC phase is dynamically stable for y 0.75. Hence, the γ -WC phase can be dynamically stabilized both by increasing the temperature and by introducing carbon vacancies.

C. δ-WC at y 1
In contrast to the γ -WC phase the C vacancy concentration in the δ-WC phase is very small. The C vacancies in the δ-WC phase can, therefore, be treated as noninteracting, i.e., in the dilute limit.
To this end, we introduce a term for the free energy without the configurational contribution, but including electronic and The effect of thermal expansion will be disregarded and hence the dependence on the volume v is left out in the above definition. The vacancy formation energy can then be obtained by computingF δ (y, T ) for two sufficiently large systems with and without a vacancy, respectively. The free energy of vacancy formation can then be written as where we have assumed equilibrium with graphite and y ref = N C /N mix . Assuming ideal mixing and using the fact that the vacancy concentration is small we obtain the standard expression for the vacancy concentration y(T ). We calculatedF δ (y, T ) for a system with N mix = 216 lattice sites, both for a system without vacancies (y = 1) and for a system with one vacancy (y ref = 215/216), treating the vibrations in the HA. Further, the difference in electronic free energy between these two systems was below the accuracy of our DFT calculations and is, therefore, regarded as zero. A workflow diagram for calculating the free energy per atom of δ-WC can be found in the Supplemental Material [61].
The C vacancy formation energy G vac δ (T ) is shown together with the corresponding vacancy concentration in Fig. 11. The vacancy formation energy decreases with increasing temperature indicating a softening of the phonon modes due to the vacancy. Furthermore, the vacancy concentration reaches a value of around 0.5% at 3000 K.

D. Phase diagram
In the previous sections, we presented how the free energy for each phase α as a function of composition y and temperature T can be calculated. Now, we need to find the optimal combination of phases as a function of overall concentration c = N C /(N C + N W ), and temperature T , as described in section II A. To this end, we construct the convex hull for different temperatures, i.e., the optimal phase combination for each c at a fixed T . The total free energy is then minimized for fixed T according to Eq. (1). Figure 12 shows the convex hull as a function of C concentration at T = 2500 K plus the mixing free energy of the δ-WC and γ -WC phases. For low C concentrations, up to about c = 32%, one obtains a two-phase region between γ -WC and bcc-W. From around c = 32% to 37% only γ -WC is stable. Further, between c = 37% and 50%, we obtain another two-phase region between δ-WC and γ -WC. For a very limited range below c = 50% only δ-WC is stable. Finally, above c = 50% the equilibrium is between δ-WC and graphite.
By constructing the convex hull for many temperatures and marking the transition concentrations, the phase diagram may be constructed by interpolating along these transition lines (Fig. 13). In the phase diagram the very narrow δ-WC homogeneity region is included, but cannot be visually re- solved at this scale. The convex hull construction is quite sensitive to small changes in the mixing free energy in regions where the mixing free energy lies very close to the convex hull. This is most noteworthy for the γ -WC + W to γ -WC and the γ -WC to γ -WC + C transitions. These transitions are therefore drawn as dashed lines in order to indicate their higher uncertainty.

IV. DISCUSSION
In this paper we have constructed a part of the binary phase diagram for the W-C system from first principles (see Fig. 13). We will now comment on the limitations in our computational modeling approach and make a comparison with the experimental phase diagram [16].

A. Exchange-correlation approximation
We have used the semilocal PBE [43] approximation to describe the exchange and correlation effects within the density functional theory (DFT). This functional has become the standard approximation used in materials science applications. To test the sensitivity of the choice of exchange-correlation functional on the results for free energy we have done the complete analysis for the stoichiometric γ -WC and δ-WC phases using the van der Waals exchange-correlation functional vdW-DFcx [63]. (For details, see Supplemental Material [61].) This functional captures nonlocal correlations [68,69] in combination with a consistent description of exchange [63]. It has been established that it provides a description of the thermophysical properties of nonmagnetic transition metals [70] that is as least on par with but usually exceeds other constraint-based functionals as PBE [43] and PBEsol [71] and, more recently, the phase behavior of a typical W containing metal-alloy system has been derived successfully using the vdW-DF-cx functional [72].

033804-10
The formation enthalpy for the hexagonal δ-WC phase is found to be somewhat better within PBE compared to vdW-DF-cx, while the opposite is the case for the formation entropy (see Table I). Using the vdW-DF-cx functional the δ-WC to γ -WC transition temperature T δ→γ for the stoichiometric phases changes from 2780 to 2890 K, where the latter is closer to the experimental value of 2993 K. The reason is that the vdW-DFcx functional improves the description of thermal expansion at high temperatures (see Supplemental Material [61]). For temperatures near these transition temperatures there is an offset in the energy difference between the γ -WC and δ-WC phases of about 20 meV/atom when comparing vdW-DF-cx to PBE. Hence, at the transition temperature for PBE, i.e., 2780 K, the free energy for γ -WC is about 20 meV/atom higher compared with δ-WC when using vdW-DF-cx. This comparison indicates the accuracy that can be obtained with DFT calculations.

B. Modeling limitations
There are also some limitations in our modeling approach. We essentially separate the electronic, vibrational, and configurational DOFs although they are coupled in reality. The relevant temperatures in this work are close to the melting point of WC and the decoupling of vibrational and configurational DOFs may become somewhat questionable.
We base the vibrational contributions in the cubic γ -WC phase on using some structures, ground-state structures as well as representative high-temperature structures. In this way we try to combine the sampling of vibrations and disordering (configurational DOFs) in a realistic and efficient way. From the discussion in section III B 2 we conclude that a reasonable accuracy that can be obtained for the vibrational contribution is ±15 meV/atom. Moreover, we are likely not able to capture the full anharmonicity near the melting temperature. It is likely that the error near the melting point is on the order of ±10 meV/atom [57]. Furthermore, the configurational free energy in γ -WC is calculated with an alloy cluster expansion which has a typical error of ±10 meV/atom (see section III B).
Our treatment of the vacation formation energy in δ-WC could be improved. It is treated in the harmonic approximation although anharmonicity likely is an important effect at high temperatures [73,74]. However, the level of accuracy needed to include anharmonicity is difficult to achieve with EHMs from AIMD simulations. Likely, anharmonicity will lower the vacancy formation energy slightly, thus increasing the vacancy concentration. However, there will only be minor changes to the phase diagram, mainly a slightly larger δ-WC homogeneity region.

C. Comparison with the experimental phase diagram
We will now compare our data with the most up to date experimental phase diagram for the W-C system [16] for compositions and temperatures where the cubic γ -WC phase is stable. This is also the part of the experimental phase diagram that is most controversial [16].
For the stoichiometric WC composition, i.e., y = N C /N W = 1, we predict that γ -WC + graphite is stable above 2700 K (see Fig. 13). Here, the γ -WC phase contains roughly 29% C vacancies as can be seen from the vertical dashed line which indicates the transition between the homogeneous γ -WC region and the γ -WC + graphite region. We may, therefore, refer to it as γ -WC 0.71 . Moreover, below 2700 K we predict that δ-WC with roughly 0.1% carbon vacancies and graphite are the stable phases.
In the experimental phase diagram [16] the corresponding phase transition happens at 2993 K, instead. Further, in this case, the γ -WC phase is stoichiometric, i.e., no graphite forms, and the δ-WC phase has roughly 1% carbon vacancies. Hence, we underestimate the δ-WC + graphite to γ -WC phase transition with about 290 K, which corresponds to about 10% of the experimental temperature. Moreover, at 2993 K our γ -WC + graphite convex hull lies only 37 meV/atom below the δ-WC + graphite convex hull. This is a small energy and could very well lie within the possible error of our method. Furthermore, the free-energy difference between the convex hull for composition y = 1 at 2780 K and stoichiometric γ -WC is about 15 meV/atom. It is, consequently, almost equally likely to have stoichiometric γ -WC as γ -WC 0.71 + graphite for composition y = 1 above 2780 K. In the former case this would move the horizontal dashed line at y = 0.71 to y = 1 which would be in agreement with the experimental phase diagram.
The stability region for γ -WC + δ-WC is predicted by the experimental phase diagram to be narrow near the δ-WC + graphite to γ -WC transition temperature (2993 K) and widen down to 2789 K where lower tungsten (semi)carbides are formed. At this temperature γ -WC + δ-WC is stable in the 0.65 < y < 1 composition range (close to y = 1 δ-WC has a tiny homogeneity region).
We predict that the γ -WC + δ-WC region has a width of 0.71 < y < 1 at the δ-WC + graphite to γ -WC + graphite transition temperature and that roughly 200 • below, it has widened to 0.6 < y < 1. We, therefore, get a rather good agreement on the width of the γ -WC + δ-WC region a bit below the δ-WC + graphite to γ -WC + graphite transition temperature.
Furthermore, there are small differences in the free energies around the triple point joining γ -WC + δ-WC, homogeneous γ -WC and γ -WC 0.71 + graphite at composition y = 0.71 and temperature 2700 K. This means that the position of the triple point and, consequently, also the width of the γ -WC + δ-WC region has a high uncertainty.
Since we do not include any tungsten semicarbides in our phase diagram our γ -WC + δ-WC region continues down to roughly 1800 K where δ-WC + bcc-W are the stable phases.
We conclude that to match the experimental phase diagram our computed data have to be shifted with up to about 40 meV/atom. We have estimated that the accuracy of the exchange-correlation approximation in the DFT calculations could be about ±20 meV/atom. The limitations in the modeling approach leads to uncertainties of the order ±15 meV/atom, ±10 meV/atom, and ±10 meV/atom, for the two vibrational contributions and the cluster expansion, respectively, which give rise to an uncertainty of the order ±20 meV/atom. This implies that the inaccuracies of the computed phase diagram certainly lie within the uncertainties related to our computational and modeling approach.

V. CONCLUSION
The hexagonal and cubic WC phases are modeled using DFT and a part of the binary W-C phase diagram is constructed. We include electronic, vibrational, configurational, and thermal expansion effects and study how they vary with the carbon content in order to qualitatively and quantitatively understand the thermodynamic stabilization of the cubic phase at high temperatures.
At 0 K, the energy of the stoichiometric cubic γ -WC phase is found to be considerably higher, about 0.5 eV/atom, than the corresponding hexagonal δ-WC phase, and moreover the γ -WC phase is dynamically unstable. Stabilization occurs above ∼800 K. Due to the considerably softer vibrational spectrum for the γ -WC phase compared to the δ-WC phase, the γ -WC phase becomes thermodynamically stable around 3000 K. The transition temperature T δ→γ between the stoichiometric δ-WC phase and the stoichiometric γ -WC phase is predicted to be 7.1% and 3.4% lower than the experimental value T = 2993 K, using the PBE and vdW-DF-cx functionals, respectively. Further, about 200 K below T δ→γ the predicted width of the γ -WC + δ-WC region is 0.6 < y < 1 in our theoretical phase diagram, similar to the corresponding experimental result 0.65 < y < 1.
Overall, we view the accuracy of the presented methodology as sufficient for studying also more complex interfacial structures in doped cemented carbides, so-called complexions [39,40], which is motivated by their technological relevance and lack of experimental information.
In conclusion, this study shows that it is feasible to theoretically derive a part of the phase diagram from first principles for tungsten carbide, including the cubic γ -WC phase, which is dynamically unstable at low temperatures.
given in the provided PAW functional (POTCAR) files, was used.
The reference energies for the alloy CE, the electronic freeenergy calculations, and forces used in HA were calculated with projection operators evaluated in the reciprocal space, i.e., LREAL=.FALSE.. Further, the AIMD simulations were done with projection operators evaluated in real space, i.e., LREAL=Auto. Moreover, after a careful analysis of the free energies from EHMs generated from AIMD snapshots, it was concluded that LREAL=Auto was sufficient, likely because of the large displacements (forces) in the snapshots. More details regarding the AIMD simulations can be found in the Supplemental Material [61].
When calculating the reference energy used for training the cluster expansions and the electronic free-energy contribution the ionic positions were optimized until all forces were below 0.01 eV/Å. Further, when calculating the vibrational free energy the ionic positions were optimized until all forces were below 0.005 eV/Å. Furthermore, all calculations were done without spin polarization since none of the systems in the study showed any magnetization when including spin polarization.