A theoretical study for thorium monocarbide (ThC)

Journal of Nuclear Materials 429 (2012) 55–69

Contents lists available at SciVerse ScienceDirect

Journal of Nuclear Materials journal homepage: www.elsevier.com/locate/jnucmat

A theoretical study for thorium monocarbide (ThC) S. Aydin, A. Tatar, Y.O. Ciftci ⇑ Gazi University, Department of Physics, Teknikokullar 06500, Ankara, Turkey

h i g h l i g h t s " We focused on high pressure behavior of ThC. " ThC is metallic, and mechanically stable. " The obtained results agree with the other available values. " ThC is hard material, and hardness increases properly with pressure.

a r t i c l e

i n f o

Article history: Received 21 November 2011 Accepted 23 May 2012 Available online 30 May 2012

a b s t r a c t The structural, mechanical, electronic and thermodynamic properties of thorium monocarbide (ThC) with NaCl-type structure have been investigated by using first-principles plane wave density functional calculations with GGA, LDA and LDA + U functionals. It is shown that calculated equilibrium structural parameters of ThC are in agreement with the experimental results. It is seen from calculated single-crystal elastic constants that ThC with NaCl-type structure is mechanically stable. And from calculated density of states and band structure, it is observed that ThC is metallic. After the properties at 0 GPa are clarified, pressure dependency of the structural parameters, the elastic properties and related mechanical properties, density of states (DOS) and hardness are studied. Furthermore, the thermodynamic properties of ThC are obtained from the quasi-harmonic Debye model (QHM) over high pressure and temperature ranges for three functionals. The results are compared to each other, and the available experimental and theoretical data. Ó 2012 Published by Elsevier B.V.

1. Introduction Lanthanides and actinide compounds with the sodium chloridetype structure (B1) have been subject to theoretical [1] and experimental studies [2] for many years. The unusual properties of these compounds [3] and the role of f electrons in the nature of bonding lead to arise a particular interest. On one hand, Söderlind et al. studied some actinide materials and complexes by DFT, DFT + U and DFT + DMFT methods [4], Ducher et al. performed DFT + U calculations due to f-electron characteristics in uranium monocarbide [5], Wang et al. calculated electronic and mechanical properties, phonon dispersion curves of NpO2 with first-principles LDA + U and GGA + U functionals [6]. Among lanthanides and actinide compounds, actinide carbides have the smallest lattice constants, corrosion resistance and high melting point. Further, the binary compounds formed by thorium with light sp elements such as carbon and boron possess high density, good thermal conductivity and some other interesting physical properties. These properties ⇑ Corresponding author. Tel.: +90 312 2021266; fax: +90 312 2122279. E-mail address: [email protected] (Y.O. Ciftci). 0022-3115/$ - see front matter Ó 2012 Published by Elsevier B.V. http://dx.doi.org/10.1016/j.jnucmat.2012.05.038

make them attractive from a fundamental point of view as well as for a variety of technical applications, for example, alternative fertile materials used in nuclear breeder systems [7–10]. However, the Th–C system has two basic phases, mono-(ThC) and dicarbide (ThC2) [9]. The cubic (B1) thorium monocarbide has a wide region of nonstoichiometry in the carbon sublattice ThC1x (0 < x < 0.33), and dicarbide exists as three polymorphic modifications a, b, and c. High-pressure X-ray diffraction studies have been performed on thorium monocarbide (ThC) powder up to 36 GPa by Gerward et al. [11]. They found that no structural phase transition has been observed. On the other hand, Olsen et al. [12] studied on uranium and thorium compounds with rock salt-structure by using highpressure X-ray diffraction, and observed that there is no phase transition up to 40–45 GPa for ThC and ThN. The electronic structure and chemical bonding of ThC have been analyzed by employing the full-potential linearized augmented-plane-wave (APW) method within a generalized gradient approximation (GGA) for the exchange correlation potential [13] and full-potential all-electron density functional theory (DFT) method with mixed basis APW + lo(LAPW) method [14]. Lim and Scuseria [15] have also studied structural parameters, bulk moduli, density of states, and

56

S. Aydin et al. / Journal of Nuclear Materials 429 (2012) 55–69

band structures of thorium carbide by using Gaussian-typeorbitals and density functional theory, using periodic boundary conditions. The chemical bonding and the physico-chemical properties of binary Th–X (X = H, B, C, N, O, P, S, As, Se, Sb) and ternary thorium compounds (carbonitrides, thorates and silicates) are reviewed by Shein and Ivanovskii [16]. Moreover, the elastic properties of all the known cubic binary phases of thorium with nonmetals [17], and the effects of carbon non-stoichiometry on the structural and electronic properties of cubic ThC [18] have been studied by the same authors. Das et al. [19] have investigated the characteristics of the electronic structures of the carbides of the transition series Sc, Ti, V, Zr and Hf, the actinides U and Th using the tight-binding linear muffin-tin orbitals methods (TB-LMTO). They studied the density of states, charge distribution, energy bands in lowest energy configuration. In this work, we have aimed to provide some additional information to the existing data on the physical properties of thorium monocarbide by using the first-principles plane wave pseudopotential method. We have investigated structural, mechanical and electronic properties of ThC as a function of hydrostatic pressure, and also studied thermodynamic properties by means of the quasi harmonic Debye model. 2. Method of calculation

ð1Þ

where E(V) is the total energy of the unit cell, PV is the constant hydrostatic pressure condition, h(V) is the Debye temperature and Avib is the vibrational Helmholtz free energy, which can be written as [33–35]

Av ib ðh; TÞ ¼ nkT

   9h h þ 3 lnð1  eh=T Þ  D 8T T

ð2Þ

where n is the number of atoms per formula unit, and D(h/T) is the Debye integral. Debye temperature h is expressed as [33]

ð3Þ

where M is the molecular mass of unit cell and BS is the adiabatic bulk modulus, which is approximated given by the static compressibility [34]

BS  BðVÞ ¼ V

" # @ 2 EðVÞ

ð4Þ

@V 2

f(r) is given by [32,33]

8 " 9  3=2  3=2 #1 =1=3 < 2 1þr 1 1þr f ð rÞ ¼ 3 2 þ : ; 3 1  2r 3 1r

ð5Þ

where r is the Poisson ratio. By mean of Eqs. (2)–(5), Avib becomes a computable quantity, and thus the non-equilibrium Gibbs function G ðV; P; TÞ can be minimized with respect to volume V as

   @G ðV; P; TÞ ¼0 @V P;T

ð6Þ

The thermal equation-of-state (EOS) V(P, T) can be obtained by solving Eq. (6). The isothermal bulk modulus BT is given by [34]

BT ðP; TÞ ¼ V

Our calculations have been performed using the plane-wave pseudopotential approach to the density-functional theory (DFT) [20] by CASTEP code [21]. The electronic wave functions were obtained by using a density-mixing minimization method for the self consistent field (SCF) calculation, and the structures were relaxed by using the Broyden, Fletcher, Goldfarb and Shannon (BFGS) method [22]. For the exchange–correlation effects, Generalized Gradient Approximation (GGA-PBE) [23], local density approximation (LDACAPZ) [24] and LDA + U (U is Hubbard parameter) [25,26] were implemented. The Vanderbilt ultrasoft pseudopotentials [27] were used for interactions between electrons and core ions. The tolerances for geometry optimization were set as the difference in total energy being within 5  107 eV/atom, the maximum ionic Hellmann–Feynman force within 0.01 eV/Å, the maximum ionic displacement within 5  104 Å, and the maximum stress within 0.02 GPa. Pseudo-atomic calculations are performed for Th: 6s26p66d27s2 and C: 2s22p2. We used two different cut-off energies and k-points sets for geometry optimizations and electronic properties: 500 eV cut-off energy and 12  12  12 k-points for geometry optimizations, and 750 eV cut-off energy and 20  20  20 k-points for electronic properties such as band structure, density of states and Mulliken population analysis, because they are very sensitive to these parameters. Additionally, geometry optimizations and elastic properties are calculated on primitive cell which consists of two atoms, while electronic structure calculations are performed on unit cell of 8-atoms. The quasi-harmonic Debye model (QHM) [28–32] has been applied to calculate the thermodynamic properties of ThC compounds as implemented by GIBBS code [28]. In this model, the non-equilibrium Gibbs function G(V, P, T) can be written as

G ðV; P; TÞ ¼ EðVÞ þ PV þ Av ib ðhðVÞ; TÞ

rffiffiffiffiffi  h Bs 2 1=2 1=3 h ¼ ½6p V n f ðrÞ k M

" # @ 2 G ðV; P; TÞ @V 2

ð7Þ

P;T

The thermodynamic quantities, e.g. the heat capacities at constant volume (CV) and at constant pressure (CP), the thermal expansion coefficient (a), the Grüneisen parameter (c) and entropy (S) can be calculated by applying the following relations [34]:

  3h=T C V ¼ 3nk 4Dðh=TÞ  h=T e 1

ð8Þ

C P ¼ C V ð1 þ acTÞ

ð9Þ

  S ¼ nk 4Dðh=TÞ  3 lnð1  eh=T Þ

ð10Þ

a¼

cC V BT V

ð11Þ

c¼

d ln hðVÞ d ln V

ð12Þ

3. Results and discussion 3.1. Structural properties ThC has face-centered cubic lattice (space group Fm-3 m, No:225) in NaCl-type (B1) structure (see Fig. 1), there are two atoms in the primitive cell, their atomic positions are (0.0, 0.0, 0.0) for Th, and (0.5, 0.5, 0.5) for C, while there are eight atoms in the unit cell. We determine optimal value of Hubbard U parameter due to the experimental lattice constant (5.335 Å [9]), for this purpose, we have performed LDA + U calculations with different U values, and the results are presented in Fig. 2a with energy difference (DE) between LDA and LDA + U. It is observed that there is a good agreement between calculated and experimental lattice constant at U = 2.3 eV. Therefore, in all calculations of LDA + U within the work, this U value is used. Also, it is note that DE increases with U. Furthermore, the total energy as a function of the cell volume is calculated for different functionals, and is shown in Fig. 2b. It is shown that GGA energies are lower than those of LDA and LDA + U. Also, LDA energies are lower than those of LDA + U. It can be seen that addition of Hubbard parameter to the standard DFT calculations do not decrease energy of the system.

S. Aydin et al. / Journal of Nuclear Materials 429 (2012) 55–69

Fig. 1. The unit cell of thorium monocarbide (ThC) (blue spheres represent thorium atoms, and gray ones are carbon atoms). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

57

The lattice parameters optimized without any constraints, bulk moduli and Th–C bond lengths for different functionals are given in Table 1 with the available theoretical [14,15,20] and experimental data [9]. In comparison with experimental measurements, our calculated lattice constant (a) for LDA and GGA is different about 1.24% and 0.1%, respectively. And, GGA result is good agreement with experiment. Moreover, the bulk modulus and its pressure derivative have been calculated by fitting the first-principles energy–volume data to the third order Birch–Murnaghan equation of state at T = 0 K, the results for bulk modulus in cases of GGA, LDA and LDA + U are 134.43, 143.41 and 145.70 GPa, respectively, they are compatible with the available values. On the other hand, as a structural property, we have studied bond length of Th–C bond, which has significant effects on the properties such as melting point, elastic constants and hardness, are also listed in Table 1. Obviously, the our results are in good agreement with experimental data and other calculations [17,36]. To further investigate the pressure-induced structural changes of ThC, the variations of normalized lattice parameter (a/a0) and volume (V/V0) versus pressure are presented in Fig. 3 for different functionals with GIBBS results, where a0 and V0 are lattice param-

Fig. 2. (a) Determination of Hubbard U parameter with respect to experimental lattice constant, (b) total energy difference (DE, relative to GGA results) versus volume in ThC for GGA, LDA and LDA + U functionals.

58

S. Aydin et al. / Journal of Nuclear Materials 429 (2012) 55–69

Table 1 Calculated equilibrium lattice parameters (a0), volume (V0), bulk modulus (B0), its pressure derivative (B00 ), and length of Th–C bond for ThC compound with other theoretical and experimental data. Material

a (Å)

V (Å3/cell) B0 (GPa)

B00

Th–C (Å)

This study (GGA) LDA LDA + U (U = 2.3 eV) Exp.

5.341 5.269 5.336 5.335– 5.344a 5.3879e 5.3878g 5.261 5.34 5.363

152.375 146.281 151.910

134.33 143.41 145.70 109b,125c

3.071 3.019 3.217 3.1b

2.671 2.635 2.668 2.673d

131.89e

2.73e

2.650f

141.09 132.21 145.79

2.80 2.88 3.18

Theory LSDAh PBEh HSEh a b c d e f g h

Ref. [9]. Ref. [11]. Ref.[8]. Ref. [36]. Ref. [14]. Ref. [17]. Ref. [20]. Ref. [15].

eter and volume at zero pressure, respectively. It is observed from Fig. 3 that all ratios decrease with pressure and calculated curves agree with the experimental data. When pressure is applied up to 50 GPa, a/a0 and V/V0 values reduce by 8.0% and 22.0% for GGA, 7.9% and 21.8% for LDA, and 7.7% and 21.3% for LDA + U, respectively. It can be seen from Fig. 3 that normalized volume values are greater than the experimental values [11]. In order to characterize changing of a/a0 and V/V0 with pressure, we fit these ratios to a quadratic function, and results are given in Table 2. Especially, experimental V/V0 curve begins to fluctuation at 30–36 GPa, and LDA + U results have also some deviations around 40 GPa. These deviations can be interpreted as disordering in the structures or approaching to a phase transition region. This idea agrees with reported results in Ref. [12] that no transition is found for ThC up to 40–45 GPa. 3.2. Band structure and density of sates (DOS) We have computed the band structure of ThC along the high symmetry directions in the first Brillouin zone and the correspond-

Table 2 Fit results of a/a0 and V/V0 to second order polynomial for GGA, LDA and LDA + U functionals. Functional

Quadratic function

GGA

a/a0 = 0.9991  2.343  103 P + 1.587  105 P2 V/V0 = 0.9968 - 6.804  103P + 5.080  105P2 a/a0 = 0.9992  2.118  103P + 1.148  105P2 V/V0 = 0.9972 - 6.191  10-3P + 3.866  10-5P2 a/a0 = 0.9995  2.105  103 P + 1.221  105P2 V/V0 = 0.9981 - 6.147  103P + 4.029  105P2

LDA LDA + U

ing partial density of the states at 0 GPa (see Figs. 4 and 5). For GGA, LDA and LDA + U functionals, it is observed that band structures are very similar. We can see that valence and conduction bands overlap at the Fermi level for all functional cases, and as a result, there is no energy gap at Ef. This suggests that ThC is metallic, and this behavior is supported by discussion of metallic nature of ThC given in Ref. [15]. Fig. 5 presents the total (TDOS) and partial (PDOS) density of states corresponding to the band structures shown in Fig. 4 for ThC with GGA, LDA and LDA + U functionals. The PDOS are very useful fundamental tool to get information about hybridization and the orbital character of the states. It is seen that there are six main peaks (V1, V2, V3, V4, C1 and C2) in the total DOS curve of ThC. The lowest valence states (V1) are essentially dominated by C-s states, with minor presence of Th-d and Th-s states. Other valence states (V2, V3, V4) are essentially dominated by C-p and Th-d states, followed by Th-f states, as also can be seen in Ref. [15]. The conduction bands (C1, C2) consist essentially of Th-d and Th-f with a minor presence of C-p states. The disappearing of the energy gap in DOS confirms the metallic nature of ThC. Additionally, GGA results are higher than LDA and LDA + U results. Unlike GGA and LDA cases, a new peak for Th-d states occurs at 4 eV in LDA + U case, and f-states and then total DOS is shifted to higher energies by 1 eV. Furthermore, we can see that there is a hybridization between C-p and Th-d states near the Fermi level, and this interaction implies that the bonding between Th and C atoms contains a covalent component at least (this situation agrees with the previous results in Ref. [4,18]). Also, due to the difference in the electronegativity between the comprising elements, some ionic character can be expected. The effect of pressure on electronic property of ThC is studied by calculating the density of states (DOS) at several pressures.

Fig. 3. Variations of normalized lattice parameter (a/a0) and volume (V/V0) as a function of pressure.

S. Aydin et al. / Journal of Nuclear Materials 429 (2012) 55–69

59

Fig. 4. The electronic band structure of ThC for GGA, LDA and LDA + U functionals.

Fig. 5. Total and partial density of states in ThC at 0 GPa for GGA, LDA and LDA + U functionals, (I) for carbon atoms, (II) for thorium atoms and (III) for total. The vertical dashed line represents the Fermi level.

The computed TDOSs of ThC at 0, 20, 40 and 50 GPa pressures are given in Fig. 6a–g for GGA, LDA and LDA + U functionals. We found that DOS curves exhibit general features and they are very similar for all of functionals. As pressure increases, the peaks in the valence states shift to the lower energies, and the peaks in the conduction states shift to the higher energies, and the gap in the valence region (between 8 eV and 4 eV) decreases with the pressure. These changes indicating that hybridization energy is increasing with small amount under high pressure. Moreover, all

peak energies decreases with increasing pressure. In addition, the value of DOS at Fermi level, N(Ef), a parameter related closely to the superconductivity of a material, and it increases with increasing pressure. From Fig. 6d–g, it is seen that LDA and GGA results are very similar at different pressures, and TDOS for LDA + U is shifted to higher energies as pressure increases. In order to reveal the bonding nature, we have plotted the charge–density distribution maps in the (1 0 0) plane for ThC in Fig. 7a–b, and have performed Mulliken population analysis in

60

S. Aydin et al. / Journal of Nuclear Materials 429 (2012) 55–69

Fig. 6. (a) Calculated total density of states of ThC at the different pressures for GGA, LDA and LDA + U functionals. (b) Comparisons of total density of states for different functionals at each pressure.

Fig. 7. Charge density distribution map in (1 0 0) plane of ThC at (a) 0 GPa and (b) 50 GPa with GGA. Electron density is high in red colored regions and is low in blue colored regions. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 2 (only GGA results are reported because GGA, LDA and LDA + U results are very similar). It is well-known that symmetric electron density centered at the atoms indicates to ionic character, while electron density along with the bond indicates to covalent character. It can be seen that electron density are centered around Th and C atoms and it’s small fraction distributes along Th–C bond, indicating dominant ionic and weak covalent character of Th–C bond (see Fig. 7a). In Fig. 7b, when the pressure increases from 0 GPa to 50 GPa, electron density between Th and C atoms increases, therefore, weak covalent character is forced with pressure (this situation is supported by density of states), but ionic character is still dominant. From Table 3, it is observed that while the charges of carbon atoms are negative, those of thorium atoms are positive. In other

words, while carbon atoms have negative charges (behave as anion), thorium atoms have positive charges (behave as cation). Firstly, due to the presence of two charge of opposite sign in the structure, ThC has an ionic nature. Secondly, owing to the fact that amount of loss valence charge by thorium atoms are equal to amount of received charge by carbon atoms for each pressure, the bonding in ThC is dominantly originated from ionic bonding component, like NaCl (its ionicity is very high). Finally, a measure of ionic character may also be obtained from the effective ionic valence (see Table 3), which is defined as the difference between the formal ionic charge and the Mulliken charge on the anion species. A value of zero indicates a perfectly ionic bond, while values greater than zero indicate increasing levels of covalency (see Ref. [37] for more details). It is seen from Table 3, when pressure is

61

S. Aydin et al. / Journal of Nuclear Materials 429 (2012) 55–69 Table 3 Calculated Mulliken atomic charges of C and Th atoms at different pressures with GGA. Pressure

Carbon (anion)

0 5 10 15 20 25 30 35 40 45 50

Thorium (cation)

s

p

Total

Charge (e)

Effective charge

s

p

d

f

Total

Charge (e)

1.56 1.56 1.55 1.55 1.55 1.54 1.54 1.54 1.53 1.53 1.53

3.22 3.22 3.22 3.22 3.23 3.22 3.22 3.22 3.22 3.22 3.21

4.79 4.78 4.78 4.77 4.77 4.76 4.76 4.76 4.75 4.75 4.75

0.79 0.78 0.78 0.77 0.77 0.76 0.76 0.76 0.75 0.75 0.75

3.21 3.22 3.22 3.23 3.23 3.24 3.24 3.24 3.25 3.25 3.25

2.36 2.36 2.35 2.35 2.35 2.35 2.36 2.36 2.36 2.36 2.36

6.19 6.15 6.1 6.06 6.03 5.98 5.95 5.92 5.89 5.86 5.83

2.21 2.24 2.28 2.3 2.32 2.35 2.36 2.37 2.39 2.4 2.42

0.45 0.47 0.5 0.51 0.53 0.55 0.57 0.59 0.61 0.63 0.64

11.21 11.22 11.22 11.23 11.23 11.24 11.24 11.24 11.25 11.25 11.25

0.79 0.78 0.78 0.77 0.77 0.76 0.76 0.76 0.75 0.75 0.75

Table 4 Calculated Mulliken bond population (M) and bond length for Th–C bond at different pressures for GGA, LDA and LDA + U functionals. Pressure

0 5 10 15 20 25 30 35 40 45 50

GGA

LDA

LDA + U

M

d (Å)

M

d (Å)

M

d (Å)

Spin

0.66 0.63 0.61 0.58 0.57 0.54 0.52 0.50 0.48 0.46 0.44

2.671 2.640 2.604 2.582 2.563 2.535 2.519 2.503 2.488 2.474 2.458

0.62 0.59 0.57 0.54 0.51 0.49 0.46 0.44 0.41 0.39 0.36

2.635 2.604 2.579 2.555 2.532 2.511 2.493 2.476 2.459 2.443 2.427

0.64 0.62 0.60 0.58 0.56 0.54 0.51 0.49 0.47 0.45 0.42

2.668 2.640 2.612 2.589 2.568 2.548 2.529 2.512 2.495 2.479 2.463

0 0 0 0 0 0 0 0 0 0 0

Table 5 Fit results of Th–C bond length to second order polynomial for GGA, LDA and LDA + U functionals. Functional

Quadratic function

GGA LDA LDA + U

d = 2.668  6.258  103 P + 4.239  105 P2 d = 2.632  5.580  103 P + 3.024  105 P2 d = 2.666  5.451  103 P + 2.866  105 P2

On the other hand, Mulliken bond population (M) can be used to describe bonding and anti-bonding interactions, and further, to determine bonding character (covalent or ionic) between two atoms. The positive and negative values of the population indicate bonding and anti-bonding states, respectively. A value of zero means an ideal ionic bond, namely, M ? 0, bond iconicity increases. And, a higher positive value indicates a high degree of covalency in the bond [38,39]. Calculated Mulliken bond population of Th–C bond, which is one characteristic bond in the structure is given in Table 4 as a function of pressure for GGA, LDA and LDA + U functionals. From Table 4 for all functionals, when pressure is increased, M approaches to zero, and ionicity of Th–C bond is increased. Moreover, it is expected that when pressure is increased, length of Th–C bond is decreased, and at 50 GPa, it is length decreases 8.0% for GGA, 7.9% for LDA, and 7.7% for LDA + U. And, bond length is fitted to second order polynomial with R2 value of 0.999 in Table 5. Furthermore, from LDA + U calculations, total spin located at the bond is zero at each pressure, and this result agrees with individual atomic spin behavior above. Moreover, density of states at the Fermi level, N(Ef), is important parameter for electronic properties of the material. In addition to used for calculation of the different parameters in first-principles such as electron–phonon coupling constant, metallicity (fm) of a material can be computed as a function of N(Ef) [40]:

fm ¼ 0:026NðEf Þ=ne increased, effective ionic valence increases, this trend indicates that ionicity of the Th–C bonds decreases, while the covalency increases. And, the bonding behavior concluded from Mulliken atomic charges and effective ionic valence analysis is supported to electron density maps discussed above. Furthermore, it is observed for LDA + U calculations that thorium and carbon atoms have not a spin, and then ThC has non-magnetic character, this situation agree with the available literature [14].

where ne is valence electron density, and is calculated as ne = N/V, N is total valence electron number, and V is volume of the cell. It is possible that N(Ef) can be calculated as partial for each orbitals of the atoms by mean of partial density of states calculations. And also, valence electron density ne can be split to orbital channels by using calculated Mulliken atomic charges above (see Table 3). Then, metallicity can be computed as partial. The calculated partial N(Ef) and corresponding partial fm are tabulated as a function of pressure in

Table 6 Calculated partial density of states at Fermi level for different pressures with GGA. Pressure

0 5 10 15 20 25 30 35 40 45 50

Carbon

Thorium

ThC

s

p

Total

s

p

d

f

Total

s

p

d

f

Total

0.01 0.01 0.02 0.02 0.02 0.03 0.03 0.03 0.04 0.04 0.04

1.86 1.95 2.03 2.07 2.08 2.11 2.12 2.16 2.16 2.16 2.12

1.87 1.96 2.05 2.09 2.10 2.13 2.15 2.20 2.19 2.20 2.17

0.02 0.03 0.03 0.04 0.04 0.04 0.05 0.05 0.06 0.06 0.06

0.41 0.41 0.41 0.42 0.41 0.40 0.40 0.41 0.41 0.41 0.40

0.83 0.96 1.10 1.18 1.24 1.33 1.35 1.41 1.43 1.48 1.52

0.59 0.67 0.74 0.77 0.82 0.90 0.93 0.98 1.01 1.06 1.11

1.85 2.06 2.28 2.40 2.50 2.67 2.73 2.85 2.90 3.01 3.09

0.03 0.04 0.05 0.06 0.06 0.07 0.08 0.09 0.09 0.10 0.10

2.27 2.36 2.44 2.48 2.48 2.51 2.52 2.57 2.57 2.57 2.53

0.83 0.96 1.10 1.18 1.24 1.33 1.35 1.41 1.43 1.48 1.52

0.59 0.67 0.74 0.77 0.82 0.90 0.93 0.98 1.01 1.06 1.11

3.72 4.03 4.33 4.49 4.60 4.80 4.88 5.05 5.09 5.20 5.26

62

S. Aydin et al. / Journal of Nuclear Materials 429 (2012) 55–69

Table 7 Calculated partial and total metallicity of ThC at the different pressures with GGA. Pressure

0 5 10 15 20 25 30 35 40 45 50

Carbon

Thorium

ThC

s

p

Total

s

p

d

f

Total

s

p

d

f

Total

0.01 0.01 0.01 0.01 0.01 0.01 0.02 0.02 0.02 0.02 0.02

0.57 0.58 0.58 0.57 0.56 0.55 0.55 0.55 0.54 0.53 0.51

0.39 0.39 0.39 0.39 0.39 0.38 0.38 0.38 0.37 0.36 0.35

0.01 0.01 0.01 0.01 0.01 0.02 0.02 0.02 0.02 0.02 0.02

0.07 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.05

0.37 0.41 0.44 0.46 0.47 0.48 0.48 0.48 0.48 0.48 0.49

1.30 1.35 1.35 1.36 1.35 1.38 1.35 1.36 1.32 1.33 1.34

0.16 0.18 0.19 0.19 0.19 0.20 0.20 0.21 0.21 0.21 0.21

0.01 0.01 0.01 0.01 0.01 0.02 0.02 0.02 0.02 0.02 0.02

0.24 0.24 0.24 0.24 0.23 0.23 0.23 0.23 0.23 0.22 0.22

0.37 0.41 0.44 0.46 0.47 0.48 0.48 0.48 0.48 0.48 0.49

1.30 1.35 1.35 1.36 1.35 1.38 1.35 1.36 1.32 1.33 1.34

0.23 0.24 0.25 0.25 0.25 0.25 0.25 0.26 0.26 0.26 0.25

Fig. 8. The variation of N(Ef) and metallicity with pressure for (a) GGA, (b) LDA and (c) LDA + U functionals.

Table 6 and Table 7, respectively (only GGA results are reported, because GGA, LDA and LDA + U results are very similar). It is seen from Table 6, the total N(Ef) of thorium atoms are higher than that of carbon atoms, except for 0 GPa with a small difference. When only carbon atoms are considered, p-orbitals are dominant on the N(Ef), while p- and f-orbitals are dominant in case of thorium atoms. However, when ThC is considered as a total, the highest contribution to N(Ef) comes from p-orbitals, and after d-orbitals. From Table 7, metallicity of ThC is originated from f-orbitals, and after d-orbitals. In order to present the variation of N(Ef) and metallicity with pressure, our calculated results are shown in Fig. 8. It is seen that N(Ef) and metallicity exhibit the same behavior for all functionals. N(Ef) increases with increasing pressure. Metallicity increases up to 40 GPa for GGA and LDA + U, but 30 GPa for LDA and after, it decreases. This trend can be interpreted as an indicator for a possible phase transition, because it is expected that a possible phase transition can change actual electronic properties of the material.

lengths and populations. The bond systems which possess negative population value are not considered in the hardness calculation, because negative population values mean anti-bonding or broken bonds [43]. For a binary bond system formed by i and j atoms, bond strength (Sij) is defined as,

Sij ¼

pffiffiffiffiffiffiffi ei ej =ðni nj dij Þ

3.3. Hardness calculation of ThC Hardness can be accurately predicted based on proper first principles results and semi-emprical data fitting. In this study, we used the method which was improved by Simunek [41,42]. A brief review of Simunek’s method is given below. Firstly, all bonds in the crystal structure are split to binary bond systems, which are formed by two atom, and in order to obtain number of different systems, these binary bond systems are grouped with respect to

Fig. 9. The variation of hardness with pressure for ThC.

ð13Þ

63

S. Aydin et al. / Journal of Nuclear Materials 429 (2012) 55–69 Table 8 Fit results of hardness to second order polynomial for GGA, LDA and LDA + U functionals. Functional

Quadratic function

GGA LDA LDA + U

H (GPa) = 5.527 + 5.555  102 P  2.564  104 P2 H (GPa) = 5.803 + 5.514  102 P  1.865  104 P2 H (GPa) = 5.506 + 5.368  102 P  2.517  104 P2

where dij is bond length, ni and nj are coordination numbers of i and j atoms, respectively. ei and ej are reference energies of i and j atoms, respectively, and reference energy of any atom is calculated with ea = Za/Ra, where Za is valence electron number, Ra is atomic radius. After bond strengths of all binary bond systems are determined, hardness of the crystal is calculated as,

H¼

C

X

" n

n Y

#1=n Nij Sij

erfe

ð14Þ

i;j¼1

where n is the number of different binary bond systems in the structure. X, Nij and fe are volume of unit-cell, number of bond each

binary bond system, and iconicity factor, respectively. C and r are the special constants in the method. In this study, we used C = 1450 and r = 2.8 [42], and took the atomic radii for different elements from Pearson text book [44]: atomic radius values of Th, and C elements are 1.798 and 0.916 Å, respectively. The calculated hardness of ThC is 5.5 GPa for GGA, 5.8 GPa for LDA, and 5.5 GPa for LDA + U at 0 GPa and increases properly with pressure (see Fig. 9). Thus, according to hardness scale, ThC is not superhard material, only hard material, and it is hardness is comparable to InP(5.1), GaSb(5.6) and AlSb(4.6) [41]. In addition, hardness can be expressed as a function of pressure by mean of fitting procedure (R2 value is 0.997) in Table 8. 3.4. Elastic properties The elastic constants of solids provide a link between the mechanical and dynamical behavior of crystals and give important information concerning the nature of the forces operating in solids. In particular, they provide information on the stability and stiffness of materials. Hence, to study the stability of ThC in B1 structure, we have calculated the second-order elastic constants (cij) at equilib-

Table 9 Elastic constants Cij (in GPa), Young’s modulus E (in GPa), Zener anisotropy factor A, Poisson’s ratio v and B/G ratio for ThC in B1 structure at different pressures for GGA, LDA and LDA + U functionals (p.cell and u.cell are abbreviations for primitive cell and unit cell, respectively). P (GPa) 0

GGA (p.cell) GGA (u.cell) LDA (p.cell) LDA (u.cell) LDA + U

5

GGA LDA LDA + U GGA (p.cell) GGA (u.cell) LDA (p.cell) LDA (u.cell) LDA + U GGA LDA LDA + U GGA (p.cell) GGA (u.cell) LDA (p.cell) LDA (u.cell) LDA + U GGA LDA LDA + U GGA (p.cell) GGA (u.cell) LDA (p.cell) LDA (u.cell) LDA + U GGA LDA LDA + U GGA (p.cell) GGA (u.cell) LDA (p.cell) LDA (u.cell) LDA + U GGA LDA LDA + U GGA (p.cell) GGA (u.cell) LDA (p.cell) LDA (u.cell) LDA + U

10

15

20

25

30

35

40

45

50

a

Ref. [17].

c11

c44

c12

B

G

E

B/G

v

A

192.6 276.4 244.6 240.9 214.8 252.2a 264.2 271.3 240.2 300.4 301.7 302.4 300.6 262.2 344.6 297.0 292.5 414.0 364.7 308.3 304.5 304.7 381.8 363.8 294.1 419.6 380.2 409.0 403.7 314.2 449.9 424.7 363.3 490.5 336.8 405.9 405.5 405.4 476.7 422.8 425.1 346.7 382.8 468.5 470.2 411.1

92.1 87.2 79.1 78.1 81.3 60.2a 99.0 76.0 79.6 82.9 82.6 73.9 73.0 76.6 88.8 70.9 74.8 111.3 91.4 68.0 66.8 72.7 79.9 67.4 69.9 102.8 90.5 64.9 62.0 67.8 88.8 59.5 67.4 87.1 38.6 52.0 52.6 65.4 78.9 46.3 60.6 37.4 37.6 44.6 47.7 53.9

100.7 99.1 98.3 96.4 87.9 96.3a 126.9 103.9 96.6 116.5 116.6 110.6 111.7 102.4 131.4 124.4 107.6 171.8 144.1 136.5 134.6 116.3 162.0 142.8 128.3 184.5 176.0 147.5 145.6 137.7 179.1 153.9 142.8 184.7 126.1 160.7 161.4 146.9 181.0 165.1 152.7 145.3 153.1 172.9 171.0 159.3

131.3 158.2 147.1 144.5 130.2 148.3a 172.7 159.7 144.5 177.8 178.3 174.6 174.7 155.6 202.5 181.9 169.2 252.6 217.7 193.7 191.3 179.1 235.3 216.5 183.6 262.9 244.1 234.7 231.6 196.5 269.4 244.2 216.3 286.7 196.3 242.4 242.8 233.1 279.6 251.0 243.5 212.4 229.7 271.4 270.7 243.2

69.7 87.8 76.7 75.7 73.6 66.7a 85.5 79.0 76.4 86.4 86.5 82.1 80.9 77.9 95.5 76.7 81.4 115.1 98.5 74.7 73.5 80.7 90.8 82.2 74.8 108.4 94.9 86.3 83.5 75.4 105.2 83.3 82.2 109.3 58.5 73.9 74.2 86.2 101.8 70.8 84.4 56.3 59.9 73.9 77.0 76.3

177.7 222.2 195.9 193.4 185.9 174.1a 220.2 203.4 194.9 223.2 223.3 212.8 210.3 200.2 247.6 201.7 210.5 299.7 256.8 198.5 195.5 210.4 241.3 219.0 197.6 286.0 252.1 230.6 223.7 200.5 279.4 224.4 218.8 290.9 159.6 201.2 202.1 230.2 272.2 194.1 227.0 155.2 165.4 203.2 211.0 207.3

1.88 1.80 1.92 1.91 1.77

0.27 0.27 0.28 0.28 0.26 0.30a 0.29 0.29 0.28 0.29 0.29 0.30 0.30 0.29 0.30 0.32 0.29 0.30 0.30 0.33 0.33 0.30 0.33 0.33 0.32 0.32 0.33 0.34 0.34 0.33 0.33 0.35 0.33 0.33 0.36 0.36 0.36 0.34 0.34 0.37 0.34 0.38 0.38 0.38 0.37 0.36

2.00 0.98 1.08 1.08 1.28 0.77a 1.44 0.91 1.11 0.90 0.89 0.77 0.77 0.96 0.83 0.82 0.81 0.92 0.83 0.79 0.79 0.77 0.73 0.61 0.84 0.87 0.89 0.50 0.48 0.77 0.66 0.44 0.61 0.57 0.37 0.42 0.43 0.51 0.53 0.36 0.45 0.37 0.33 0.30 0.32 0.43

2.02 2.02 1.89 2.06 2.06 2.13 2.16 2.00 2.12 2.37 2.08 2.19 2.21 2.59 2.60 2.22 2.59 2.63 2.45 2.42 2.57 2.72 2.77 2.61 2.56 2.93 2.63 2.62 3.36 3.28 3.27 2.70 2.75 3.55 2.89 3.77 3.83 3.67 3.52 3.19

64

S. Aydin et al. / Journal of Nuclear Materials 429 (2012) 55–69

Table 10 Calculated density, q (g/cm3), longitudinal, transverse and average elastic wave velocities (m/s) and Debye temperature (K) for ThC as a function of pressure for GGA, LDA and LDA + U functionals. P 0

5

10

15

20

25

30

35

40

45

50

GGA LDA LDA + U GGA LDA LDA + U GGA LDA LDA + U GGA LDA LDA + U GGA LDA LDA + U GGA LDA LDA + U GGA LDA LDA + U GGA LDA LDA + U GGA LDA LDA + U GGA LDA LDA + U GGA LDA LDA + U

q (g/cm3)

vt (m/s)

vl (m/s)

vm (m/s)

hD (K)

10.638 11.082 10.671 11.015 11.473 11.018 11.470 11.818 11.374 11.776 12.152 11.680 12.038 12.482 11.968 12.437 12.792 12.249 12.671 13.079 12.523 12.914 13.353 13.050 13.149 13.626 12.789 13.381 13.903 13.308 13.653 14.168 13.555

2559.8 2630.1 2627.0 2786.1 2623.7 2633.4 2745.2 2635.0 2616.7 2847.9 2512.3 2640.4 3092.0 2445.9 2596.3 2701.8 2535.5 2471.6 2925.3 2568.2 2453.3 2854.8 2497.6 2509.4 2883.0 2328.5 2596.2 2757.5 2256.6 2518.2 2031.1 2283.8 2372.7

4591.6 4743.1 4626.1 5101.5 4806.4 4728.5 5054.8 4901.9 4776.4 5292.4 4835.7 4877.1 5807.5 4847.5 4894.3 5352.7 5049.5 4809.8 5670.5 5171.0 4870.0 5632.4 5157.8 4997.0 5734.4 5002.1 5216.4 5570.7 4984.4 5172.4 4589.0 5109.9 5044.6

2850.2 2929.7 2920.6 3107.1 2926.1 2932.4 3062.8 2942.1 2917.5 3179.5 2811.7 2946.6 3454.7 2742.3 2901.6 3029.2 2843.7 2768.0 3275.3 2882.2 2750.9 3199.9 2806.9 2814.4 3233.1 2622.1 2913.3 3095.2 2544.4 2829.2 2292.4 2576.5 2670.6

458.3 477.5 470.1 505.4 482.5 477.0 505.0 489.9 479.7 528.8 472.6 488.7 578.8 465.0 485.2 513.1 486.2 466.5 558.2 496.4 467.0 548.8 486.8 484.4 557.9 457.9 498.1 537.2 447.3 490.1 400.5 455.8 465.5

rium lattice parameter by using the ‘‘stress–strain’’ technique. Other hand, the physical properties such as bulk (B), shear (G) and Young modulus (E), and Possion’s ratio can be written as functions of the elastic constants within Voigt–Reuss–Hill approximation [45]:

BV ¼ BR ¼ ðc11 þ c12 Þ=3

ð15Þ

GV ¼ ðc11  c12 þ 3c44 Þ=5;

and GR

¼ 5ðc11  c12 Þc44 =½4c44 þ 3ðc11  c12 Þ

ð16Þ

B ¼ ðBV þ BR Þ=2; G ¼ ðGV þ GR Þ=2

ð17Þ

E ¼ 9BG=ðG þ 3BÞ

ð18Þ

v ¼ ð3B  2GÞ=½2ð3B þ GÞ

ð19Þ

and Zener anisotropy factor (A) [46,47],

A ¼ 2c44 =ðc11  c12 Þ

ð20Þ

V subscript stands for Voigt, and is corresponding to the upper bound of the modulus, while R is Reuss, and is corresponding to the lower bound of the modulus. For cubic systems, there are three independent elastic constants (c11, c12 and c44), they have to satisfy following conditions for the mechanical stability:

c11  c12 > 0;

c11 > 0;

c44 > 0

ð21Þ

Calculated elastic constants and related physical properties are listed in Table 9 for GGA, LDA and LDA + U functionals. As shown in Table 6, calculated elastic constants satisfy all these stability criteria for each pressure. Therefore, we can say that this compound is elastically/mechanically stable in B1 structure over 0–40 GPa pressure range. At 0 GPa, c11 for LDA, c44 and c12 for GGA are greater than the others. c11 and c12 increases with pressure, while c44 decreases. Additionally, it is observed that GGA results dramatically decrease after 40 GPa. This situation can be originated from a possible phase transition. At this point, in order to remove possible suspicions, and to force this suggestion, we performed the calculations with GGA-PW91 at each pressure. It is seen that GGA-PW91 results support the GGA-PBE results. At the same time, this idea is supported by metallicity discussion above. At 40 GPa, c11 changes 154% for GGA, 66% for LDA, and 89% for LDA + U relative to 0 GPa values. Unlike LDA and LDA + U, it is values decreases by 30% at 50 GPa for GGA relative to 40 GPa. c12 changes 84% for GGA, 63% for LDA, and 67% for LDA + U at 40 GPa. Unlike LDA and LDA + U, it is values decreases by 27% at 50 GPa for GGA relative to 40 GPa. And, c44 changes 5% for GGA, 34% for LDA, and 20% for LDA + U at 40 GPa. Unlike LDA and LDA + U, it is values decreases by 130% at 50 GPa for GGA relative to 40 GPa. As a general result, the pressure sensitivity of c11 is higher than that of c12 and c44.

Fig. 10. (a) The variation of volume with pressure, and (b) with temperature for various isotherms and isobars.

S. Aydin et al. / Journal of Nuclear Materials 429 (2012) 55–69

However, c11 represents the elasticity in length, and a longitudinal strain produces a change in c11. c12 and c44 are related to the elasticity in shape, which are shear constants, and a transverse strain causes to change in shape. Moreover, the bulk modulus is a measure of resistance to volume change by applied pressure, whereas the shear modulus is a measure of resistance to reversible deformations upon the shear stress [48]. And, Young modulus E (the ratio between stress and strain) is required to provide information about the measure of the stiffness of the solids. It is shown from Table 9 that the mechanical properties such as bulk, shear and Young modulus increase with pressure. At 40 GPa, the values of bulk, shear and Young modulus are 286.7, 109.3 and 290.9 GPa, respectively, the deviations of them are 54%, 36% and 39% for GGA. It is seen that B exhibits larger deviation than the others, and pressure dependency is very high. Further, the ratio of B/G is frequently used to indicate the ductility of the compound. B/G is smaller than 1.75 for the brittle compound (0.8 for diamond), and is greater than

65

1.75 for the metallic compound (2.74 for Al) [49]. The calculated values of the B/G in Table 9 are higher than 1.75, and increase with increasing pressure which means that pressure improve ductility of the material. On the other hand, the Zener anisotropy factor A is a measure of the degree of elastic anisotropy in solids. The A takes the value of 1 for a completely isotropic material. If the value of A smaller or greater than unity it shows the degree of elastic anisotropy. Before our results are presented, we want to give a technical note about interesting observation in elastic constant calculations. As mentioned in method of calculation section of the work, we performed elastic constants calculations on primitive cell, but elastic constants calculated on primitive cell with GGA functional did not generate acceptable results for Zener anisotropy factor (2.00 at 0 GPa), thus in case of GGA functional, we recalculated elastic constants on unit cell, and obtained good results (0.98 at 0 GPa), this value agrees with the available data [14], and indicating that ThC can be characterized as an isotropic material. Furthermore, it is

Fig. 11. Relative volume ratio to GGA results versus pressure for (a) LDA and (b) LDA + U, versus temperature for (c) LDA and (d) LDA + U.

Fig. 12. (a) The variations of bulk modulus with pressure at different temperatures and (b) with temperature at different pressures for ThC.

66

S. Aydin et al. / Journal of Nuclear Materials 429 (2012) 55–69

seen that elastic constants calculated on the primitive and unit cell exhibit same trend under pressure in case of GGA, this situation/ observation do not affect our pressure investigation. However, our LDA results on the primitive cell were tested by those of calculated on the unit cell, and we showed that results are in good agreement, that is, we did not observed a case like in GGA. To present a general aspect, results for the unit cell are also listed in Table 9. In the wide range of applied pressure, our calculations show that the value of A decreases to 0.37 at 50 GPa, and anisotropy of thorium carbide increases. Additionally, the Poisson’s ratio is very important property for industrial applications since provides more information about the characteristics of the bonding forces than the elastic constants. The lower limit and upper limit of v are given 0.25 and 0.50 for central forces solids, respectively [50]. Calculated v value is equal to 0.27 for ThC at 0 GPa, and agrees with other the-

oretical values [17]. It is concluded that the interatomic forces in the ThC are central. From Table 9, Poisson’s ratio increases, and central nature of the forces is destroyed with applied pressure. The Debye temperature (hD) is known as an important fundamental parameter closely related to many physical properties such as specific heat and melting temperature. At low temperatures, the vibrational excitations arise more solely than acoustic vibrations. Hence, at low temperatures, the Debye temperature calculated from elastic constants is the same as that determined from specific heat measurements. We have calculated the Debye temperature, hD, from the elastic constants data using the average sound velocity, vm, by the following common relation [51],

hD ¼

  1=3  3n NA q h vm k 4p M

Fig. 13. Relative bulk modulus to GGA results versus pressure for (a) LDA and (b) LDA + U, versus temperature for (c) LDA and (d) LDA + U.

Fig. 14. (a) The variations of Grüneisen parameter with pressure at different temperatures and (b) with temperature at different pressures for ThC.

ð22Þ

S. Aydin et al. / Journal of Nuclear Materials 429 (2012) 55–69

where  h is Planck’s constant, k is Boltzmann’s constant, NA Avogadro’s number, n is the number of atoms per formula unit, M is the molecular mass per formula unit, q(=M/V) is the density, and vm is obtained from [52],

vm ¼

  1=3 1 2 1 þ 3 3 3 vt vl

ð23Þ

where vl and vt are the longitudinal and transverse elastic wave velocities, respectively, which are obtained from Navier’s equations [53]:

vl

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3B þ 4G ¼ 3q

vt ¼

q

The calculated density, longitudinal, transverse and average elastic wave velocities, and then Debye temperature for ThC are given in Table 10 for GGA, LDA and LDA + U functionals. No other theoretical or experimental data are exist for comparison with the present values. It is observed that the behaviors of the longitudinal, transverse and average wave velocities are nearly same. However, the Debye temperature increases up to 40 GPa, and after starting to decrease for all functionals. 3.5. Thermodynamic properties

ð24Þ

and

sffiffiffiffi G

67

ð25Þ

The thermal properties are determined in the temperature range of 0–1000 K for ThC, where the quasi-harmonic Debye model remains fully valid. The pressure effect is studied over 0–50 GPa range. The relationships between volume and pressure at different temperatures are shown in Fig. 10a, and those of between volume and temperature at different pressures are shown in Fig. 10b for

Fig. 15. Relative Grüneisen parameter to GGA results versus pressure for (a) LDA and (b) LDA + U, versus temperature for (c) LDA and (d) LDA + U.

Fig. 16. The variation of Cv with temperature at different pressures for ThC.

68

S. Aydin et al. / Journal of Nuclear Materials 429 (2012) 55–69

ThC. It is seen that when the pressure increases from 0 to 50 GPa, the volume decreases, this trend is convenient to the general expect. However, the volume is almost fixed with increasing temperature at constant pressure, it is changing ratio is very small. It can be expressed that pressure effects are more dominant on volume than temperature effects. In cases of LDA and LDA + U cases, results are higher than GGA results as pressure increases at constant temperature (see Fig. 11a and b), and difference is bigger at higher temperatures. Similarly, LDA and LDA + U results are higher than GGA results as temperature increases at constant pressure (see Fig. 11c and d), and difference is bigger at higher pressures. The variation of bulk modulus (B) with pressure at 0, 300, 700, 1000 K temperatures and the variation of bulk modulus with temperature at 0, 15, 30, 50 GPa pressures are shown in Fig. 12a and b for ThC, respectively. It can be seen that bulk modulus increases with pressure at a given temperature and decreases with temperature at a given pressure. And, the bulk modulus rapidly increases

– almost linear – with pressure, and effect of the temperature on the isothermal bulk modulus is very small. Furthermore, LDA and LDA + U results approaches to GGA results up to 25 GPa at constant temperature (see Fig. 13a and b), and after, difference between the results increase. As temperature increases, difference is higher. At the same time, LDA and LDA + U results are higher than GGA results as temperature increases at constant pressure (see Fig. 13c and d), and difference is bigger at higher temperatures. The calculated material properties at different temperatures are very sensitive to the vibrational modes. In the quasi-harmonic Debye model, the Grüneisen parameter, c(T) is a key quantity. The variation of Grüneisen parameter with pressure at various temperatures and with temperature at various pressures is given in Fig. 14a and b, respectively. It can be seen from Fig. 14 that Grüneisen parameter decreases with increasing pressure, while it increases with the increasing temperature. Moreover, in cases of LDA and LDA + U cases, results are being higher than GGA results

Fig. 17. (a) The variation of thermal expansion coefficient (a) with pressure at different temperatures and (b) with temperature at different pressures for ThC.

Fig. 18. Relative thermal expansion to GGA results versus pressure for (a) LDA and (b) LDA + U, versus temperature for (c) LDA and (d) LDA + U.

S. Aydin et al. / Journal of Nuclear Materials 429 (2012) 55–69

as pressure increases (after 25 GPa) at constant temperature (see Fig. 15a and b), and difference is bigger at higher pressures. Similarly, LDA and LDA + U results are higher than GGA results as temperature increases at constant pressure (see Fig. 15c and d), and difference decreases as temperature increases. The temperature-dependent behavior of the constant-volume heat capacity, Cv, at different pressures is shown in Fig. 16 for ThC (only GGA results are reported, because GGA, LDA and LDA + U results are very similar). At high temperatures, Cv tends to the Dulong–Petit limit [54]. At sufficiently low temperatures, Cv is proportional to T3 [54]. However, at intermediate temperatures, the temperature dependence of Cv is governed by details of vibrations of atoms. It is seen from Fig. 16 that when T < 300 K, Cv increases very rapidly with temperature, which is due to the anharmonic approximation of the Debye model. When T > 400 K the anharmonic effect on Cv is suppressed, and Cv increases slowly with temperature and it almost approaches to a constant, called as Dulong–Petit limit at higher temperatures. The thermal expansion coefficient (a) can be obtained from the temperature derivative of lattice constant. The variations of a with temperature at different pressures, and the variations of a with pressure at various temperatures are presented in Fig. 17a and b for ThC, respectively. In a–P graphs (see Fig. b), it is seen that a exponentially decreases with increasing pressure. In a–T graphs (see Fig. b), it is noted that a rapidly increases with T at low temperatures, and achieves a saturation at high temperatures. Furthermore, in cases of LDA and LDA + U cases, results are lower than GGA results up to 25 GPa at constant temperature (see Fig. 18a and b), and difference is bigger at higher pressures. Similarly, LDA and LDA + U results are lower than GGA results as temperature increases at 0 and 15 GPa pressures (see Fig. 18c and d), and are higher than GGA results at 30 and 50 GPa pressures. 4. Summary and conclusion In summary, the first principles calculations have been performed to obtain the structural, elastic, electronic, bonding, and hardness properties of thorium carbide with GGA, LDA and LDA + U functionals. Calculated elastic constants satisfy the traditional mechanical stability conditions at 0–50 GPa pressures. From band structure and density of states, thorium carbide is metallic at 0 GPa, and metallicity of ThC is originated from f-orbitals, and after d-orbitals, and it increases up to 40 GPa, and after, it decreases for all functionals. Other hand, considering of the Mulliken population analysis, electron density distribution maps and density of states together, it is revealed that the bonding in ThC is a mixture of covalent, ionic and metallic, and ionic component is dominate. The mechanical data such as Zener anisotropy factor, Poisson’s ratio, Young modulus, shear modulus are determined as a function of pressure. As a special interest, calculated hardness of ThC is 5.6 GPa at 0 GPa, thus thorium carbide is not superhard material, only hard material, and hardness properly increases with pressure. It is observed considerable changes in structural properties, elastic constants, hardness and metallicity around 40 GPa, and these may be an indicator for a possible phase transition. Furthermore, thermodynamical properties of the ThC compound have been presented in details for the first time by using quasi-harmonic Debye model in cases of GGA, LDA and LDA + U, and it is seen that this model generates acceptable results. Acknowledgements This work is supported partly by the State of Planning Organization of Turkey under Grant No. 2001K120590.

69

References [1] C.P. Mallett, J. Phys. C: Solid State Phys. 15 (1982) 6361. [2] L. Gerward, J.S. Olsen, U. Benedict, J.-P. Itie, J.C. Spirlet, J. Appl. Cryst. 19 (1986) 308. [3] J.-M. Fournier, R. Troc, in: A.J. Freeman, G.H. Lander (Eds.), Handbook on the Physics and Chemistry of the Actinides, vol. 2. p. 29 [4] P. Söderlind, G. Kotliar, K. Haule, P.M. Oppeneer, D. Guillaumont, MRS Bull. 35 (2010) 883. [5] R. Ducher, R. Dubourg, M. Barrachin, A. Pasturel, Phys. Rev. B 83 (2011) 104107. [6] B.-T. Wang, H. Shi, W. Li, P. Zhang, Phys. Rev. B 81 (2010) 045119. [7] A.J. Freeman, G.H. Lander (Eds.), Handbook on the Physics and Chemistry of the Actinides, North-Holland, Amsterdam, 1985. [8] R. Benz, A. Naoumidis, Thourium compounds with nitrogen, Gmelin Handbook of Inorganic Chemistry, eighth ed., vol. C3, Springer, Berlin, 1987 (Thorium supplement). [9] H. Kleykamp, Thorium Carbides, Gmelin Handbook of Inorganic and Organometallic Chemistry, eighth ed., vol. C6, Springer, Berlin, 1992 (Thorium supplement). [10] A.H.M. Evensen, R. Catherall, P. Drumm, P. Van Duppen, O.C. Jonsson, E. Kugler, J. Lettry, O. Tengblad, V. Tikhonov, H.L. Roun, Nucl. Instrum. Meth. Phys. Res. B 126 (1997) 160. [11] L. Gerward, J.S. Olsen, U. Benedict, J.P. Itie, J.C. Spirlet, J. Appl. Cryst. 19 (1986) 308. [12] J.S. Olsen, S. Steenstrup, L. Gerward, U. Benedict, J. -P. Itie, Physica 139+140B (1986) 308. [13] I.R. Shein, K.I. Shein, G.P. Shveikin, A.I. Ivanovskii, Dokl. Phys. Chem. 407 (2006) 106. [14] I.R. Shein, K.I. Shein, A.L. Ivanovskii, J. Nucl. Mater. 353 (2006) 19. [15] I.S. Lim, G.E. Scuseria, Chem. Phys. Lett. 460 (2008) 137. [16] I.R. Shein, A.L. Ivanovskii, J. Struct. Chem. 49 (2008) 348. [17] I.R. Shein, K.I. Shein, A.L. Ivanovskii, Tech. Phys. Lett. 33 (2007) 128. [18] I.R. Shein, A.L. Ivanovskii, Solid State Sci. 12 (2010) 1580. [19] T. Das, S. Deb, A. Mookerjee, Physica B 367 (2005) 6. [20] P. Hohenberg, W. Kohn, Phys. Rev. B 136 (1964) 864–871. [21] S.J. Clark, M.D. Segall, C.J. Pickard, P.J. Hasnip, M.J. Probert, K. Refson, M.C. Payne, Z. Kristallogr. 220 (2005) 567. [22] B.G. Pfrommer, M. Cote, S.G. Louie, M.L. Cohen, J. Comput. Phys. 131 (1997) 133–140. [23] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. [24] J.P. Perdew, A. Zunger, Phys. Rev. B 23 (1981) 5048–5079. [25] V. Anisimov, J. Zaanen, O.K. Andersen, Phys. Rev. B 44 (1991) 943. [26] M. Cococcioni, S. de Gironcoli, Phys. Rev. B 71 (2005) 035105. [27] D. Vanderbilt, Phys. Rev. B 41 (1990) 7892–7895. [28] M.A. Blanco, E. Francisco, V. Luano, Comput. Phys. Commun. 158 (2004) 57. [29] F. Peng, H.Z. Fu, X.L. Cheng, Physica B 400 (2007) 83. [30] F. Peng, H.Z. Fu, X.D. Yang, Solid State Commun. 145 (2008) 91. [31] F. Peng, H.Z. Fu, X.D. Yang, Physica B 403 (2008) 2851. [32] H. Fu, D. Li, F. Peng, T. Gao, X. Cheng, J. Alloy. Compd. 473 (2009) 255–261. [33] M.A. Blanco, A.M. Pendàs, E. Francisco, J.M. Recio, R. Franko, J. Mol. Struct. (Theochem) 368 (1996) 245. [34] M. Flórez, J.M. Recio, E. Francisco, M.A. Blanco, A.M. Pendàs, Phys. Rev. B 66 (2002) 144112. [35] E. Francisco, J.M. Recio, M.A. Blanco, A.M. Pendàs, A. Costales, J. Phys. Chem. 102 (1998) 1595. [36] U. Benedict, J. Less-Common Met. 128 (1987) 7. [37] M.D. Segall, R. Shah, C.J. Pickard, M.C. Payne, Phys. Rev. B 54 (1996) 16317. [38] L. Guan, B. Liu, L. Jin, J. Guo, Q. Zhao, Y. Wang, G. Fu, Solid State Commun. 150 (2010) 2011. [39] K. Haddadi, A. Bouhemadou, L. Louail, J. Alloy. Compd. 504 (2010) 296–302. [40] Y. Li, Y. Gao, B. Xiao, T. Min, Z. Fan, S. Ma, L. Xu, J. Alloy. Compd. 502 (2010) 28. [41] A. Simunek, J. Vackar, Phys. Rev. Lett. 96 (2006) 085501. [42] A. Simunek, Phys. Rev. B 75 (2007) 172108. [43] J. Sun, X.-F. Zhou, Y.-X. Fan, J. chen, H.-T. Wang, Phys. Rev. B 73 (2006) 045108. [44] W.B. Pearson, The Crystal Chemistry and Physics of Metals and Alloys, Wiley, NewYork, 1972. p. 151 (Table 4-4). [45] Z.J. Wu, E.J. Zhao, H.P. Xiang, X.F. Hao, X.J. Liu, J. Meng, Phys. Rev. B 76 (2007) 054115. [46] K. Haddadi, A. Bouhemadou, L. Louail, S. Bin-Omran, Solid State Commun. 150 (2010) 1995–2000. [47] A. Hao, X. Yang, X. Wang, R. Yu, X. Liu, W. Xin, R. Liu, Comp. Mater. Sci. 48 (2010) 59–64. [48] A.F. Young, C. Sanloup, E. Gregoryanz, S. Scandolo, R.E. Hemley, H.K. Mao, Phys. Rev. Lett. 96 (2006) 155501. [49] Y. Li, Y. Gao, B. Xiao, T. Min, Z. Fan, S. Ma, L. Xu, J. Alloy. Compd. 502 (2010) 28– 37. [50] J.P. Watt, L. Peselnick, J. Appl. Phys. 51 (1980) 1525. [51] J.R. Christman, Fundamentals of Solid State Physics, Wiley, New York, 1988. [52] O.L. Anderson, J. Phys. Chem. Solids 24 (1963) 909. [53] E. Screiber, O.L. Anderson, N. Soga, Elastic Constants and their Measurements, McGraw-Hill, New York, 1973. [54] A.T. Petit, P.I. Dulong, Ann. Chim. Phys. 10 (1819) 395.