Thermodynamics of Plutonium Monocarbide from Anharmonic and Relativistic Theory

: Thermodynamics of plutonium monocarbide is studied from ﬁrst-principles theory that includes relativistic electronic structure and anharmonic lattice vibrations. Density-functional theory (DFT) is expanded to include orbital-orbital coupling in addition to the relativistic spin-orbit interaction for the electronic structure and it is combined with anharmonic, temperature dependent, lattice dynamics derived from the self-consistent ab initio lattice dynamics (SCAILD) method. The obtained thermodynamics are compared to results from simpler quasi-harmonic theory and experimental data. Formation enthalpy, speciﬁc heat, and Gibbs energy calculated from the anharmonic model are validated by direct comparison with a calculation of phase diagram (CALPHAD) assessment of PuC and sub-stochiometric PuC 0.896 . Overall, the theory reproduces CALPHAD results and measured data for PuC rather well, but the comparison is hampered by the sub-stoichiometric nature of plutonium monocarbide. It was found that a bare theoretical approach that ignores spin-orbit and orbital-orbital coupling (orbital polarization) of the plutonium 5f electrons promotes too soft phonons and Gibbs energies that are incompatible with that of the CALPHAD assessment of the experimental data. The investigation of PuC suggests that the electronic structure is accurately described by plutonium 5f electrons as “band like” and delocalized, but correlate through spin polarization, orbital polarization, and spin-orbit coupling, in analogy to previous ﬁndings for plutonium metal. thermodynamic properties of ideal-stoichiometry PuC as well as PuC 0.869 (speciﬁc heat). In the case of the Gibbs energy, we only consider the ideal-stoichiometry PuC for an appropriate comparison with the ﬁrst-principles results.


Introduction
The physics of the actinides and their compounds are fascinating but also somewhat controversial. The controversy is primarily focused on the nature of the actinide 5f electrons and their degree of correlation. Strong electron correlation manifests itself as localization of the 5f electrons on the actinide atom and for the actinide-oxide compounds this localization leads to band gaps in the electronic structure [1,2]. On the other hand, one finds weaker electron correlation for the early elemental actinide metals, thorium through plutonium. Particularly for the first four, there is now consensus that they are appropriately described by band-like 5f electrons [3,4]. In terms of theoretical approaches, weak or intermediate electron correlation implies that density-functional theory is a viable starting point and methods aimed for stronger electron correlations, assuming explicitly an intra-atomic Coulomb interaction with a Hubbard U parameter (density-functional theory (DFT) + U), is not necessary. When it comes to plutonium metal the electronic structure is still debated in the literature. There are views that the 5f electrons are not strongly correlated and that no Hubbard U (U = 0) [5,6], or merely a small value (U~1 eV) [7], is appropriate. However, there are also those who contrarily believe that the 5f electrons are strongly correlated, nearly localized, and suitably characterized by a very large U parameter (U~4-4.5 eV) [8,9]. Some of these models are convoluted and require assumptions of parameters such as the debated U and others. One cannot ignore, however, that a more transparent and less complicated model that relies on no uncertain parameters describes plutonium very well [5], including its magnetic profile [10]. Of course, such a straightforward theory is preferred if it sufficiently explains the experimental observations.
Regarding the actinide mononitrides and monocarbides, they appear to lie between the oxides and the elemental metals in terms of the 5f-electron correlation strength [11]. In the nitrides and carbides, one does not encounter band gaps as in the oxides and they are also not insulators but metallic. Analogous with the actinide oxides, they are magnetic, but the detailed magnetic structures are not known for all of them.
The actinide mononitrides and monocarbides form in the rock-salt (B1 or NaCl) structure, i.e., a cubic phase with the two atoms in the (000) and ( 1 2 1 2 1 2 ) positions, respectively. Their simple crystal structure and the fact that the 5f electrons are less correlated than in their oxide counterparts make these compounds well suited for DFT studies. One complication for modeling, however, is that vacancies (of nitrogen or carbon) tend to form in the materials, which, thus, become nitride or carbide deficient. Consequently, they form as AnN 1−x or AnC 1−x (An is actinide) compounds, where x is small (~0.1) but significant. In recent literature we find that the electronic structures for several AnC compounds have been studied with quantum-chemical methods [12].
In addition to the fundamental-science interest in the actinide mononitrides and monocarbides, they are compelling from a practical and applications perspective. The technological interest arises from their potential use as advanced nuclear fuels for fast-breeder reactors. These materials have good mechanical characteristics, but also possess superior thermophysical properties [13] such as high melting temperatures, high density of heavy atoms, and high thermal conductivity. In spite of the interest, there have, unfortunately, not been as many studies on them as there have been on the actinide oxides. Robust modeling, particularly at high temperature, is therefore much needed. In the present report we focus on ground-state and high-temperature thermodynamical properties of plutonium monocarbide, PuC, from first-principles theory.
Recently, we studied on uranium mononitride [14] where the thermodynamic modeling was based on DFT electronic structure combined with lattice dynamics that accounts for anharmonic lattice vibrations. For UN, we compared and validated our first-principles model directly with experiments in addition to results from a calculation of phase diagram (CALPHAD) assessment of the available measured data. In the present work for PuC, we adopt the same modeling approach but recognize that the plutonium 5f electrons pose a greater challenge for the theory. Therefore, we go beyond our previous treatment of the electronic structure in UN and now include relativistic effects and an extension of DFT that addresses orbital polarization, a mechanism that is known to be important for plutonium [5,15]. It turns out that these additional electron correlations are necessary for realistic PuC Gibbs energies. For comparison, and to confirm our first-principles model for PuC, we carry out CALPHAD calculations of the Gibbs energies, heat capacities, and formation enthalpy, utilizing a published thermodynamic database [16] and the Thermo-Calc software [17,18].
In the following Sections 2-4 we briefly detail the DFT implementations and CALPHAD method and continue by showing our results and provide context in a summary and discussion section.

Electronic-Structure Methods
We apply three methods, each with their own strengths, for calculating the electronic structure of PuC. Two of them are "all-electron" approaches (i.e., no core-electron approximation), that in one case is implemented with a "full potential" (i.e., no geometrical or structural approximations). The other all-electron method is constructed with a Green's function technique that allows for realistic alloy and disorder treatment. The third method is not all-electron but a plane-wave pseudopotential method that is efficient and fast for calculating forces on atoms in large cells. All three approaches rely on DFT and the generalized gradient approximation (GGA) [19] for the electron exchange and correlation energy functional.
First, we spotlight the all-electron full-potential linear muffin-tin orbital (FPLMTO) method for which implementation-details can be found the literature [20]. It adopts no approximations for the core states that exist at deeper energy levels than the valence states. This is a more accurate treatment than that of the plane-wave methods where the core electrons are replaced by a pseudopotential.
Some quantities are expanded in the series (basis functions, electron densities, and potentials) of spherical harmonics inside non-overlapping spheres centered at each atomic position. The radial part of the basis functions inside these spheres are calculated from a wave equation that addresses all relativistic corrections including spin-orbit coupling for d and f states but not for the p states, as is appropriate for actinides [6]. The orbital-polarization mechanism, considered here for the d and f states [6,21], is not included in conventional density-functional theory but is known to be important for some f-electron systems, particularly plutonium [5,15]. Generally, the set-up parameters of the present calculations for PuC are close to those applied for plutonium metal [5]. The FPLMTO method is employed to calculate the PuC formation enthalpy, the elastic constants, and the so-called "cold curve", i.e., the total-energy variation with atomic volume. The latter provides the fundamental input to the Debye-Grüneisen quasi-harmonic simulations (see below). The elastic constants are calculated, adopting conventional strains for cubic crystals while the shear modulus is the Voigt-Hill-Reuss average [22] of the single-crystal moduli. Furthermore, the presented electronic density of states and the related simulated photoelectron spectra are obtained from FPLMTO. Lastly, the Racah 5f parameters, that enter in the orbital-polarization formalism [21], are self-consistently calculated with FPLMTO for use as fixed parameters (~0.05 eV) in our plane-wave calculations [6].
For the computationally more demanding supercell calculations that couple to self-consistent ab initio lattice dynamics (SCAILD) [23], we utilize an efficient electronic-structure approach. Namely, the pseudopotential plane-wave Vienna ab initio simulation package (VASP) with the projector-augmented-wave method with the plane-wave basis set, as implemented in VASP [24][25][26]. The computational set-up is defined by an energy cut-off of 400 eV and an energy convergence of 100 eV. The VASP approach furthermore includes non-collinear magnetism with spin-orbit coupling and orbital polarization, as implemented recently [6].
PuC is sub-stoichiometric, as mentioned, and the sensitivity to the deviation of the (number of atoms) Pu/C ratio from unity is explored by VASP as well as utilizing a technique that incorporates efficient alloy theory within the coherent-potential approximation (CPA) [27]. The exact muffin-tin orbital (EMTO) method relies on Green's function formalism where the one-electron potential is represented by optimized overlapping muffin-tin potential spheres [28]. The EMTO-CPA [29] is well suited to explore compositional disorder and sub-stoichiometric conditions. It can also be used to study the influence of randomly distributed vacancies on a sub-lattice (see the Summary and Discussion section below). In this report, we restrict ourselves to EMTO calculations without spin-orbit coupling.

Lattice-Dynamics Methods
Our main tool for lattice dynamics is the SCAILD methodology [23] that we often refer to as the self-consistent phonon method. The idea behind the scheme is to employ the small-displacement method to calculate the phonons in a first step and then apply thermally induced (fixed finite temperature) "frozen" phonons on the atoms and calculate corresponding DFT atomic forces. The next step is to recalculate new phonons informed by these DFT forces and repeat until convergence [23]. This self-consistent phonon approach is appealing because it correlates displacements of all atoms with each other (unlike the frozen-phonon method) and, therefore, it can account for strong phonon-phonon coupling and anharmonic behavior. We have used SCAILD for uranium mononitride and the cubic γ phases of uranium and plutonium in the past [5,14,30,31] with satisfying results. For SCAILD, one needs to define a supercell of the crystal structure and for PuC we specifiy a 3 by 3 by 3 supercell for a total of 54 atoms. The DFT forces are gathered from VASP calculations and the self-consistent scheme requires about 200 iterations to converge the lattice-vibration contribution to the free energy to the meV level. We perform several simulations: For each temperature (T = 600, 800, 1000, 1200, 1500, and 2000 K), we choose three to five (depending on temperature) lattice constants in order to determine the equilibrium volume on the basis of the Gibbs-energy minimum. To assess the importance of 5f-electron correlations beyond conventional DFT-GGA, we employ two approaches; one full treatment that includes spin-orbit coupling and orbital polarization "DFT+SO+OP" and one bare treatment that excludes these interactions "DFT".
In Figure 1 we show results from SCAILD for PuC at a temperature of 800 K, for the two levels of electronic-structure theory. The lattice constants are held fixed, corresponding to the zero-temperature equilibrium volumes (16.2 and 15.6 Å 3 , see Table 1 in Section 3.2).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 4 of 16 contribution to the free energy to the meV level. We perform several simulations: For each temperature (T = 600, 800, 1000, 1200, 1500, and 2000 K), we choose three to five (depending on temperature) lattice constants in order to determine the equilibrium volume on the basis of the Gibbsenergy minimum. To assess the importance of 5f-electron correlations beyond conventional DFT-GGA, we employ two approaches; one full treatment that includes spin-orbit coupling and orbital polarization "DFT+SO+OP" and one bare treatment that excludes these interactions "DFT".
In Figure 1 we show results from SCAILD for PuC at a temperature of 800 K, for the two levels of electronic-structure theory. The lattice constants are held fixed, corresponding to the zerotemperature equilibrium volumes (16.2 and 15.6 Å 3 , see Table 1 below).   We notice in Figure 1 that the SCAILD vibrational free energy is well converged and that the energy from the lesser theory (DFT) lies below that of the full treatment (DFT+SO+OP). This result is related to softness of some phonon branches in the DFT approach that we will discuss further in the Results section.
One can judge the importance of anharmonic phonons by comparing calculated thermal properties from anharmonic and lower-level quasi-harmonic theory. For this reason, we conduct a limited set of calculations employing two implementations of the Debye-Grüneisen quasi-harmonic approach [32]. The Debye-Grüneisen scheme is very efficient and only requires the DFT cold curve as input in addition to the atomic mass of the atoms in the material. There are, however, some We notice in Figure 1 that the SCAILD vibrational free energy is well converged and that the energy from the lesser theory (DFT) lies below that of the full treatment (DFT+SO+OP). This result is related to softness of some phonon branches in the DFT approach that we will discuss further in the Results section.
One can judge the importance of anharmonic phonons by comparing calculated thermal properties from anharmonic and lower-level quasi-harmonic theory. For this reason, we conduct a limited set of calculations employing two implementations of the Debye-Grüneisen quasi-harmonic approach [32]. The Debye-Grüneisen scheme is very efficient and only requires the DFT cold curve as input in addition to the atomic mass of the atoms in the material. There are, however, some assumptions within this model that must be made. First, Moruzzi et al. [32] suggested that the Debye temperature can be derived from this simple relationship: where r is the atomic radius and M the atomic mass in atomic units, B is the bulk modulus in kbar, and Const a scaling factor. The Const coefficient can be chosen to optimize agreement with experimental data or calculated assuming an isotropic material. Second, there are two philosophies regarding the formulation of the Grüneisen parameter. One that is supposedly [32] better at higher temperatures, γ HT , and one that is more suitable when comparing to experimental specific-heat data at lower temperatures, γ LT . They are referred to as Slater (γ HT ) and Dugdale-MacDonald (γ LT ) [33,34], respectively, and are defined as: here, P is the pressure and V the atomic volume. Third, the necessary analytical representation of the DFT-calculated cold curve is an approximation in itself. Moruzzi et al. [32] applied a Morse function [35] for this purpose but other forms can be considered, and the actual choice will influence the results of the Debye-Grüneisen model to an extent. In our quasi-harmonic calculations of the specific heat we utilize the Gibbs2 package [36] and our own implementation [37]. For Gibbs2, referenced below as "quasi-harmonic (a)", we use Const = 58.03, the low-temperature γ LT , and a third-order Birch-Murnaghan analytical form [38]. For isotropic materials one can calculate Const directly assuming that the Poisson's ratio is known [37] and for Gibbs2, ν = 0.25 is assumed. In the alternative treatment [37], referred to as "quasi-harmonic (b)", Const = 57.23, as dictated by our theoretical Poisson's ratio (ν = 0.256, see Table 1 in Section 3.2). Our implementation [37], furthermore, assumes the high-temperature expression for the Grüneisen parameter (γ HT ) and a Morse function for representing the DFT total energies.

CALPHAD Method
We apply the thermodynamic CALPHAD technique to compute the heat of formation, heat capacity, and Gibbs energy of PuC as functions of temperature. Our first-principles thermodynamic data are compared with the CALPHAD results and, in the case of the specific heat and formation enthalpy, with the experimental data too. Importantly, CALPHAD is able to, also, interrogate the significance of the sub-stoichiometric formation of PuC 1−x . Generally, the fundamental objective of the CALPHAD scheme is to model the Gibbs energy of individual phases pertaining to binary and ternary systems to optimally reproduce carefully reviewed phase diagrams and thermodynamic properties. From the Gibbs energy one computes phase stability and thermodynamic properties of multi-component systems [39][40][41].
Here, the CALPHAD data is compared to the first-principles modeling, but, additionally, we understand that the interaction between CALPHAD and theory improves the thermodynamic-modeling capability particularly for materials with many unknown variables. For example, ab initio results such as heats of formation can directly provide important constraints to the CALPHAD modeling framework in the absence of experimental data. Furthermore, we know that optimization of parameters and minimizing errors within in the CALPHAD technique is an inverse problem with infinite degrees of freedom [42]. As a consequence, many combinations of parameters chosen by the user can produce coinciding phase diagrams. The use of DFT-predicted properties related to the CALPHAD assessment constrains the optimization and certifies the resulting thermodynamical database, both in terms of phase stability and energetics. CALPHAD has been successful in predicting complex phase transformation in plutonium alloys [43,44] and first-principles-informed CALPHAD assessments for actinide systems have become customary in recent years [45][46][47][48].
Specifically, for the plutonium-carbide system, we utilize the CALPHAD assessment performed by Guéneau et al. [16] that integrates a critical review of available experimental data and previous CALPHAD assessment and apply it to calculate thermodynamic properties of ideal-stoichiometry PuC as well as PuC 0.869 (specific heat). In the case of the Gibbs energy, we only consider the ideal-stoichiometry PuC for an appropriate comparison with the first-principles results.

Electronic Structure
As mentioned in the introduction, one of the central questions regarding the electronic structure of any actinide compound is the nature of the actinide 5f electrons. For the elemental metals, up to americium, the 5f electrons can be regarded as bonding and itinerant and not localized as they are in americium and the following, heavier, actinides. For plutonium there is still a debate on this, but our view is that the 5f electrons are more delocalized than not and that opinion is based on a wealth of evidence [5]. Specifically, for plutonium monocarbide, it has been argued that the 5f electrons are less correlated than they are in both plutonium oxides and nitrides [11]. No theoretical method is currently able to predict with certainty the nature of the 5f electrons in these systems but model comparisons with experimental data provide clues.
In Figure 2 we show the (DFT+SO+OP) total electronic density of states (e-DOS) for PuC.
Specifically, for the plutonium-carbide system, we utilize the CALPHAD assessment performed by Guéneau et al. [16] that integrates a critical review of available experimental data and previous CALPHAD assessment and apply it to calculate thermodynamic properties of ideal-stoichiometry PuC as well as PuC0.869 (specific heat). In the case of the Gibbs energy, we only consider the idealstoichiometry PuC for an appropriate comparison with the first-principles results.

Electronic Structure
As mentioned in the introduction, one of the central questions regarding the electronic structure of any actinide compound is the nature of the actinide 5f electrons. For the elemental metals, up to americium, the 5f electrons can be regarded as bonding and itinerant and not localized as they are in americium and the following, heavier, actinides. For plutonium there is still a debate on this, but our view is that the 5f electrons are more delocalized than not and that opinion is based on a wealth of evidence [5]. Specifically, for plutonium monocarbide, it has been argued that the 5f electrons are less correlated than they are in both plutonium oxides and nitrides [11]. No theoretical method is currently able to predict with certainty the nature of the 5f electrons in these systems but model comparisons with experimental data provide clues.
In Figure 2 we show the (DFT+SO+OP) total electronic density of states (e-DOS) for PuC. The PuC e-DOS is dominated by 5f states and shows substantial weight and even a peak in the vicinity of the Fermi level (highest occupied energy state that is shifted to zero). This behavior is quantitatively comparable to that of both the α and δ phases of plutonium [5] because these phases also display significant occupation of 5f states at the Fermi level. The degree of 5f itineracy can, however, only be judged knowing how well theory reproduces measured photoelectron spectra.
The photoelectron spectra for sub-stoichiometric PuC0.85 has been measured by Gouder et al. [49] and in Figure 3 we compare that result with our e-DOS for PuC that has been broadened and convoluted to simulate instrumental resolution and photon lifetimes [50]. Gouder et al. utilized two different photon energies; HeI at 21.22 eV and HeII at 40.81 eV. Here we show the latter (40.81 eV) because it enhances the 5f photoexcitation cross section that makes the 5f contribution dominate the The PuC e-DOS is dominated by 5f states and shows substantial weight and even a peak in the vicinity of the Fermi level (highest occupied energy state that is shifted to zero). This behavior is quantitatively comparable to that of both the α and δ phases of plutonium [5] because these phases also display significant occupation of 5f states at the Fermi level. The degree of 5f itineracy can, however, only be judged knowing how well theory reproduces measured photoelectron spectra.
The photoelectron spectra for sub-stoichiometric PuC 0.85 has been measured by Gouder et al. [49] and in Figure 3 we compare that result with our e-DOS for PuC that has been broadened and convoluted to simulate instrumental resolution and photon lifetimes [50]. Gouder et al. utilized two different photon energies; HeI at 21.22 eV and HeII at 40.81 eV. Here we show the latter (40.81 eV) because it enhances the 5f photoexcitation cross section that makes the 5f contribution dominate the spectra [49]. We also include, for comparison, the corresponding simulated photoemission for δ-plutonium [50].
Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 16 spectra [49]. We also include, for comparison, the corresponding simulated photoemission for δplutonium [50].  [49] and the simulation for -Pu is reproduced from reference [50]. The curves are arbitrarily shifted relative to each other to make the plot more readable.
The experimental PuC0.85 photoelectron result actually resembles our simulation for PuC quite well. One would not necessarily expect a perfect agreement because the stoichiometries are different and the calculated spectral density do not include photoexcitation matrix elements. Nevertheless, the favorable comparison supports the interpretation that the 5f electrons, that dominate the spectra for the plutonium-carbide system, are indeed delocalized as treated by the theory. We further notice in Figure 3 that the result for -plutonium is quantitatively comparable to PuC. The estimated 5f-band occupation is marginally lower in PuC (4.95) than in -plutonium (5.31) [51].
Plutonium monocarbide (PuC1−x) is known to be anti-ferromagnetic (AF) at temperatures below 100 K, but, to date, the size of the magnetic moments has not been accurately determined but estimated to be ~0.7 μB [52]. Previous DFT calculations within the local-density approximation [53] reproduce the AF configuration for PuC0.75 but for ideal-stoichiometry PuC, a ferromagnetic (FM) configuration was predicted with considerably larger magnetic moments ~2 μB. The present calculations for PuC, that are based on GGA and the aforementioned spin-orbit and orbitalpolarization correlations, predict the FM state but it is only weakly stable over the AF configuration (~1 mRyatom −1 ). Both configurations have total magnetic moments that are smaller (~0.1 μB) because of the effective compensation between the spin and anti-parallel orbital components (both close to 3.7 μB in absolute magnitude). Orbital polarization tends to enhance the orbital moment in metallic plutonium systems [5] and ignoring this interaction in the model produces a much smaller orbital moment (~2.7 μB). The total moment from our best calculation is somewhat smaller than the experimental estimate but there are uncertainties in experimental and theoretical assumptions. For Figure 3. Total electronic density of states, from our best theory (DFT+SO+OP) for PuC (red) and δ-Pu (blue) that have been convoluted due to instrumental resolution and photon lifetime broadening. The experimental data points for PuC 0.85 are from photoemission (40.81 eV photon energies) by Gouder et al. [49] and the simulation for δ-Pu is reproduced from reference [50]. The curves are arbitrarily shifted relative to each other to make the plot more readable.
The experimental PuC 0.85 photoelectron result actually resembles our simulation for PuC quite well. One would not necessarily expect a perfect agreement because the stoichiometries are different and the calculated spectral density do not include photoexcitation matrix elements. Nevertheless, the favorable comparison supports the interpretation that the 5f electrons, that dominate the spectra for the plutonium-carbide system, are indeed delocalized as treated by the theory. We further notice in Figure 3 that the result for δ-plutonium is quantitatively comparable to PuC. The estimated 5f-band occupation is marginally lower in PuC (4.95) than in δ-plutonium (5.31) [51].
Plutonium monocarbide (PuC 1−x ) is known to be anti-ferromagnetic (AF) at temperatures below 100 K, but, to date, the size of the magnetic moments has not been accurately determined but estimated to be~0.7 µ B [52]. Previous DFT calculations within the local-density approximation [53] reproduce the AF configuration for PuC 0.75 but for ideal-stoichiometry PuC, a ferromagnetic (FM) configuration was predicted with considerably larger magnetic moments~2 µ B . The present calculations for PuC, that are based on GGA and the aforementioned spin-orbit and orbital-polarization correlations, predict the FM state but it is only weakly stable over the AF configuration (~1 mRyatom −1 ). Both configurations have total magnetic moments that are smaller (~0.1 µ B ) because of the effective compensation between the spin and anti-parallel orbital components (both close to 3.7 µ B in absolute magnitude). Orbital polarization tends to enhance the orbital moment in metallic plutonium systems [5] and ignoring this interaction in the model produces a much smaller orbital moment (~2.7 µ B ). The total moment from our best calculation is somewhat smaller than the experimental estimate but there are uncertainties in experimental and theoretical assumptions. For example, carbon vacancies enhance the total magnetic moment (see Section 4). The magnetic moments and the bonding properties of the AF and FM configurations are nearly identical for PuC so that the atomic volume and bulk modulus are essentially the same. Because the energetics and bonding between these two magnetic states are so similar and since the magnetic ordering occurs below 100 K, we only consider the FM state as we proceed and focus on thermodynamics at high temperatures.

Ground-State Properties and Thermodynamics
Here, we discuss our predicted ground-state properties, including the elastic constants, and our lattice-dynamics results. In Table 1 we show atomic volumes and elastic constants for PuC from the full theoretical treatment (DFT+SO+OP) and for the simpler theory that excludes spin-orbit coupling and orbital polarization (DFT). First, we recognize that the two theoretical approaches overestimate the atomic volume relative to both Zachariasen's and Storms's experimental values (14.8 and 15.4 Å 3 , respectively) [54,55]. The same trend, with a predicted too-large volume for PuC, was observed in a hybrid-DFT study [11] where they argued that the discrepancy was due to an overestimation of 5f localization. We believe the reason is that the approximation for the electron exchange and correlation, GGA, is not ideal for PuC. It is known that GGA is accurate for the light actinide metals and Pu [56,57], but for carbon (graphite), the local-density approximation (LDA) is better for the structural properties [58]. Consequently, neither LDA nor GGA is perfectly suitable for PuC and, therefore, the GGA atomic volume is imperfect.
Next, we compare our elastic constants in Table 1 with results obtained from Born-Mayer and Coulomb model potentials that Vermal et al. reported for PuC [59]. Their zero-temperature moduli are considerably smaller (C 11 = 63.9, C 12 = 27.2, C 44 = 27.2 GPa) than ours.
We expect that our elastic constants for PuC in Table 1 are reasonable because they are rather accurately calculated, using the same theoretical framework, for the known phases of plutonium metal [60]. The Poisson's ratio, ν, in Table 1 is obtained from the standard formula for cubic crystals: Unfortunately, there are no published experimental data on the single-crystal elastic constants (C ij ) that could help distinguish between the present and previous [59] modeling for PuC. For the isothermal bulk modulus (B), however, an experimental value has been reported (118 GPa) [61] that is far closer to our best theory (125 GPa) than that of the model potential (39.4 GPa) [59]. Table 1. Ground-state properties obtained from the full theoretical treatment (DFT+SO+OP) and for a treatment that ignores spin-orbit coupling and orbital polarization (DFT). The atomic volume (V) is given in units of Å 3 , while the bulk (B), shear (G), and elastic (C ij ) moduli are given in units of GPa. The shear modulus is a Voigt-Reuss-Hill average and the Poisson's ratio, ν, is obtained from B and G (see main text). In regard to our first-principles results in Table 1, we find some differences between the full theory (DFT+SO+OP) and the simpler scalar-relativistic approach (DFT). The former produces a somewhat larger atomic volume and a correspondingly softer bulk modulus, while the shear modulus is about the same. Most notably is the fact that the tetragonal shear constant, C = (C 11 − C 12 )/2, is very small for the scalar-relativistic approximation (27 vs. 70 GPa). The small C indicates that the crystal is close to a mechanical instability and this is linked to softer phonons for the related phonon branches. Softer phonon modes imply more entropy and a greater phonon contribution to the free energy, and that is also observed in our calculated free energies shown in Figure 1.

Method
Moving on to thermodynamics and our lattice-dynamics predictions, we show in Figure 4 our phonon density of states (p-DOS) at 1200 K. The p-DOS clearly suggests that including spin-orbit coupling and orbital polarization (DFT+SO+OP) in the electronic structure stiffens the lattice dynamics resulting in weight of the phonon modes shifting to higher energies.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 16 The consequence of the results shown in Figures 1 and 4 is that when the temperature dependent part of the electronic structure (Fermi-Dirac distribution and electronic entropy) is added to the Gibbs energies of the two models, the "DFT" is below that of the full "DFT+SO+OP" theory because the magnitude of the lattice-vibration contribution is greater.
In Figure 5 we show the Gibbs energies obtained from adding the electronic and lattice-vibration contributions at constant volumes (16.2 Å 3 and 15.6 Å 3 ) together with our CALPHAD Gibbs energy for PuC. The CALPHAD data begin at 298 K, but an extrapolation to zero K and a shift so that the Gibbs energy vanish at zero temperature makes the comparison appropriate. The lesser (DFT) theory lies below the full (DFT+SO+OP) treatment and appears to be in slightly better agreement with the CALPHAD result.  The consequence of the results shown in Figures 1 and 4 is that when the temperature dependent part of the electronic structure (Fermi-Dirac distribution and electronic entropy) is added to the Gibbs energies of the two models, the "DFT" is below that of the full "DFT+SO+OP" theory because the magnitude of the lattice-vibration contribution is greater.
In Figure 5 we show the Gibbs energies obtained from adding the electronic and lattice-vibration contributions at constant volumes (16.2 Å 3 and 15.6 Å 3 ) together with our CALPHAD Gibbs energy for PuC. The CALPHAD data begin at 298 K, but an extrapolation to zero K and a shift so that the Gibbs energy vanish at zero temperature makes the comparison appropriate. The lesser (DFT) theory lies below the full (DFT+SO+OP) treatment and appears to be in slightly better agreement with the CALPHAD result. The consequence of the results shown in Figures 1 and 4 is that when the temperature dependent part of the electronic structure (Fermi-Dirac distribution and electronic entropy) is added to the Gibbs energies of the two models, the "DFT" is below that of the full "DFT+SO+OP" theory because the magnitude of the lattice-vibration contribution is greater.
In Figure 5 we show the Gibbs energies obtained from adding the electronic and lattice-vibration contributions at constant volumes (16.2 Å 3 and 15.6 Å 3 ) together with our CALPHAD Gibbs energy for PuC. The CALPHAD data begin at 298 K, but an extrapolation to zero K and a shift so that the Gibbs energy vanish at zero temperature makes the comparison appropriate. The lesser (DFT) theory lies below the full (DFT+SO+OP) treatment and appears to be in slightly better agreement with the CALPHAD result.  This is an accidental better agreement because important contributions to the Gibbs energy are still missing in the first-principles modeling. First, removing the constraint of a fixed atomic volume, i.e., by accounting for thermal volume expansion, lowers the Gibbs energies a substantial amount. It comes with a significant computational cost to compute this contribution, however, because SCAILD needs to be repeated for two or more volumes (for each temperature) so that the equilibrium volume can be determined. In addition, since we are studying magnetically disordered PuC (over 100 K [52]), there is a simple magnetic-entropy contribution due to magnetic disorder [62,63] that we shall include: in this equation, k B is the Boltzmann constant and µ the total (spin and orbital) magnetic moment. Because orbital polarization, that enhances the magnitude of the orbital moment [5], is addressed in the model the total magnetic moment is small (~0.1 µ B ) as is the F mag contribution. As mentioned, the small total moment is due to a near complete compensation between the spin and orbital components. Lastly, there is an electron-phonon-coupling term that we are not considering for the free energy. Accounting for all energy excitations from electrons and phonons and their distributions in a universal fashion is difficult and to our knowledge there is no efficient procedure to accurately determine this quantity. When we add the missing terms (except the electron-phonon term) to the Gibbs energy, the scalar-relativistic (DFT) energy is significantly below and incompatible with the CALPHAD results (not shown). The full theoretical treatment (DFT+SO+OP), on the other hand, produces energies that are above but near that obtained from CALPHAD. In Figure 6 we show our best theory of the Gibbs energy together with CALPHAD for PuC. They are quite close, with the first-principles data slightly above CALPHAD, consistent with the fact that electron-phonon interaction is neglected in the ab initio model. This is an accidental better agreement because important contributions to the Gibbs energy are still missing in the first-principles modeling. First, removing the constraint of a fixed atomic volume, i.e., by accounting for thermal volume expansion, lowers the Gibbs energies a substantial amount. It comes with a significant computational cost to compute this contribution, however, because SCAILD needs to be repeated for two or more volumes (for each temperature) so that the equilibrium volume can be determined. In addition, since we are studying magnetically disordered PuC (over 100 K [52]), there is a simple magnetic-entropy contribution due to magnetic disorder [62,63] that we shall include: in this equation, kB is the Boltzmann constant and μ the total (spin and orbital) magnetic moment. Because orbital polarization, that enhances the magnitude of the orbital moment [5], is addressed in the model the total magnetic moment is small (~0.1 µ B) as is the Fmag contribution. As mentioned, the small total moment is due to a near complete compensation between the spin and orbital components. Lastly, there is an electron-phonon-coupling term that we are not considering for the free energy. Accounting for all energy excitations from electrons and phonons and their distributions in a universal fashion is difficult and to our knowledge there is no efficient procedure to accurately determine this quantity. When we add the missing terms (except the electron-phonon term) to the Gibbs energy, the scalar-relativistic (DFT) energy is significantly below and incompatible with the CALPHAD results (not shown). The full theoretical treatment (DFT+SO+OP), on the other hand, produces energies that are above but near that obtained from CALPHAD. In Figure 6 we show our best theory of the Gibbs energy together with CALPHAD for PuC. They are quite close, with the first-principles data slightly above CALPHAD, consistent with the fact that electron-phonon interaction is neglected in the ab initio model. Figure 6. Free energies from our best first-principles theory and CALPHAD (same as in Figure 5). The solid line is guide to the eye only.
In order to make direct connection with experimental data we also calculate the specific heat at constant pressure, Cp. We compare anharmonic theory (DFT+SO+OP combined with SCAILD) with Figure 6. Free energies from our best first-principles theory and CALPHAD (same as in Figure 5). The solid line is guide to the eye only.
In order to make direct connection with experimental data we also calculate the specific heat at constant pressure, C p . We compare anharmonic theory (DFT+SO+OP combined with SCAILD) with lower-level quasi-harmonic theory within the Debye-Grüneisen model. This comparison allows us to evaluate the possible importance of anharmonicity at higher temperatures. In addition, we calculate CALPHAD C p for ideal PuC and PuC 0.869 for a better understanding of the significance of sub-stoichiometry. Three sets of experimental data [64][65][66] are collected for comparison with theory. These data sets are shown in Figure 7 where we also plot the free-electron contribution obtained from the following simple Sommerfeld assumption [67]: here, D(E F ) is the e-DOS at the Fermi level that is predicted by the first-principles model (DFT+SO+OP).
There is an implicit temperature dependence in D(E F ) because the electronic structure is subject to a Fermi-Dirac thermal distribution that influences the e-DOS. From the slope of the free-electron specific heat we obtain the Sommerfeld coefficient 39.6 mJmol −1 K −2 . We did not find any published experimental PuC Sommerfeld coefficient but a linear interpolation of the specific heat between 14 K and 19 K [68] gives a value of 40.6 mJmol −1 K −2 for PuC 0.8 . The electronic contribution to the heat capacity is included in the anharmonic theory as well as the two examples of quasi-harmonic theory in Figure 7 and it is an essential component, particularly at higher temperatures.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 16 lower-level quasi-harmonic theory within the Debye-Grüneisen model. This comparison allows us to evaluate the possible importance of anharmonicity at higher temperatures. In addition, we calculate CALPHAD Cp for ideal PuC and PuC0.869 for a better understanding of the significance of sub-stoichiometry. Three sets of experimental data [64][65][66] are collected for comparison with theory. These data sets are shown in Figure 7 where we also plot the free-electron contribution obtained from the following simple Sommerfeld assumption [67]: here, D(EF) is the e-DOS at the Fermi level that is predicted by the first-principles model (DFT+SO+OP). There is an implicit temperature dependence in D(EF) because the electronic structure is subject to a Fermi-Dirac thermal distribution that influences the e-DOS. From the slope of the freeelectron specific heat we obtain the Sommerfeld coefficient 39.6 mJmol −1 K −2 . We did not find any published experimental PuC Sommerfeld coefficient but a linear interpolation of the specific heat between 14 K and 19 K [68] gives a value of 40.6 mJmol −1 K −2 for PuC0.8. The electronic contribution to the heat capacity is included in the anharmonic theory as well as the two examples of quasi-harmonic theory in Figure 7 and it is an essential component, particularly at higher temperatures Perhaps the most glaring aspect of Figure 7 is that the three experimental data sets deviate strongly above 1000 K. We mention that all points above 1300 K in the Kruger and Savage [64] data set are extrapolated by them. To a degree, the differences in the data sets may be explained by different stoichiometry of the samples, but two of samples have the same stoichiometry. Contradictory experimental data complicates the CALPHAD assessments and the corresponding results become more uncertain. In the CALPHAD assessment of the PuC1−x phase [16], more weight was given to the experimental data of Kruger and Savage [64] relative to that of references [65,66] where the sharp increase of Cp above 1000 K leads to erroneous stability of the solid in the liquid . Measured and calculated specific heats at constant pressure, C p . The solid square (red), diamond (blue), and circle (black) are experimental data sets from references [64][65][66]. Dashed lines refer to results from CALPHAD. The solid black, blue, and red lines without symbols refer to our best anharmonic theory and two parameterizations of Debye-Grüneisen quasi-harmonic theory. The solid black line with triangles shows the calculated free-electron contribution to C p .
Perhaps the most glaring aspect of Figure 7 is that the three experimental data sets deviate strongly above 1000 K. We mention that all points above 1300 K in the Kruger and Savage [64] data set are extrapolated by them. To a degree, the differences in the data sets may be explained by different stoichiometry of the samples, but two of samples have the same stoichiometry. Contradictory experimental data complicates the CALPHAD assessments and the corresponding results become more uncertain. In the CALPHAD assessment of the PuC 1−x phase [16], more weight was given to the experimental data of Kruger and Savage [64] relative to that of references [65,66] where the sharp increase of C p above 1000 K leads to erroneous stability of the solid in the liquid regime. The CALPHAD data show, consistent with experiments, that the stoichiometry relatively weakly influences C p . At 700 K the CALPHAD C p varies from 28.2 (PuC) to 26.8 (PuC 0.869 ) Jmol −1 K −1 , corresponding to a decrease of~5% (1.4 Jmol −1 K −1 ) due to sub-stoichiometry. This result agrees with the data from Kruger and Savage [64], which show a small difference in C p of 1.25 Jmol −1 K −1 (~4.5%) between PuC 0.869 and PuC at 700 K.
Our anharmonic theory agrees fairly well with all experimental and CALPHAD data up to 1000 K, but only with the Kruger and Savage [64] and the PuC-specific CALPHAD results above 1000 K. Not surprisingly, it agrees somewhat less favorably with CALPHAD for sub-stoichiometric PuC 0.869 . Furthermore, Figure 7 suggests that the quasi-harmonic treatments, both the (a) and (b) variants (see Computational methods, Section 2), deviate from CALPHAD and anharmonic theory for PuC above 1000 K. This behavior is consistent with our findings for uranium mononitride where the quasi-harmonic approach was shown to be increasingly inaccurate at temperatures above 1000 K [14].

Summary and Discussion
We calculate thermodynamical properties; lattice dynamics, free energies, and heat capacities for ideal-stoichiometry plutonium monocarbide. The highest level of theory includes spin-orbit coupling and orbital polarization for the DFT-GGA electronic structure. This approach assumes delocalized (band) 5f electrons on plutonium and a direct comparison with the experimental photoelectron spectroscopy confirms that this assumption is appropriate for PuC. Furthermore, this interpretation of the 5f electrons is consistent with the calculated Pu-C and Pu-Pu distances (2.52 and 3.57 Å, respectively). Namely, the Pu-Pu distance is close to the Hill limit (~3.4 Å) [69] that is an approximate criteria for 5f-band formation in the actinides.
The thermal properties are obtained from combining our advanced DFT to a temperature-dependent self-consistent phonon method that includes strong anharmonic lattice dynamics. We show that this level of electronic-structure theory and lattice dynamics are necessary to reproduce the CALPHAD Gibbs energy and produce consistent specific-heat data. The quasi-harmonic approaches, on the other hand, predict heat capacities that do not compare as well with CALPHAD and the most believable experimental C p (Kruger and Savage [64]). Still, the simplicity and efficiency of the Debye-Grüneisen technique are appealing for general thermodynamic modeling. However, we suggest that the best results are obtained if the inherent model-parameters are chosen so that experimental data (heat capacity) or the higher-order anharmonic theory (lattice-vibration energy) are reproduced as closely as possible. On the other hand, if the anharmonic effects are substantial, a good fit to either the heat capacity or the Gibbs energy may not be feasible.
The consistency between our first-principles modeling, CALPHAD, and experiments is further validated by comparing PuC formation enthalpies. The zero-temperature formation enthalpy for PuC is calculated, with our best theory (DFT+SO+OP), as the total-energy difference between the compound and its solid-form constituents, i.e., carbon (α-C, graphite) and α-plutonium; E(PuC)-E(α-Pu)-E(α-C) = −21 kJmol −1 . Here, the structure of all phases is carefully relaxed (not a trivial task for α-plutonium that is a complex monoclinic phase [5]). The first-principles formation enthalpy is in satisfactory agreement with CALPHAD data (−26.4 kJmol −1 ) and the experimental values, that range from −15.4 to −25.3 kJmol −1 [16,66]. Hence, there is favorable coherence in the thermodynamical properties between modeling and experiments.
The high-temperature upturn in the specific heat found in two experiments [65,66] could possibly be a magnetic effect but is more likely due to the defects of the material that are also affecting the comparisons with theory. Namely, the correlation between theoretical and experimental specific heats is obscured by the fact that plutonium monocarbide naturally forms with a deficiency of carbon atoms, i.e., PuC 1−x , where x is~0.1.
One can capture the effect of carbon deficiency in supercell calculations by removing carbon atoms and creating vacancies on the carbon sites. Specifically, by removing one carbon atom from a 16-atom supercell, we study the Pu 8 C 7 compound, i.e., PuC 1−x , where x = 0.125. The corresponding calculation is performed with VASP and the full DFT+SO+OP treatment, allowing for structural relaxation of the supercell with the vacancy. It turns out that the relaxation effects are small, but the atomic volume increases significantly accompanied by a softening of the bulk modulus and a small enhancement (~0.2 µ B ) of the total moment. Complementary to these supercell calculations, we conduct an EMTO-CPA investigation of the disordered PuC 1−x system where the number-of-atoms ratio, C/Pu = 0.90, and 10% of the carbon atoms are replaced by vacancies but without structural relaxation. In other words, the carbon atoms and vacancies are randomly distributed on the carbon-type sub-lattice with 10% probability of being a vacancy and 90% probability of being a carbon atom. Consistent with the VASP supercell results, our disorder EMTO-CPA model predicts an increase in atomic volume and a decrease of the bulk modulus when carbon is eliminated. This behavior mirrors a recent theoretical study [70] of the stoichiometry in PuC.
From the computed bonding energetics of the carbon-deficient system we apply the quasi-harmonic treatment [37] and compare that with an analogous calculation for ideal-stoichiometric PuC for VASP and EMTO (not shown). The difference in the specific heats, due to the sub-stoichiometry, proves to be very small for both methods and it cannot explain even the relatively weak sensitivity to the stoichiometry reflected in the experimental heat capacity.