Magneto-elastic effects and thermodynamic properties of ferromagnetic hcp Co

Physica B 441 (2014) 72–79

Contents lists available at ScienceDirect

Physica B journal homepage: www.elsevier.com/locate/physb

Magneto-elastic effects and thermodynamic properties of ferromagnetic hcp Co Fang-Guang Kuang, Xiao-Yu Kuang n, Shu-Ying Kang, Ai-Jie Mao Institute of Atomic and Molecular Physics, Sichuan University, Chengdu 610065, China

art ic l e i nf o

a b s t r a c t

Article history: Received 18 November 2013 Received in revised form 21 January 2014 Accepted 10 February 2014 Available online 18 February 2014

Using first principles projector augmented wave (PAW) potential method, the magneto-elastic effects and thermodynamic properties of ferromagnetic hcp Cobalt at high pressure and temperature are investigated. The calculated elastic constants from PBE þU method demonstrate a noticeable improvement with regard to experimental data. Various physical quantities under high pressure also present significant improvements, such as the bulk modulus, shear modulus, Young's modulus, Debye temperature, various sound velocities and the normalized acoustic velocities in the meridian plane. That is due to the fact that Cobalt system possesses large correlation effects. Meanwhile, the phonon dispersion curves are in excellent agreement with experimental data. It is not observed any anomaly or instability under compression. However, according to the E2g-phonon frequencies, the obtained pressure variation of C44 elastic modulus also suggests that the system has miraculous magneto-elastic effects. Moreover, the pressure and temperature dependence of thermodynamic properties are derived within the quasiharmonic approximation for the first time. The obtained Grüneisen ratio, Anderon–Grüneisen parameter and the volume dependence of Grüneisen ratio display manifestly temperature and pressure dependences. & 2014 Elsevier B.V. All rights reserved.

Keywords: First-principles Magneto-elastic effects Phonon dispersion quasi-Harmonic approximation Cobalt

1. Introduction The magnetic 3d transition metal under compression and expansion is an outstanding model system to understand the composition and structure of the Earth's core better. For instance, the study of hexagonal close-packed (hcp) Fe under high pressure and temperature is of great geophysical and high-pressure physical interest. However, the hcp Fe has been considered to be nonmagnetic and no single crystal of hcp Fe has been grown to date. Since studying of elasticity under compression on ferromagnetic hcp Fe are inundated with challenges, more and more attention are focused on the metal adjacent to Fe in the Periodic Table, i.e. hcp Cobalt. The study of Co provides a number of advantages over Fe to study the hcp structure of 3d transition metals. Because the pressure–temperature (P–T) phase diagram of Co shows that the ferromagnetic hcp Co has no phase transition until temperature goes higher than 695 K at ambient pressure and can be stable up to 100 GPa at room temperature. Recent experiment has proved that the magnetic moment of hcp Co can be even stable up to higher temperature under certain chemical conditions [1].

n

Corresponding author. Tel./fax: þ 86 2885405515. E-mail address: [email protected] (X.-Y. Kuang).

http://dx.doi.org/10.1016/j.physb.2014.02.017 0921-4526 & 2014 Elsevier B.V. All rights reserved.

Elastic constants directly describe the response of the crystal to external stresses and are essential for many practical applications related to the mechanical properties of materials, such as load deflection, thermo-elastic stress, fracture toughness, etc. Since modern high-pressure advances in diamond-anvil-cell techniques [2], there were lots of experimental techniques on the high-pressure elasticity of single hcp Co, such as inelastic x-ray scattering (IXS) [3], impulsive stimulated light scattering (ISLS) [4], and radial x-ray diffraction (XRD) [5] and so on. Nevertheless, these studies did not reach the critical pressure region above 60 GPa. Using ISLS and Raman scattering methods to investigate the elastic and vibrational properties of hcp Co up to 120 GPa, Goncharov et al. discovered that the Raman active transverse optic (TO) phonon mode showed a change in slope near 60 GPa [6]. The anomalous behavior was confirmed by Antonangeli et al. through determining the longitudinal acoustic phonon dispersion of polycrystalline Co by IXS [7]. It is observed that the slope of the linear variation of sound velocities immediately as the pressure above 75 GPa. But this change occurred well below the phase transition pressure, it was attributed to magneto-elastic effects, which were supported by the onset of a loss of magnetic moment in hcp Co predicted by the density functional theory (DFT) computations [8–10]. Generally, first-principles calculations have been improved in reliability and accuracy and have been widely used to predict the high pressure–temperature behavior, deepening the fundamental

F.-G. Kuang et al. / Physica B 441 (2014) 72–79

understanding of the experiments. However, direct DFT calculations are often difficult to represent the magnetic properties of Co and thus cause a considerable change in aggregate elastic properties. In this case, the high pressure elastic properties agree qualitatively with experimental data, while they have greatly deviations from experimental values in the quantitative. In this paper, using different computational methods based on first-principles, the problem of how to obtain accurately elastic properties of ferromagnetic hcp Co is solved. Despite that the thermodynamic properties of ferromagnetic hcp Co were studied experimentally in considerable details, there was less work on thermodynamic properties of hcp Co using increasingly pervasive DFT. Most of theoretical thermodynamic studies on ferromagnetic hcp Co were performed by semiempirical theory. They were used to choose a model that was particularly useful for the treatment of the thermodynamic properties of the 3d ferromagnetic elements. The model was always based upon the phenomenological description of the magnetic and the nonmagnetic contributions to the Gibbs energy of the substance. Recently, the supercell method is established as a good method for calculating the vibrational and thermodynamic properties. Therefore, investigating the magneto-elastic effects from lattice vibrations is another aim of our work. Furthermore, under quasi-harmonic approximation (QHA), the study of thermodynamic properties at high temperature and high pressure is also the fundamental understanding of the experiments. This paper is organized as follows: in Section 2, we are outlining the method used for the first-principles calculations, along with practical details of the calculations. In Section 3, we apply PAW method based on DFT at high pressure and high temperature of ferromagnetic hcp Co by expanding on previous computational work. The calculated results of elastic, vibrational and thermodynamic properties of ferromagnetic hcp Co are reported and discussed. Finally, in Section 4, we present our summary and conclusions.

73

and the anisotropies of the wave polarized perpendicular to the basal plane (S1) and the polarized one in the basal plane (S2) by followed: ΔS1 ¼(C11 þ C33 2C13)/4C44, ΔS2 ¼ 2C44/(C11  C12). 2.1.2. Thermodynamic properties Under quasi-harmonic approximation (QHA), the Gibbs function G(V; p, T) takes the form GðV ; p; TÞ ¼ E0 ðV Þ þ F vib ðT; VÞ þ F elec ðT; VÞ þ pV:

ð3Þ

In Eq. (3), the E0(V) is the static contribution to the internal energy at volume V and is easily obtained from standard DFT calculations. The pV corresponds to the constant hydrostatic pressure condition, Felec(T;V) is the thermal free energy arising from electronic excitations that is also evaluated from firstprinciples calculations directly, and Fvib(T;V) is the vibrational free energy that comes from the phonon contribution. Within the QHA, Fvib(T;V) is given by    hωq;λ ðVÞ : ð4Þ F vib ðT; V Þ ¼ kB T∑In 2 sinh 4π kB T q;λ In Eq. (4), the sum is over all three phonon branches λ and over all wave vectors q in the first Brillouin zone. kB is the Boltzmann constant, h is the Planck constant and ωq,λ(V) is the frequency of the phonon with wave vector q and polarization λ, evaluated at constant volume V. Once the phonon spectrum is obtained, the temperature dependence of the entropy Svib ¼  (∂Fvib/∂T)V and the internal energy Uvib ¼ Fvib þTSvib at constant volume are easily computed. Thenceforth the specific heat capacity is Cvib V ¼(∂Uvib/ ∂T)V. The electronic specific heat at low temperature is obtained from C elec ¼ γ elec T; V

ð5Þ

where γelec is the electronic heat capacity coefficient. The total elec specific heat at constant volume as CV ¼ Cvib V þCV . The specific heat at constant pressure Cp is computed by using the relation

2. Theoretical and computational details

C p ¼ C V þ αV BT VT;

2.1. Theory

where the coefficient of volume thermal expansion coefficient αV and the isothermal bulk modulus BT are obtained from Vinet equation of state (EOS) [16] fitting.

2.1.1. Elastic properties Aggregate elastic properties, such as the bulk modulus B, the shear modulus G and the Young's modulus Y can be derived from the single-crystal elastic constants by Voigt–Reuss–Hill average [11,12]. Combining with the density ρ, the aggregate compressional sound velocity vL, shear sound velocity vT and bulk sound velocity vB are readily obtained from standard relations. The average sound velocity vm can be accurately estimated [13]. Once the vm is determined, the Debye temperature ΘD can be expressed [14]. The acoustic vibrational modes of a solid in the longwavelength limit are obtained from the Christoffel equation [15]

ρω2 ui  C ijkl kj kk kl ¼ 0;

ð1Þ

where ω is the vibrational angular frequency, ui is the displacement amplitudes and summation convention is used for repeated Cartesian indices, Cijkl is the elastic constant, and kj is the wave vector components of the vibrational wave. In case of a hexagonal crystal, the acoustic anisotropy can be described as

Δi ¼ Mi ½nk =Mi ½1 0 0

ð2Þ

where M ¼ ρω and hence ω ¼ Cn(k)k with Cn(k) the sound velocity for polarization n and direction k, and i denotes the three types of elastic wave (one longitudinal and two polarizations of the shear wave). By solving the Cristoffel equation for hcp Co, the anisotropy of the compressional wave (P) can be obtained from ΔP ¼C33/C11, 2

ð6Þ

2.2. Details of the calculations The static total energy and elastic properties were calculated using the projector augmented wave (PAW) potential method with three different exchange correlation potentials, namely, the local density approximation (LDA) in the scheme of Ceperley and Alder [17] as parameterized by Perdew and Zunger [18], generalized gradient approximation (GGA) of Perdew–Wang (PW91) [19] and Perdew–Burke–Ernzerhof (PBE) [20], as implemented in the Viena Ab Initio Simulation Package (VASP) [21–23]. All pseudopotetials based on the projector augmented wave method explicitly included the Co 3d84s1 electrons in the valence states. To properly describe the magnetic behavior of hcp Co phase, an accurate treatment of the electron correlation in the localized d orbital is crucial. Thus, several typical interactions were employed, such as spin polarization (SP), spin–orbit coupling (SOC) and Hubbard U correction. Hubbard U correction is shown to provide a quite accurate description of the electronic structure and properties of several transition metal and rare-earth systems [24–26]. In present work, the Dudarev implementation with on-site coulomb interaction U¼2.8 eV and on-site exchange interaction J ¼1.0 eV were used which was close to the value used for bulk Co [27]. A plane wave basis set with a cut-off of 500 eV was used and specific

74

F.-G. Kuang et al. / Physica B 441 (2014) 72–79

K-points were selected using a 13  13  9 grid for integration within the Brillouin zone. Tolerances for geometry optimization were set as the difference in total energy being within 10  5 eV per atom. The elastic constants were carried out by computing the change in energy under finite volume-conserving strains. For the hcp structure, it had five independent constants, namely, C11, C12, C13, C33, and C44. They were determined by applying a slight strain to an optimized unit cell through five different strain modes. The maximum strain amplitude was set from  0.015 to 0.015 in steps of 0.003, and the force on each atom was less than 0.005 eV per angstrom. Under high pressure, it was favorable to separate the effect of hydrostatic pressure. So the elastic constants were defined as relating stress and strain deviations relative to the hydrostatic state. The other elastic moduli (e.g. the shear modulus G, the Young's modulus Y and the bulk modulus B) of hcp cobalt were estimated according to the Voigt–Reuss–Hill approximation. For Eq. (3), the PBE exchange-correlation functional with spin polarization as implemented in the VASP code was employed to compute E0(V) and Felec(T;V). The Fvib(T;V) was obtained from firstprinciples phonon calculations by using Phonopy [28,29], which could support the VASP interface to calculate force constants directly through finite displacement method. In comparison, the 2  2  2 supercell with 16-atom was determined to calculate phonon dispersions. The total Gibbs free energy at temperature points were calculated with a step of 10 K from 0 to 1500 K at 12 volume points. At each temperature point, the equilibrium volume V0T and isothermal bulk modulus B0T were gotten by minimizing Gibbs free energy with respect to volume by fitting the integral form of the Vinet EOS at constant pressure.

3. Results and discussion 3.1. Elastic properties The elasticity of transition metals under compression is currently attracting much attention both in experiments and theory due to its significant contribution to the study of the internal structure of the Earth and other planets. The elastic moduli of the single crystal hcp Co under high pressures were studied by using IXRS [3], ISLS [4], XRD [5] and ultrasonic waves (US) [30]. Although these measurements have not reached above 60 GPa, the wealth of experimental values provide useful references for the theoretical calculations. The calculated structural properties and elastic moduli at zero pressure are listed in Table 1, together with the available experimental and theoretical data. The LDA potentials are the first choice due to its superiority for uniform electron gas system like bulk material of metal. Nevertheless the calculated elastic constants largely deviated from the experimental values (specifically, compared with Ref. [30], elastic constants are overestimated  35% in average). These disappointing results may be due to the greatly underestimated equilibrium volume. Based on the widely applicability in the first principles calculations, the GGA method both within PW91 and PBE potentials are employed. Compared with experiment [31], the V0 and c/a are better than LDA and the calculated elastic constants have more positive results. In particular, the off-diagonal constants by PBE method (i.e. C12 and C13) only showed  2% errors compared with experiments [3,30]. However, both the longitudinal (C11 and C33) and the shear (C44 and C66) elastic constants are still overestimated  20%. In order to reduce these deviations, the spin-order coupling is considered within PBE method, named PBE þSOC. Compared with preceding methods, the axial ratio of 1.622 is closer to the experimental result [31]. Confusingly, the elastic constants are still substantial

Table 1 The elastic constants (in GPa), atom equilibrium volume V0 (in Å3), axial ratio c/a and magnetic moment M (in μB) of ferromagnetic hcp Co at 0 GPa by various methods, together with experimental measurements and others calculations. Method

V0

c/a

M

C11

C12

C13

C33

C44

C66

Ref.

LDA PW91 PBE PBEþ SOC PW91þ U PBEþ U IXS ISLS US LAPW-LDA LAPW-GGA

10.02 10.86 10.88 10.88 10.82 10.86 – – – 9.32 10.37

1.611 1.615 1.616 1.622 1.617 1.619 – – – – –

1.48 1.57 1.58 1.59 1.70 1.71 – – – – –

472 368 369 378 320 304 293 303 307 515 440

215 172 169 151 151 166 143 161 165 245 210

160 120 105 123 114 116 90 107 102 175 140

510 411 414 368 351 359 339 – 358 575 485

118 99 99 94 82 83 78 72.9 75 160 125

128 97 100 114 84 69 75 71.2 71 135 115

[3] [4] [30] [9] [9]

deviations. The comparison between these values, the correction of spin-order coupling cannot essentially change the accuracy of the elastic constants. For the purpose of reducing the errors of elastic moduli, a special approach is employed, i.e. Hubbard U correction. Hubbard U correction is considered within PW91 and PBE functionals for a more comprehensive discussion, named PW91 þU and PBEþ U, respectively. It is observed that V0 and c/a have no significant change. The concern is the value of magnetic moment (  1.70 μB). It is a satisfactory agreement with DFT calculation using the linearized augmented plane wave method (LAPW) [9] and is also consistent with experiment [32]. In this case, the obtained elastic constants are quite close to experimental data, especially, results by PBEþU method. In summary, considering the Hubbard U correction can better describe the elasticity of ferromagnetic hcp Co. Why do the obtained elastic constants by PBE þ U show better than other methods? The mainly reason is that the introduction of U parameter gives rise to an electronic redistribution in the d-band structure, which leads to variations in the magnetic properties. The electronic density of states provides an alternative important way to understand this problem. Based on PBE and PBE þ U methods, the total, Co-4s and Co-3d density of states are shown in Fig. 1. Upon inclusion of the on-site Hubbard electron correlation correction to PBE parameter, the calculated value of the Fermi level (Ef) decreases by 0.03 eV and a shift in the energy of the up spin states toward more negatives values is obtained. Obviously, the inclusion of the U parameter leads to a variation in the 3d band structure, giving rise to a redistribution of the density of states in the subbands, leading to an increase in the spin down density that corroborates the enhancement observed in the magnetic properties. And thus cause a proper change in aggregate elastic properties. In order to investigate the pressure-dependence behaviors of the elastic properties, we concentrate our investigation on the pressure range below 100 GPa. The pressure dependences of elastic moduli are displayed in Fig. 2. Both results of PBE and PBEþU method are shown, along with the reported measurements by IXS and ISLS as well as other calculations. Both PBE and PBEþU method are well qualitative description of the relationship between elastic constants and pressure, while the PBE þU method is better quantitative description. So the following analysis utilizes results by PBEþU method. The pressure derivatives of the various Cij except C66 by PBEþ U are correctly described by a linear law. For the longitudinal (C11 and C33) constants, it is noted that C11 oC33 at zero pressure and their pressure derivatives have a similar relationship (i.e. ∂C11/∂P E5.87, ∂C33/∂PE 6.81). The calculated C11 are consistent with the measurements by ISLS, but they are larger than the experimental data by IXS. For the off-diagonal (C12 and C13) constants, the pressure derivative of C12 (4.02 70.03) is within

F.-G. Kuang et al. / Physica B 441 (2014) 72–79

the error range of 3.3 70.94 by ISLS while it is larger than the results of IXL (3.3 70.4). The consequence of ∂C12/∂P 4∂C13/∂P is in line with ISLS measurements, but the IXS measurements show ∂C12/∂P o∂C13/∂P. For the shear elastic modulus (C44 and C66), the shear anisotropy ratio (C44/C66) here (C44/C66 41) is in excellent agreement with ultrasonic measurements and IXS measurements. The pressure derivatives of these two elements have harmonious relationship with the ISLS measurements. In addition, a change in slope near 75 GPa is unobserved of PBEþU method but it is

Fig. 1. Total, Co-4s, and Co-3d density of states for ferromagnetic hcp Co. Dashed line indicates the Fermi level. All curves are represented for spin up and down contributions.

75

showed above 70 GPa by PBE method. The cause of this change is the sudden decreased magnetic moment. However, the magnetic moment by PBE þU method is equal to a constant of 1.70 μB up to 100 GPa. According to Vogit–Reuss–Hill approximation, the obtained elastic properties are presented in Fig. 3(a), including bulk (B), shear (G) and Young's (Y) modulus under compression. Apparently, the bulk modulus, shear modulus and Young's modulus increase monotonously as pressure goes up. At zero pressure, the obtained bulk modulus B of 196 GPa is in good agreement with previous experimental value B ¼199 76 GPa [33]. Our predicted value of Young's modulus at ambient condition is 219 GPa, which is slightly larger than the experimental value of 211 GPa. All of these indicate that the results by PBEþU method show improvements. It is also observed that the pressure has an important influence on the bulk modulus, shear modulus and Young's modulus. The type of comparative study of compressional and shear aggregate sound velocity is essential in providing quantitative estimates for the various proposed mechanisms of the observed elastic anomaly, including magneto-elastic effects of polycrystalline hcp Co. The compressional (vL) and shear (vT) aggregate sound velocities are plotted as a function of density in Fig. 3(b), together with results from ISLS [6], IXS [7] measurements and LAPW-GGA calculations [9]. Directly seen from Fig. 3(b), our calculated results are more in accordance with to the experimental measured values than LAPW calculations. However, the softening of the compressional and shear aggregate sound velocities is not observed. The reason is that the magnetic moments are not collapse up to 100 GPa when Hubbard U correction is considered. In Fig. 3(c), the bulk (vB) and average (vm) aggregate sound velocities as well as the Debye temperature (ΘD) as a function of density are plotted. All of them increase monotonously with pressure increasing. In addition, the result of ΘD at 0 GPa is equal to 447 K that agrees with the reported value of 445 K. By solving the Cristoffel equation for hcp Co, it is obtained that the values of ΔS1 goes away from 1.30, while ΔP tends to approach 1.15 gradually and the value of ΔS2 shows complex changes with pressure increasing from 0 to 100 GPa. The determination of full tensor allows us to derive the complete velocity profile of hcp Co

Fig. 2. The pressure dependence of the elastic modulus Cij.

76

F.-G. Kuang et al. / Physica B 441 (2014) 72–79

Fig. 3. (a) The bulk (B), shear (G) and Young's (Y) modulus of hcp Co as a function of pressure. (b) The compressional (vL) and shear (vT) aggregate sound velocities as a function of density. (c) The bulk (vB) and average (vm) aggregate sound velocities as well as the Debye temperature (ΘD) as a function of density.

ferromagnetic hcp Co. Our results of the sound dispersion yield almost parallel to the IXS measurements with a 4% error at 451. The transverse modes are displayed in Fig. 4(b). The experimentally data and our results of vT1 keep good consistency. The maximum error of  6% at 451 is larger than FAPW-GGA calculations at 35 GPa but our results of 0 GPa agree with them. For vT2, there is a disagreement in the case of ferromagnetic hcp Co. Theoretical calculations predict that it is a sigmoidal shape with a higher sound velocity along the c axis, but the IXS results show a sigmoidal shape with a higher sound velocity along the a axis. A direct cause is that C44 is still larger than C66 at 35 GPa by PBEþU calculations, while the IXS result gives C44 oC66. It should be point out that a consistent conclusion with IXS experiment can be obtained by PBE method. 3.2. Phonon dispersions

Fig. 4. The normalized acoustic sound velocities as a function of the propagation direction in the meridian plane of hcp Co, from c to a axis. (a) quasi-Longitudinal acoustic sound velocity; (b) quasi-Shear and pure shear acoustic sound velocities. IXS measurements are obtained at 36 GPa from Ref. [3]; FAPW-GGA results are computed at 36 GPa form Ref. [9].

crystal. In Fig. 4, the normalized acoustic velocities of the quasilongitudinal wave (vL), quasi-shear wave (vT1), and pure shear wave (vT2) as a function of the angle φ of the propagation direction with respect to the c axis are plotted. The IXS measurements [3] and FAPW-GGA calculations [9] are also shown for comparison. Our results are given at 0 and 35 GPa, while IXS measurements and FLAW-GGA calculations are reported at 36 GPa. From Fig. 4(a), it is observed that there is a sigmoidal shape of vL for

Vibrational properties can be understood in terms of the excitation of the noninteracting phonon. The phonon dispersion curves of ferromagnetic hcp Co along high-symmetry directions at separate volumes are displayed. Seen from Fig. 5(a), the computed results at equilibrium volume are in good agreement with the inelastic neutron scattering data [34]. It is also observed that the phonon dispersions of ferromagnetic hcp Co increase with volume decreasing. The dispersion curves for the applied volumes do not show any anomaly or instability. These conclusions can be also seen from phonon density of states in Fig. 5(b). The elastic properties can be also explained through studying the dynamical behavior of hcp metals. The primitive cell of hcp Co has two atoms, so there are six phonon modes for a given wave vector, three of which are acoustic branches. According to group theory, the optical phonons at Γ-point are belonging to the following irreducible representations: Γopt ¼2E2g þ E1u, where E2g mode corresponding to the transverse optical mode (TO) registers as Raman active and E1u mode corresponding to the longitudinal optical mode (LO) registers as infrared active. The E2g transverse mode is greatly concerned in the Raman experiments. Because the C44 elastic modulus can be calculated from E2g-phonon frequencies using the relation, C44 ¼2π2M[√3c/6a2]ω2 [35], where M is the atomic mass, c and a are the cell parameters, and ω is the E2g-phonon frequency at Γ point. At ambient pressure, our

F.-G. Kuang et al. / Physica B 441 (2014) 72–79

77

Fig. 5. (a) Phonon dispersion curves of ferromagnetic hcp Co at various volumes (V0 is the equilibrium volume); (b) the phonon density of states in the corresponding volume.

Table 2 The equation of state parameters for ferromagentic hcp Co at various pressures. All of our calculations are performed at 300 K, and results of experiments are obtained at room temperature.

This study

Experiment [31] Experiment [33] Experiment [36]

Pressure (GPa)

V0 (Å3)

B0 (GPa)

B00

0 10 20 60 Ambient pressure Ambient pressure Ambient pressure

11.01 10.52 10.14 9.01 11.07 11.10 11.00

198 247 293 450 202 199 199

5.04 4.69 4.40 3.67 3.9 3.6 3.6

Fig. 6. The pressure variation of C44 elastic modulus, along with experimental data (Ref. [6]) and PWSCF calculations (Ref. [10]).

calculated value of 71.7 GPa for C44 is in good agreement with the experimental value 72.9 GPa by ISLS measurements [4] and 70 GPa estimated by plane-wave self-consistent field (PWSCF) method [10]. In Fig. 6, the pressure variation of calculated C44 elastic modulus is presented along with Raman spectroscopy data of Goncharov et al. [6] and theoretical calculations of Modak et al. by PWSCF method [10]. The results of C44 elastic moduli are close to PWSCF results as the pressure is less than 40 GPa and become close to the experimental measurements as the pressure is larger than 40 GPa. The calculated results can be made the linear fitting then a slope change near 75 GPa. The occurred change of slope reflects the stronger magneto-elastic effects of the ferromagnetic hcp Co. 3.3. Thermodynamic properties Thermal properties of solids such as thermal expansion depend on their lattice dynamical behavior. The phonon contribution Fvib as functions of volume V and temperature T can be derived from Eq. (4). Then Gibbs free energies at the corresponding volume and temperature are obtained. Owing to Vinet EOS's versatility and high accuracy, Gibbs free energies of different volumes of single crystal hcp Co are fitted at each temperature. The Vinet EOS

Fig. 7. The fitted Vinet EOS parameters V0(T), B0(T), and B0 0(T) as functions of temperature.

parameters at 300 K are listed in Table 2. Our calculations at 0 GPa and 300 K agree with current experimental measurements [31,35,36] at ambient pressure and room temperature. According

78

F.-G. Kuang et al. / Physica B 441 (2014) 72–79

to the reported literatures, the ferromagnetic state of hcp Co is stable to 100 GPa at ambient temperature [36]. And it remains stable at high temperatures (up to 1000 K) under CO or H2 þCO mixtures [1]. So it is necessary to investigate the thermodynamic properties of hcp Co at higher temperatures and pressures. The room temperature equations of state parameters by fitting Vinet EOS for hcp Co at different pressures are showed in Table 2. The temperature dependences of the equations of state parameters V0(T), B0(T), and B00 (T) are plotted in Fig. 7, which show typical features of transition metals: (a) thermal expansion, decrease of bulk modulus and increase of its pressure derivative with increasing temperature; (b) volume contraction, increase of bulk modulus and decrease of its pressure derivative with increasing pressure.

The other thermal equations of state properties can be derived from Gibbs free energy [37]. The thermal expansion coefficient αV is αV ¼1/V(∂V/∂T)p. A comparison between the calculated and reported values is showed in Fig. 8(a). The calculated αV varies with temperature at 0 GPa agree well with the data obtained by Owen [38]. However, the calculated results are deviation from the values given by Arbuzov [37]. The discrepancy between the calculated thermal expansion coefficients and reported values may be ascribed to the use of the exchange-correlation function. Moreover, at a given pressure, the thermal expansion coefficients increase with temperature, and then converge to a nearly constant value at high temperature. The thermal expansion drop rapidly with increasing pressure (Fig. 8(a)), and this characterize can also

Fig. 8. (a) The thermal expansion coefficient αV as a function of temperature; (b) the Anderson–Grüneisen parameter δT as a function of temperature.

Fig. 9. (a) The Grüneisen ratio γ as a function of temperature; (b) the temperature dependences of the parameter q.

F.-G. Kuang et al. / Physica B 441 (2014) 72–79

observe from the Anderon–Grüneisen parameter δT ¼(∂lnαV/∂lnV)T (Fig. 8(b)). The calculated δT drop rapidly with increasing pressure at a given temperature, and show its complex temperature dependences. It is further shown that the δT is a function of the exponents of the attractive and repulsive terms in Mie's potential energy equation. And it is thus reasonable to expect that it will be independent of temperature. In many cases, the δT is taken as a constant at a fixed pressure. The Grüneisen ratio γ is a most important parameter in understanding the relationship between the thermal and elastic properties, which is particularly for reduction of shock data. The γ can be obtained from: γ ¼V(∂p/∂U)V ¼ αVBTV/CV, where U is the internal energy. Numerous different techniques can determine the Grüneisen ratio γ. By Raman spectroscopy measurements, Goncharov et al. [6] have showed the density dependence of the Grüneisen ratio γ. In this paper, we show the temperature dependence of γ at a given pressure in Fig. 9(a). It is observed that our calculated γ first speedily increase with temperature at low temperature (p¼0, 10 and 20 GPa), and then slowly increase with temperature as T4150 K. By contrast, the γ first drop rapidly at low temperature at p¼ 60 GPa, and then drop slowly as T4150 K. It is found that the pressure and temperature dependence of γ is complex. The volume dependence of the Grüneisen ratio is defined by the parameter q, q¼∂lnγ/∂lnV, and q is usually treated as a constant. However, Fig. 9(b) shows that q strongly depends on both the temperature and pressure, and it even becomes negative at some temperature and pressure. 4. Conclusions In conclusion, the elastic, vibrational and thermodynamic properties of ferromagnetic hcp Cobalt are studied based on first principles projector augmented wave potential method. From the above studies, we have the following conclusions: (a) The PBE þU method can well promote the accuracy of elastic moduli of ferromagnetic hcp Co with regard to experimental data. The bulk modulus, shear modulus, Young's modulus, all the sound velocities and Debye temperature as a function of pressure (or a function of density) are also displayed beneficial improvements. The elastic anisotropies are discussion and the normalized acoustic sound velocities as a function of the propagation direction from the c to the a axis are presented by solving the Cristoffel equation. (b) The computed phonon dispersion curves of hcp Co agree well with experiments at 0 GPa. Under compression, the dispersion curves do not show any anomaly or instability. The C44 elastic modulus from the Raman frequencies of hcp Co, i.e. E2g-phonon frequencies, shows a slope change near 75 GPa. This change also implies the stronger magneto-elastic effects of hcp Co. (c) Under quasi-harmonic approximation, both the calculated thermal expansion coefficient and the Grüneisen ratio are good in agreement with available experimental data and other calculations. The variation in the Grüneisen ratio with volume, i.e. the parameter q, and the Anderson–Grüneisen parameter δT present manifestly temperature and pressure dependences.

79

Acknowledgment This project was supported by the National Natural Science Foundation of China (Grant no. 11274235), the Young Scientists Fund of the National Natural Science Foundational of China (no. 11104190) and the Doctoral Education fund of Education Ministry of China (nos. 20100181110086 and 20110181120112).

References [1] V.A. de la PeñaO'shea, P. Ramírez de. la Piscina, N. Homs, G. Aromí, J.L.G. Fierro, Chem. Mater. 21 (23) (2009) 5637. [2] R.J. Hemley, H.K. Mao, in: G.L. Trigg. (Ed.), Encyclopedia of Applied Physics, VCH Publishers, New York, 1997, p. 555. [3] D. Antonangeli, M. Krisch, G. Fiquet, D.L. Farber, C.M. Aracne, J. Badro, F. Occelli, H. Requardt, Phys. Rev. Lett. 93 (21) (2004) 215505. [4] J.C. Crowhurst, D. Antonangeli, J.M. Brown, A.F. Goncharov, D.L. Farber, C.M. Aracne, Appl. Phys. Lett. 89 (11) (2006) 111920. [5] D. Antonangeli, S. Merkel, D.L. Farber, Geophys. Res. Lett. 33 (24) (2006) L24303. [6] A.F. Goncharov, J. Crowhurst, J.M. Zaug, Phys. Rev. Lett. 92 (11) (2004) 115502. [7] D. Antonangeli, M. Krisch, G. Fiquet, J. Badro, D.L. Farber, A. Bossak, S. Merkel, Phys. Rev. B: Condens. Matter 72 (13) (2005) 134303. [8] G. Steinle-Neumann, L. Stixrude, R.E. Cohen, Phys. Rev. B: Condens. Matter 60 (2) (1999) 791. [9] G. Steinle-Neumann, Phys. Rev. B: Condens. Matter 77 (10) (2008) 104109. [10] P. Modak, A.K. Verma, R.S. Rao, B.K. Godwal, R. Jeanloz, Phys. Rev. B: Condens. Matter 74 (1) (2006) 012103. [11] D.W. Voigt, Lehrbuch der Kristallphysik, Taubner, Leipzig, 1928. [12] R. Hill, Proc. Phys. Soc. London, Sect. A 65 (5) (1952) 349. [13] J.P. Poirier, Introduction to the Physics of the Earth's Interior, Cambridge University Press, Cambridge (2000) 11. [14] O.L. Anderson, J. Phys. Chem. Solids 24 (7) (1963) 909. [15] L.D. Landau, E.M. Lifschitz, Theory of Elasticity, Course of Theoretical Physics, vol. 7, Pergamon Press, New York, 1980 (Chap. 3). [16] P. Vinet, J.H. Rose, J. Ferrante, J.R. Smith, J. Phys. Condens. Matter 1 (11) (1989) 1941. [17] D.M. Ceperley, B.J. Alder, Phys. Rev. Lett. 45 (7) (1980) 566. [18] J.P. Perdew, A. Zunger, Phys. Rev. B: Condens. Matter 23 (10) (1981) 5048. [19] J.P. Perdew, J.A. Chevary, S.H. Vosko, Phys. Rev. B: Condens. Matter 46 (11) (1992) 6671. [20] J.P. Perdew, K. Burke, Y. Wang, Phys. Rev. B: Condens. Matter 54 (23) (1996) 16533. [21] G. Kresse, J. Furthmüller, Comput. Mater. Sci. 6 (1) (1996) 15. [22] G. Kresse, J. Furthmüller, Phys. Rev. B: Condens. Matter 54 (16) (1996) 11169. [23] G. Kresse, J. Hafner, Phys. Rev. B: Condens. Matter 47 (1) (1993) 558. [24] C. Lu, X.Q. Yang, C.Y. Zhu, X.Y. Kuang, Dalton Trans. 41 (32) 9781. [25] C. Loschen, J. Carrasco, K.M. Neyman, F. Illas, Phys. Rev. B: Condens. Matter 75 (3) (2007) 035115. [26] H.J. Kulik, M. Cococcioni, D.A. Scherlis, N. Marzari, Phys. Rev. Lett. 97 (10) (2006) 103001. [27] V.A. de la Peña O'Shea, I. de P.R. Moreira, A. Roldán, F. Illas, J. Chem. Phys. 133 (2) (2010) 024701. [28] A. Togo, L. Chaput, I. Tanaka, G. Hug, Phys. Rev. B: Condens. Matter 81 (17) (2010) 174301. [29] A. Togo, F. Oba, I. Tanaka, Phys. Rev. B: Condens. Matter 78 (13) (2008) 134106. [30] H.J. McSkimin, J. Appl. Phys. 26 (4) (1955) 406. [31] D. Antonangeli, L.R. Benedetti, D.L. Farber, G. Steinle-Neumann, A.L. Auzende, J. Badro, M. Hanfland, M. Krisch, Appl. Phys. Lett. 92 (11) (2008) 111911. [32] C. Kittel, Introduction to Solid State Physics, Wiley, New York, 2005. [33] H. Fujihisa, K. Takemura, Phys. Rev. B: Condens. Matter 54 (1) (1996) 5. [34] N. Wakabayashi, R.H. Scherm, H.G. Smith, Phys. Rev. B: Condens. Matter 25 (8) (1982) 5122. [35] J.C. Upadhyaya, D.K. Sharma, D. Prakash, S.C. Upadhyaya, Can. J. Phys. 72 (1–2) (1994) 61. [36] C.S. Yoo, H. Cynn, P. Soderlind, V. Iota, Phys. Rev. Lett. 84 (18) (2000) 4132. [37] O.L. Anderson, Equations of State of Solids for Geophysics and Ceramic Science, Oxford University Press, New York, 1995. [38] E.A. Owen, D.M. Jones, Proc. Phys. Soc. London, Sect. B 67 (6) (1954) 456.