Pressure and temperature dependent ab-initio quasi-harmonic thermoelastic properties of tungsten

We present the ab-initio temperature and pressure dependent thermoelastic properties of body-centered cubic tungsten. The temperature dependent quasi-harmonic elastic constants (ECs) are computed at several reference volumes including both the phonon and the electronic excitations contribution to the free energy and interpolated at different temperatures and pressures. Good agreement with the experimental ECs on a single crystal at ambient pressure is found. The pressure and temperature dependence of the shear sound velocity measured on polycrystalline tungsten by Qi et al is also in agreement with theory. Some discrepancies are found instead for the compressional velocity at high temperature and this is attributed to the temperature derivative of the bulk modulus, higher in theory than in experiment. These conclusions are reached both by PBE and by PBEsol functionals. The two give elastic properties with a similar pressure and temperature dependence although the latter is closer to experiment at 0 K.


Introduction
The study of elastic constants (ECs) of materials has a long history in mechanics, thermodynamics, and acoustics [1].EC describe the stresses induced in a solid by a strain as well as its mechanical stability.Moreover, the velocity of the sound waves can be derived from the EC and the solid density.Therefore ECs are key physical quantities for the design of any engineering application of materials.
The EC and the bulk modulus of tungsten are known from ultrasonic experiments at room pressure up to 2073 K [27][28][29] and have been calculated, at zero temperature, for several pressures [26].However information on high pressure and high temperature thermoelastic properties of tungsten is still incomplete.The pressure derivatives of the ECs at room temperature are known, [30] but only recently density as well as compressional and shear sound velocities have been measured by Qi et al [21] in polycrystalline tungsten up to 1073 K and 105 kbar using ultrasonic interferometry and x-ray diffraction.These measurements give experimental values of the adiabatic bulk and shear moduli at high pressures and temperatures.So far these data have not been compared with ab-initio calculations, but such comparison is timely both to test the ab-initio methods at high pressure and to further support the experimental measurement.Unfortunately, temperature dependent ECs (TDECs) within the quasi-harmonic approximation are numerically heavy to compute and the examples in the literature are focused mainly on materials of geophysical interest [31,32].To our knowledge, no quasi-harmonic calculation of the TDECs of tungsten exists and theoretical results are presently limited to molecular dynamics based on embedded atom method potentials at room pressure [33].
Recently, Malica and Dal Corso [34,35] presented a method to calculate ab-initio the quasi-harmonic TDECs in an automatic way, a method that is now implemented in thermo_pw [36].Applications, so far, to metals (Al, Cu, Ag, Au, Pd, Pt) and semiconductors (Si and BAs) are limited to room pressure, where the TDEC are computed only for a few geometries close to equilibrium.In this paper we report an ab-initio investigation of the thermoelastic properties of tungsten extending this computational scheme to high pressures.We compute the temperature dependent isothermal ECs from the second strain derivatives of the Helmholtz free energy, including both the quasi-harmonic (QHA) vibrational term and the contribution of electronic thermal excitations.These calculations, being carried out on several reference geometries at high and low pressures, can be used to interpolate the isothermal ECs at an arbitrary volume and therefore to study their temperature and pressure dependence.Adiabatic TDECs are calculated by thermodynamic relationships for comparison with experiments.From these ECs we derive the elastic parameters of polycrystalline tungsten, and compare with experiment.Finally, from the density and the polycrystalline elastic parameters we determine the compressional and shear sound velocities that are compared with the results of Qi et al [21].We confirm some of the experimental results, but we also find some discrepancies that might require further theoretical and/or experimental investigations.

Theory and computational parameters
In this paper we use thermo_pw [36] to calculate the thermodynamic properties.The QHA, as implemented in thermo_pw, has been discussed in previous publications [34,35,[37][38][39][40][41].Here we summarize the main formula and illustrate the finite pressure approach.Within the QHA, the Helmholtz free energy F(V, T) of a cubic solid is a function of temperature T and (unit cell) volume V.It can be written as the sum of three contributions: where U(V) is the static energy, F ph (V, T) is the vibrational free energy, and F el (V, T) is the electronic thermal excitations contribution.U(V) is computed via density functional theory (DFT), [42] F ph (V, T) is written in terms of the phonon frequencies ω η (q, V): Here h is the reduced Planck's constant and N is the number of cells of the solid (equal to the number of q points used in the sum).β = 1 kBT , where k B is the Boltzmann constant, q are the phonon wavevectors and η indicates the different modes.
The contribution of the electronic thermal excitations to the free energy is calculated within the rigid bands approximation as F el = U el − TS el .U el is electronic excitation contribution to the energy given by where E F is the Fermi energy, N(E) is the electronic density of states, f(E, T, µ) are the Fermi-Dirac occupations and µ is the chemical potential, and S el is the electronic entropy given by: These free energies are computed for a set of volumes V i , i = 1, N V .U(V) is interpolated by a fourth order Birch-Murnaghan equation while the vibrational and electronic free-energies are interpolated by a fourth-degree polynomial in V.In the supplementary material, we summarize the calculation of the volume thermal expansion β V , of the isobaric heat capacity C P and of the isoentropic bulk modulus B S .The isothermal ECs are calculated from the second strain derivatives of the free energy: correcting for finite pressure effects to obtain the stress-strain ECs [43]: The second derivatives of the free energy are calculated as described in [34] taking a subset of the volumes V i as equilibrium configurations.The ECs at any other volume at temperature T and pressure p are obtained by interpolation by a fourth-degree polynomial.Adiabatic ECs are calculated from the isothermal ones as: where b ij are the thermal stresses: and C V is the isochoric heat capacity (see the supplementary material).For a cubic system, the linear thermal expansion tensor is diagonal and from the volume thermal expansion β V we get: α kl = δ kl β V /3.With the knowledge of the adiabatic ECs, using the Voigt-Reuss-Hill approximation, we compute the polycrystalline average of the bulk modulus B S , of the shear modulus G S , of the Young's modulus E S , and of the Poisson's ratio ν S (see the thermo_pw manual for the expressions used).Finally, the compressional and the shear sound velocities are given by: where ρ is the density.Note that these equations also hold with pressure as long as B S and G S are computed from the stressstrain ECs [44].
The calculations of the TDECs presented in this work were done by using DFT implemented in Quantum ESPRESSO [45,46] with the PBE [47] exchange and correlation functional.In addition, we present also some results using the LDA [48] and PBEsol, [49] a functional that, modifying the PBE exchange, gives a better description of solids [50] at the expense of the accuracy in molecules.We employ the projector augmented wave (PAW) method [51] and a planewave basis with pseudopotentials generated by us starting from those available in pslibrary [52,53].We called these pseudopotentials W.pz-spn-kjpaw_psl.1.0.1.UPF, W.pbesol-spn-kjpaw_psl.1.0.1.UPF, and W.pbe-spn-kjpaw_psl.1.0.1.UPF for LDA, PBEsol, and PBE, respectively and added them to pslibrary.They can be obtained from the web page given in [53].They have 5s, 5p, 5d, and 6s valence states, while the 4f states are frozen in the core and accounted for by the nonlinear core correction [54].For the wave functions cut-offs, we use 70 Ry, 90 Ry, 90 Ry while for the charge density we use 280 Ry, 360 Ry, 360 Ry, for LDA, PBE, and PBEsol, respectively.The Fermi surface has been dealt with by a smearing approach [55] with a smearing parameter σ = 0.02 Ry.With this smearing, the Brillouin zone integrals converge with a 40 × 40 × 40 k-point mesh.Density functional perturbation theory (DFPT) [56,57] is used to calculate the dynamical matrices on a 8 × 8 × 8 q-point grid.These dynamical matrices have been Fourier interpolated on a 200 × 200 × 200 q-point mesh to evaluate the free-energy and thermodynamic quantities.For quasi-harmonic calculations, the free energy was calculated in N V = 15 geometries with lattice constants ranging from a 0 − 1.2 a.u. to a 0 + 0.2 a.u. in steps of ∆a = 0.1 a.u., where a 0 is the equilibrium 0 K lattice constant (see table 1).The ECs C 11 , C 12 , and C 44 are calculated by using six strained configurations for each type of strain (see [34]) with δε = 0.005.Therefore, each reference configuration requires the phonon dispersions and the electronic density of states for 18 geometries, 6 with a body-center cubic lattice, 6 with centered tetragonal lattice, and 6 with a rhombohedral lattice.Phonons are calculated on a 8 × 8 × 8 q-point grid, sampling the Brillouin zone by a 45 × 45 × 45 k-point mesh.Among the N V = 15 geometries used for anharmonic calculations, we computed the QHA ECs for six reference geometries: 1, 5, 8, 12, 13, and 14.The phonon frequencies of geometries 1, 5, 8, and 12 have been computed on Marconi100 at CINECA with an accelerated GPU version of thermo_pw optimized for calculations with a dense grid of k-points [63].

Results and discussion
We start in table 1 by comparing the LDA, PBEsol, and PBE equilibrium lattice constants, bulk moduli, and pressure derivatives of the bulk moduli obtained from the interpolation of U(V) with selected previous calculations and experiments.Using a exp = 5.972 a.u. as the 0 K experimental lattice constants obtained subtracting 0.009 a.u., the thermal expansion contribution, to the reported 300 K value of [20], LDA underestimates the lattice constant by 0.7% while PBE overestimates it by 0.8%.PBEsol is closer to the experiment underestimating it by 0.1%.This agrees with the all-electron LAPW calculation of [61].Using B = 3142 kbar as the 0 K experimental bulk modulus (the value given in [27]), the LDA and PBEsol overestimations are about 8% and 4%, while PBE underestimation is about 2%.The three functionals give similar values for the pressure derivative of the bulk modulus which are in reasonable agreement with the experiment.In the supplementary material we show, as a reference, the phonon dispersions, the thermal EOS, the volume thermal expansion, the isobaric heat capacity, the bulk modulus, and the average Grüneisen parameter calculated with LDA, PBEsol, and PBE and compare them with the available experiments and previous calculations.We refer to these data when needed in the following.
The PBE and PBEsol thermal EOS in the range of pressures and for the temperatures measured in the experiment of Qi et al [21] (298 K, 473 K, 673 K, 873 K, and 1073 K) are shown in the upper part of figure 1.Since the PBE functional overestimates the volume, its predicted densities are all below experiment.The PBEsol densities are much closer to the experimental data although slightly higher.In order to facilitate the comparison of the temperature and pressure dependence of the theoretical densities and experiment we superimpose the two sets of curves in the lower part of figure 1 making a upper shift of ∆ρ = 0.475 g cm −3 of all the PBE curves and of ∆ρ = −0.06g cm −3 of all the PBEsol curves so that the two curves at 298 K overlap with experiments at the point at 53 kbar.This comparison shows that the slope of the lines and the distance between lines are very similar to experiment for both functionals as expected from the bulk moduli and thermal expansion.The pressure slope of the PBEsol density is slightly smaller than the PBE one (in agreement with the different bulk moduli), but the difference is not significant in the comparison with experiment in this range of pressures.
The LDA, PBEsol, and PBE ECs calculated at the 0 K equilibrium volume are reported in table 2 and compared with experiments and selected earlier calculations.The experimental values of Lowrie and Gonas [28] refer to 273.15 K, but they can be extrapolated at 0 K using our PBE differences between 0 K and 273 K.We obtain a good agreement with the 0 K values of [27] and take these extrapolated values as experimental reference.The LDA errors are 361 kbar (7%), 154 kbar (8%), and −118 kbar (−7%) for C 11 , C 12 , and C 44 , respectively.PBEsol errors are 155 kbar (3%), 84 kbar (4%), and −119 kbar (−7%) while PBE errors are −181 kbar (−3%), −50 kbar (−2%), and −194 (−12%).In table 2 we also report the bulk, shear, and Young's moduli of polycrystalline tungsten derived from the ECs.The bulk moduli derived from the ECs are about 1% smaller than those reported in table 1 derived from the Birch-Murnaghan interpolation.This shows the accuracy of our ECs calculations.The difference is within the numerical uncertainty of the Birch-Murnaghan fitting.
Pressure dependent ECs calculated with the three functionals are reported in figure 2 and compared with the previous PBE calculation [26].On the scale of this figure the three functionals are almost equivalent and in good agreement with earlier results.The numerical effort to calculate the temperature and pressure dependent ECs is large and we could afford only two functionals, so we chose PBEsol and PBE which have smaller errors than LDA.Presently, we do not include phonon-phonon anharmonic effects in the free energy, so our results are expected to be reliable only up to 1500 K, a temperature for which, as shown in the supplementary material, the quasi-harmonic approximation is sufficient to describe accurately the other thermodynamic properties of this system.In figure 3 we report the temperature dependent isothermal and adiabatic ECs as a function of temperature at room pressure.As is well known from theory, in cubic systems adiabatic and isothermal ECs differ only for C 11 and C 12 while they coincide for symmetry reasons for C 44 .Experimentally ECs are derived from sound velocities and are adiabatic.We compare our data with Lowrie and Gonas [28] which provided analytic fits of their experiments on single crystal tungsten.The points shown in the figure have been obtained from these fits.The temperature dependence of the ECs is well reproduced, theory and experiment differ only for a rigid shift.From 297 K to 2073 K the experimental values decrease by 993 kbar (19%), −44 kbar (−2%), and 250 kbar (15%) for C 11 , C 12 , and C 44 while the theoretical ones decrease by 953 kbar (18%), −12 kbar (−0.6%) and 289 kbar (20%) with PBEsol and by 874 kbar (17%), 17 kbar (0.8%), and 263 kbar (16%) with PBE.The behavior of C 12 is quite unusual.In experiment C 12 is almost constant with a slight enhancement with temperature.We find an almost constant curve, slightly increasing with PBEsol and slightly decreasing with PBE.
The bulk modulus, shear modulus, and Young's modulus of polycrystalline tungsten calculated by the Voigt-Reuss-Hill  approximation are shown against temperature in figure 4 and compared with the analytic fits of the experimental data of [28].Again except for the absolute positions, the temperature dependence of B, E, and G is reproduced reasonably well.
From 297 K to 2073 K, these quantities decrease by 10%, 21%, and 23% both in experiment and in our PBE calculation.With PBEsol these figures become 10%, 22%, and 24%, but the absolute position is closer to experiment.The same is found for the Poisson's ratio ν, shown in the inset where the increase is 10% in experiment and 9% (PBE) and 10% (PBEsol) in theory.At 300 K and 0 kbar, the pressure derivative of the adiabatic bulk modulus is dBS dp  PBEsol (green lines) and PBE (blue lines) temperature dependent isothermal (dashed) and adiabatic (solid) elastic constants calculated within the quasi-harmonic approximation at zero pressure.The gold and red circles indicate the adiabatic experimental data of [29] and of [28], respectively.For C 44 isothermal and adiabatic ECs coincide.
In figure 5 we show the pressure dependent longitudinal modulus that derives from our ECs.This modulus is computed as L = B S + 4 3 G S .A similar calculation is done also with the experimental data of [21].Also in these curves it is easier to compare the two functionals and experiment by doing a shift that removes the differences due to the T = 0 K ECs.The figure with unshifted curves is presented in the supplementary material.We make a shift of 437 kbar for PBE and of 128 kbar for PBEsol so that the 298 K curves and the experimental point at 53 kbar and 298 K coincide.The theoretical data reproduce accurately the pressure dependence of the experimental data and both functionals give curves with the same slope, however our calculation predicts a larger variation of L with temperature than experiment.This difference is mainly due to the bulk modulus as can be deduced by comparing the shear modulus (see figure 6).After a shift of 204 kbar for PBE and of 122 kbar for PBEsol (again to make the curve at 298 K to pass through the point at 53 kbar), the temperature and pressure dependence PBEsol (green lines) and PBE (blue lines) temperature dependent isothermal (dashed) and adiabatic (solid) polycrystalline averages of the macroscopic elastic properties (bulk modulus B, Young's modulus E, and shear modulus (G) calculated within the quasi-harmonic approximation at zero pressure.The red circles are the adiabatic data of [28].In the inset we compare the isothermal and adiabatic Poisson's ratio with experiment. of the shear modulus follows with experiment with the PBEsol shear modulus that decreases slightly faster than the PBE one with temperature.
Finally, in figures 7 and 8 we present a comparison between the sound velocity data and our results (equations ( 9) and ( 10)).We show our theoretical PBE results and the same curves shifted by 146 m s −1 (figure 7) and 149 m s −1 (figure 8) so that the 298 K curves are above the experimental data at 53 kbar.In the same figures we show also the PBEsol data shifted by  65 m s −1 (V P ) and 113 m s −1 (V G ).As expected from previous analysis, we find that the pressure dependence is well reproduced for both the compressional and the shear sound velocities.The temperature dependence of the shear sound velocity is also very well accounted for, while the distances between compressional velocities at different temperatures agree with our curves at low temperatures, but the experimental points at 873 K and those at 1073 K have higher values than our theoretical curves.These findings are due to the values of dBS dT and dGS dT , the latter in quite good agreement with experiment, while the former only in fair agreement.

Conclusions
We presented the quasi-harmonic temperature and pressure dependent thermoelastic properties of tungsten calculated by thermo_pw using PAW pseudopotentials.The temperature dependent quasi-harmonic adiabatic ECs as well as the derived polycrystalline bulk, shear, and Young's moduli and Poisson's ratio have been compared with experiment at ambient pressure.The compressional and shear sound velocities, as well as the longitudinal and shear moduli, have been calculated as a function of pressure at the temperatures measured in the experiment of Qi et al [21].We find that the quasi-harmonic theory reproduces well the temperature dependent ECs of tungsten at room pressure.The PBE and the PBEsol functionals give similar temperature dependences with PBEsol closer to experiment than PBE.Good agreement is also found for the temperature dependence of the polycrystalline elastic moduli, again with PBEsol closer to experiment than PBE.
The pressure and temperature dependence of the shear modulus of polycrystalline tungsten in the range of pressures measured in the experiment of Qi et al [21] is reproduced very well by both PBE and PBEsol that are practically indistinguishable when a shift is used to facilitate the comparison with experiment.The pressure and temperature variations of the longitudinal module are instead described less accurately and also after a shift PBEsol and PBE predict a larger decrease with temperature than experiment.The PBEsol temperature derivative of the longitudinal module is slightly larger than the PBE one.These discrepancies are attributed to the temperature derivative of the adiabatic bulk modulus larger in the calculation than in experiment.As a consequence, the pressure and temperature dependences of the shear sound velocity are well reproduced by theory, while the compressional sound velocity as a function of pressure is well reproduced at low temperatures, but larger and larger differences are found at high temperatures.These conclusions hold for both PBE or PBEsol.Different functionals change the position of the elastic moduli or sound velocity curves, but not their slope nor their temperature variation.The reasons for the discrepancies between theory and experiment are presently unclear and might require an improvement of the theory or a revision of the experiment.
In this paper we could compare the elastic moduli and the sound velocities of polycrystalline tungsten at high pressure and high temperature with experiment and we demonstrated that these calculations are now routinely feasible in modern supercomputers.These parameters which are quite critical for high pressure and high temperature applications of materials are still missing from the literature for many metals.We hope that our paper will stimulate further measurements and theoretical calculations in this direction.

Figure 1 .
Figure 1.Upper: pressure dependent PBEsol (dashed lines) and PBE (solid lines) density calculated for several temperatures (red 298 K, green 473 K, blue 673 K, yellow 873 K, and pink 1073 K) compared with the experimental data of Qi et al [21].Lower: the PBE density is shown after a shift of ∆ρ = 0.475 g cm −3 (indicated by the red arrow on the line at 298 K) to facilitate the comparison with experiment.The dashed lines are the PBEsol curves after a shift of −0.06 g cm −3 also shown by a red arrow.

Figure 2 .
Figure 2.Pressure dependent ECs calculated at 0 K obtained from the second derivatives of the energy with respect to strain for the N V geometries used for the quasi-harmonic approximation.ECs are calculated within LDA (red line), PBEsol (green line), and PBE (blue line).The red circles are the PBE data of[26].

Figure 3 .
Figure 3.PBEsol (green lines) and PBE (blue lines) temperature dependent isothermal (dashed) and adiabatic (solid) elastic constants calculated within the quasi-harmonic approximation at zero pressure.The gold and red circles indicate the adiabatic experimental data of[29] and of[28], respectively.For C 44 isothermal and adiabatic ECs coincide.

Figure 4 .
Figure 4.PBEsol (green lines) and PBE (blue lines) temperature dependent isothermal (dashed) and adiabatic (solid) polycrystalline averages of the macroscopic elastic properties (bulk modulus B, Young's modulus E, and shear modulus (G) calculated within the quasi-harmonic approximation at zero pressure.The red circles are the adiabatic data of[28].In the inset we compare the isothermal and adiabatic Poisson's ratio with experiment.

Figure 5 .
Figure 5. Temperature dependent longitudinal modulus calculated within the quasi-harmonic approximation as a function of pressure for several temperatures (red 298 K, green 473 K, blue 673 K, yellow 873 K, and pink 1073 (K).The circles (with the same color code) are the data of [21], measured at the same temperatures and calculated as L = B S + 4/3G S from the B S and G S in their table III.To facilitate the comparison with experiment, theoretical lines translated by 437 kbar, as indicated by the red arrows, are also shown.PBEsol results are shown with dashed lines after a shift of 128 kbar.

Figure 6 .
Figure 6.Temperature dependent PBE shear modulus calculated within the quasi-harmonic approximation as a function of pressure for several temperatures (red 298 K, green 473 K, blue 673 K, yellow 873 K, and pink 1073 K).The circles (with the same color code) are the data of [21], measured at the same temperatures (from their table III).To facilitate the comparison with experiment, theoretical lines translated by 204 kbar, as indicated by the red arrows, are also shown.The PBEsol results are shown with dashed lines after a shift of 122 kbar.

Figure 7 .
Figure 7. Temperature dependent compressional sound velocity (V P ) calculated within the quasi-harmonic approximation as a function of pressure for several temperatures (red 298 K, green 473 K, blue 673 K, yellow 873 K, and pink 1073 K).The circles (with the same color code) are the data of [21], measured at the same temperatures.To facilitate the comparison with experiment, theoretical lines translated by 146 m s −1 , as indicated by the red arrows, are also shown.PBEsol results are shown by dashed lines after a shift of 65 m s −1 .

Figure 8 .
Figure 8. Temperature dependent shear sound velocity (V G ) calculated within the quasi-harmonic approximation as a function of pressure for several temperatures (red 298 K, green 473 K, blue 673 K, yellow 873 K, and pink 1073 K).The circles (with the same color code) are the data of [21], measured at the same temperatures.To facilitate the comparison with experiment, theoretical lines translated by 149 m s −1 , as indicated by the red arrows, are also shown.PBEsol results are shown by dashed lines after a shift of 113 m s −1 .

Table 1 .
The equilibrium lattice constants (a 0 ), the bulk moduli (B T ) and the pressure derivatives of the bulk moduli (B ′ T ) of tungsten calculated in this work compared with selected previous calculations and with experiment.
a These data are used to calculate the equations of state reported in figure3of the supplementary material.b Ultrasonic experiment.c Shock wave experiment.

Table 2 .
The 0 K elastic constants calculated with the different functionals compared with experiment and selected previous calculations.B, E, G, and ν are the bulk modulus, the Young's modulus, the shear modulus, and the Poisson's ratio, respectively.