First principle study of elastic and thermodynamic ...

First principle study of elastic and thermodynamic properties of FeB4 under high pressure Xinyu Zhang, Jiaqian Qin, Jinliang Ning, Xiaowei Sun, Xinting Li et al. Citation: J. Appl. Phys. 114, 183517 (2013); doi: 10.1063/1.4829926 View online: http://dx.doi.org/10.1063/1.4829926 View Table of Contents: http://jap.aip.org/resource/1/JAPIAU/v114/i18 Published by the AIP Publishing LLC.

Additional information on J. Appl. Phys. Journal Homepage: http://jap.aip.org/ Journal Information: http://jap.aip.org/about/about_the_journal Top downloads: http://jap.aip.org/features/most_downloaded Information for Authors: http://jap.aip.org/authors

JOURNAL OF APPLIED PHYSICS 114, 183517 (2013)

First principle study of elastic and thermodynamic properties of FeB4 under high pressure Xinyu Zhang,1,a) Jiaqian Qin,2,1,a) Jinliang Ning,1 Xiaowei Sun,1 Xinting Li,1 Mingzhen Ma,1 and Riping Liu1,a)

1 State Key Laboratory of Metastable Materials Science and Technology, Yanshan University, Qinhuangdao 066004, People’s Republic of China 2 Metallurgy and Materials Science Research Institute, Chulalongkorn University, Bangkok 10330, Thailand

(Received 23 September 2013; accepted 28 October 2013; published online 13 November 2013) The elastic properties, elastic anisotropy, and thermodynamic properties of the lately synthesized orthorhombic FeB4 at high pressures are investigated using first-principles density functional calculations. The calculated equilibrium parameters are in good agreement with the available experimental and theoretical data. The obtained normalized volume dependence of high pressure is consistent with the previous experimental data investigated using high-pressure synchrotron x-ray diffraction. The complete elastic tensors and crystal anisotropies of the FeB4 are also determined in the pressure range of 0–100 GPa. By the elastic stability criteria and vibrational frequencies, it is predicted that the orthorhombic FeB4 is stable up to 100 GPa. In addition, the calculated B/G ratio reveals that FeB4 possesses brittle nature in the range of pressure from 0 to 100 GPa. The calculated elastic anisotropic factors suggest that FeB4 is elastically anisotropic. By using quasi-harmonic Debye model, the compressibility, bulk modulus, the coefficient of thermal expansion, the heat capacity, and the Gr€uneisen parameter of FeB4 are successfully obtained in the C 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4829926] present work. V

I. INTRODUCTION

Superhard materials are useful in a variety of industrial applications, such as abrasives, cutting tools, and coatings. So great effort is currently focused on the synthesis and characterization of novel superhard materials exhibiting simultaneously very low compressibility, wide thermodynamic ranges of chemical stability, and high scratch resistance as well as surface durability.1–16 Traditionally, it is commonly accepted that superhard materials (diamond, cBN, BC2N, etc.) are those strongly covalent bonded compounds formed by light elements (LE), namely, B, C, N, and O. Recent attempts to design new intrinsically superhard materials concentrated on the introduction of light elements forming strong bonds (B, C, N, and O) into transition metal (TM) with high elastic moduli.4–8,11,12,14–16 Several TM-LE compounds have been synthesized (e.g., OsB2,4 ReB2,6 PtN2,17 etc.). However, further calculation and experiment indicate they are not superhard materials.8,16 Recently, Gu et al.7 synthesized WB4 and the measured average hardness is very high exceeding 46.2 GPa, comparable to that of cBN (45–50 GPa). Niu et al.12 synthesized CrB4, and predicted the Vickers hardness of 48 GPa. Very recently, Gou et al.14 synthesized iron tetraboride FeB4 and measured the nanoindentation hardness as 65(5) GPa. Zhang et al.15 calculated the elastic moduli and ideal strength of the new synthesized FeB4, which reveal that FeB4 is a hard material. Although there are some studies, still the mechanical properties of FeB4 under high pressure and thermodynamic properties are limited both experimentally and theoretically as far as we are a)

Authors to whom correspondence should be addressed. Electronic addresses: [email protected]; [email protected]; and [email protected]

0021-8979/2013/114(18)/183517/7/$30.00

aware. Especially, the elastic constants and anisotropy under pressure are very important to determine the response of the crystal to external forces, as characterized by the bulk and shear module and can provide a deeper insight into mechanical behavior and hardness of materials. Moreover it is essential for many practical applications, such as load deflection, themoelastic stress, fracture toughness, anisotropic character of the bonding, and structural stability. On the other hand, no report of the thermodynamic properties for the new synthesized FeB4 has been found in literatures. Therefore, in this paper, the elastic and thermodynamic properties of FeB4 under pressure up to 100 GPa are investigated using first-principles density functional calculations. The elastic properties of orthorhombic FeB4 under high pressure are studied for the first time, from which the elastic anisotropy are also determined. To further study the FeB4, the thermodynamic properties, such as the heat capacity, thermal expansion, velocity, Gr€unesisen parameters, and so on are estimated by the quasi-harmonic Debye model. II. COMPUTATIONAL METHODS

In the present work, our first-principles calculations on FeB4 are performed with the VASP code18 using the plane-wave pseudopotential method and employing the generalized gradient approximation (GGA).19 The all-electron projector augmented wave (PAW) method20 was employed with a plane-wave cutoff energy of 600 eV. The k-point samplings with 9  8  13 in the Brillouin zone were performed using the Monkhorst-Pack scheme.21 The total energy convergence tests showed that convergence to within 1 meV/atom was achieved with the above calculation parameters. Single crystal elastic constants were calculated via a strain-energy approach,

114, 183517-1

C 2013 AIP Publishing LLC V

183517-2

Zhang et al.

J. Appl. Phys. 114, 183517 (2013)

i.e., by applying a small strain to the equilibrium lattice of orthorhombic unit cell and fitting the dependence of the resulting change in energy on the strain.22 The nine distortion types necessary for computing the nine elastic constants of an orthorhombic structure are given in Table I, in which V0 is the equilibrium volume. All internal atomic coordinates were fully optimized for every given strain state. The bulk modulus, shear modulus, Young’s modulus, and Poisson’s ratio were determined by using the Voigt-Reuss-Hill approximation.23 The phonon calculations were carried out by using a supercell approach as implemented in the PHONOPY code.24 The Debye temperature can be estimated from the elastic constants using the average sound velocity  m, by the following equation:25    1 h 3n qNA 3 H¼ m ; (1) k 4p M where h is Planck’s constant, k is Boltzmann’s constant, NA is Avogadro’s number, n is the number of atoms per formula unit, M is the molecular mass per formula unit, and q is the density. The average sound velocity  m is given by   1 1 2 1 3 m ¼ þ ; (2) 3  3t  3t where  t and  l are the transverse and longitudinal elastic wave velocity of the polycrystalline materials and are given by Navier’s equation, sffiffiffiffi G ; (3) t ¼ q sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3B þ 4G l ¼ ; (4) 3q where q is the density. TABLE I. Deformation types, expressed in the Voigt notation (xx ¼ 1, yy ¼ 2, zz ¼ 3, yz ¼ 4, xz ¼ 5, and xy ¼ 6) that are used to calculate the nine independent single-crystal elastic constants of an orthorhombic crystal. Parameters (unlisted ei ¼ 0)

Relationship between DE=V0 and elastic constant

1

e1 ¼ d

2

e2 ¼ d

3

e3 ¼ d

4

e4 ¼ d

5

e5 ¼ d

6

e6 ¼ d

7

e1 ¼ d; e2 ¼ d

DE 1 ¼ C11 d2 V0 2 DE 1 ¼ C22 d2 V0 2 DE 1 ¼ C33 d2 V0 2 DE 1 ¼ C44 d2 V0 2 DE 1 ¼ C55 d2 V0 2 DE 1 ¼ C66 d2 V0 2   DE C11 C22 2 ¼ þ C12 þ d V0 2 2   DE C22 C33 2 ¼ þ C23 þ d V0 2 2   DE C11 C33 2 ¼ þ C13 þ d V0 2 2

Strain

8

e2 ¼ d; e3 ¼ d

9

e1 ¼ d; e3 ¼ d

In order to investigate the thermodynamic properties, we performed the quasi-harmonic Debye model,26 in which the phononic effect is considered. However, it should be noted that the Debye model of phonon density of state is essentially a linear extrapolation of the sound speed of acoustic branches but not optical branches in order to get density states. This model is very quick and easy, and it has been successfully applied to predict the thermodynamic properties of some materials.27–29 The nonequilibrium Gibbs function G*(V; P, T) takes the form of G ðV; P; TÞ ¼ EðVÞ þ PV þ AVib ðHðVÞ; TÞ;

(5)

where HðVÞ is the Debye temperature, and the vibrational term AVib can be written as   9H þ 3 ln ð1  eH=T Þ  DðH=TÞ ; (6) AVib ðH; TÞ ¼ nKB T 8T here DðH=TÞ represents the Debye integral, n is the number of atoms per formula unit. For an isotropic solid, H is expressed as rffiffiffiffiffiffiffiffiffi BS H ¼ hð6p2 V 1=2 nÞ1=3 f ðrÞ 2 ; (7) k M where V, M, n, and f ðrÞ is the molar volume, molar mass, the number of atom per formula unit, and a scaling function that depend on Poisson’s ratio of the isotropic solid, respectively. The BS is the adiabatic bulk modulus, which is equal to the isothermal bulk modulus BT in the quasi-harmonic Debye model, leading to the following equation:   d2 E : (8) BS ¼ BT ¼ V dV 2 Here E is the total energy of the crystal at 0 K. The non-equilibrium Gibbs function G ðV; P; TÞ as a function of (V; P, T) can be minimized with respect to volume V,    @G ðV; P; TÞ ¼ 0: (9) @V P;T By solving Eq. (9), one can get the thermal equation of state (EOS) V (P, T). The isothermal bulk modulus BT , heat capacity CV, and the thermal expansion coefficient a are given by   @ 2 G ðV; P; TÞ ; (10) BT ðP; TÞ ¼ V @V 2 P;T   3H=T ; (11) C ¼ 3nk 4DðH=TÞ  H=T 1 e CP ¼ CV ð1 þ a c TÞ; a¼

c CV ; BT V

(12) (13)

where c is the Gr€uneisen parameter defined as c¼

d ln HðVÞ : d ln V

(14)

183517-3

Zhang et al.

J. Appl. Phys. 114, 183517 (2013)

FIG. 1. The total energy per formula unit as a function of volume for FeB4, and the inset shows the crystal structure of FeB4.

FIG. 2. The normalized lattice parameters and volume as a function of pressure.

III. RESULTS AND DISCUSSION

ambient conditions and B-Fe bonding along b-direction that has been clarified in our previous work.15

A. The equilibrium lattice structure and equation of state

To obtain equilibrium lattice parameters for the orthorhombic FeB4, we first make the total energy electronic structure calculations by the first principles method. Figure 1 shows the total energy per formula unit as a function of volume for FeB4, and the inset shows the crystal structure of FeB4. The total energy is calculated by varying the volume for orthorhombic FeB4. The calculated E–V data are fitted to the third-order Birch-Murnaghan EOS, and the calculated equilibrium structure parameters, bulk modulus, and its pressure derivative are given in Table II together with the experimental results for comparison. The calculated lattice parameters are in well agreement with the experimental results, and the mismatch of lattice parameter is within 0.02%. Furthermore, the pressure dependence of the normalized parameters a/a0, b/b0, c/c0 and volume V/V0 as a function of pressure for the FeB4 is plotted in Fig. 2, where a0, b0, c0 and V0 are the values at T ¼ 0 K and P ¼ 0 GPa. From Fig. 2, it is clear that the structure of FeB4 is most incompressible along the b-direction, while softest along the a-axis. The calculated results agree well with the high pressure X-ray diffraction results14 and the ideal tensile results.15 This indicates the clear elastic anisotropy of the FeB4 crystal. The smallest compression along the b-direction might be ˚ ) at attributed to the strong covalent B-B bonding (1.707 A ˚ ) and TABLE II. The calculated equilibrium lattice constants a0, b0, c0 (A ˚ 3), EOS fitted bulk modulus B0 (GPa), and its presequilibrium volume V0 (A sure derivative B00 for the orthorhombic FeB4 at 0 K and 0 GPa.

This work Experimentala Theoreticalb a

Reference 14. Reference 30.

b

a0

b0

c0

V0

B0

B00

4.524 4.5786 4.521

5.284 5.2981 5.284

3.004 2.9991 3.006

71.810 72.752 71.810

266 252

4.1 3.5

B. Elastic properties

The nine calculated elastic stiffness coefficients (Cij) are listed in Table III, along with the theoretical bulk modulus, shear modulus and Young’s modulus at 0 GPa.31 It can be seen clearly that the present data are in excellent agreement with the data reported by Bialon et al.31 When compared with the latest synthesized CrB4, which can be synthesized under ambient pressure,12 the obtained bulk modulus (B ¼ 277 GPa) of FeB4 is larger than that of CrB4 (B ¼ 265 GPa), whereas the shear modulus (G ¼ 186 GPa) of FeB4 is lower than that of CrB4 (G ¼ 261 GPa). Because the shear modulus and Pugh’s ratio32 (k ¼ G/B ¼ 0.671) are thought of as two important elastic properties which strongly correlated to hardness,12 the lower shear modulus and Pugh’s ratio reveal that FeB4 cannot be harder than CrB4. The elastic constants of FeB4 under high pressure up to 100 GPa are plotted in Fig. 3. Unfortunately, there are not yet experimental data or theoretical calculations to compare with our predicted results for the pressure derivative of elastic properties. Therefore, we believe that, our results could serve as a prediction for future studies. It can be seen that the elastic constants Cij increase monotonically with the increasing applied pressure P. Moreover, these nine independent elastic constants Cij still satisfy the well known Born stability criteria up to 100 GPa, indicating that FeB4 is still mechanical stable at high pressure of 100 GPa. As shown in Table III, the calculated elastic moduli also increase monotonically with increasing pressure as expected. The ratio between the bulk and the shear modulus B/G is used to predict the brittle or ductile behavior of materials. Here, we analyzed the ductility and brittle nature at pressure according to the Pugh criterion.32 The ductile behavior is predicted when B/G > 1.75, otherwise the materials behaves in a brittle manner. According to Table III, the ratio of B/G increases with pressure and reaches 1.61 at 100 GPa. The results indicate that the FeB4 is prone to brittleness in the presently studied

183517-4

Zhang et al.

J. Appl. Phys. 114, 183517 (2013)

TABLE III. The elastic constants Cij (GPa), bulk modulus B (GPa), shear modulus G (GPa), Young’s modulus E (GPa), and the B/G ratio of FeB4 under pressure. P

C11

C22

C33

C12

C23

C13

C44

C55

C66

B

G

E

B/G

0

408

754

448

165

154

160

219

141

229

491 529 597 610 655 665 743 796 844 862

854 890 965 1027 1054 1099 1172 1243 1289 1341

521 554 602 652 664 710 749 786 830 857

177 194 238 250 282 282 325 352 368 376

195 223 269 291 343 358 389 416 437 458

194 220 268 288 332 353 402 440 482 504

246 264 281 298 312 328 342 356 370 383

180 203 222 241 252 267 282 295 308 320

261 284 307 327 346 364 382 399 415 431

186 187a 217 233 247 262 266 279 292 304 315 324

456 457a 534 572 616 653 670 702 740 774 805 829

1.49

10 20 30 40 50 60 70 80 90 100

277 274a 325 353 405 430 468 487 536 574 608 630

1.50 1.52 1.63 1.64 1.76 1.75 1.84 1.89 1.93 1.95

a

Reference 31.

pressure range. Moreover, the partial and total densities of states (DOSs) of FeB4 at 0 GPa, 50 GPa, 100 GPa are calculated as shown in Fig. 4 to examine the evolution of electronic structure of FeB4 under pressures. It can be seen that the DOSs at 0 GPa agree well with previous calculation30 and they show similar profiles under pressures, which indicates the electronic stability of FeB4 up to 100 GPa. Additionally, a stable crystalline structure requires all phonon frequencies to be positive. Figure 5 shows the full phonon dispersion curves of FeB4 at 50 GPa and 100 GPa. As shown in Fig. 5, no imaginary phonon frequency was found in the whole Brillouin zone, indicating the dynamical stabilities of FeB4 up to 100 GPa. C. Elastic anisotropy

It is known that superhard materials should preferably be isotropic, otherwise it would deform preferentially in a given direction.33 That is to say, microcracks may be induced in materials due to the significant elastic anisotropy. Furthermore, our previous ideal tensile strength15 and also the present pressure dependence of the normalized

FIG. 3. Pressure dependence of the elastic stiffness coefficients (Cij) of FeB4 at 0 K.

parameters show clear elastic anisotropy of FeB4. Hence it is important to study elastic anisotropy to shed light on their mechanical durability. Elasticity describes the response of a crystal under external strain and provides key information about the bonding characteristics between adjacent atomic planes and the anisotropic character of the solid.34 The shear anisotropic factors provide a measure of the degree of anisotropy in the bonding between atoms in different planes. The shear anisotropic factor for the {100} shear planes between the h011i and h010i directions is A1 ¼

4c44 ; c11 þ c33  2c13

(15)

for the {010} shear planes between the h101i and h001i directions is A2 ¼

4c55 ; c22 þ c33  2c23

(16)

and for the {001} shear planes between the h110i and h010i directions is

FIG. 4. Total and partial density of states for FeB4 at 0 GPa, 50 GPa, and 100 GPa.

183517-5

Zhang et al.

J. Appl. Phys. 114, 183517 (2013)

FIG. 5. Phonon dispersion curves for FeB4 at 50 GPa and 100 GPa.

A3 ¼

4c66 : c11 þ c22  2c12

(17)

For an isotropic crystal the factors A1, A2, and A3 must be 1.0, while the deviation from one is a measure of the degree of the elastic anisotropy. In addition, since the FeB4 is orthorhombic, not cubic, the shear anisotropic factors are not sufficient to describe the elastic anisotropy. Thus, the anisotropy of the linear bulk modulus is also considered. The directional bulk modulus along different crystallographic axis can be defined as Bi ¼ i(dP/di)(i ¼ a, b, and c).34 Using the relations mentioned above, the parameters about elastic anisotropy are estimated and presented in Table IV. It is seen that FeB4 is elastic anisotropic. At 0 GPa, the shear anisotropy results indicate that the elastic anisotropy for the {100} shear planes between the h011i and h010i directions is close to that of the {010} shear planes between the h101i and h001i directions, but larger that of the {001} shear planes between the h110i and h010i directions. Furthermore, the shear anisotropic factor for the {001} shear planes between the h110i and h010i directions is closer to 1.0 than other two factors, which can be explained by the unique corrugated B6 units ring along [100] directions.15 These results also reveal TABLE IV. The shear anisotropy factors A1, A2, A3 and elastic anisotropy index Au and the directional bulk modulus Ba, Bb, and Bc of FeB4 under pressure. P 0 10 20 30 40 50 60 70 80 90 100

A1

A2

A3

Au

Ba

Bb

Bc

1.636 1.577 1.642 1.695 1.738 1.905 1.961 1.988 2.028 2.085 2.155

0.631 0.731 0.814 0.863 0.879 0.977 0.977 0.987 0.986 0.990 0.998

1.099 1.053 1.102 1.131 1.150 1.209 1.213 1.208 1.196 1.188 1.188

0.451 0.349 0.339 0.350 0.375 0.451 0.483 0.492 0.522 0.551 0.604

617.9 744.1 816.0 966.0 964.9 1082.0 1069.5 1265.3 1390.4 1498.6 1514.5

1466.6 1608.9 1709.0 1961.3 2097.1 2299.5 2304.7 2494.8 2656.3 2665.5 2774.2

705.3 842.9 923.5 1028.4 1158.6 1220.6 1379.7 1423.0 1491.7 1613.5 1706.8

that the {100} and {010} shear planes are easier to be the cleavage planes among these principal planes. Furthermore, the lowest ideal shear strength (100) h001i (38.3 GPa) could be explained using this highest anisotropy of {100} shear planes between the h011i and h010i directions. Moreover, it can be seen that the A1 and A2 increase with increasing pressure, and A3 increases in the pressure range 0–20 GPa, and then remains nearly invariant. Meanwhile, the directional bulk modulus Bb along the b axis is the largest when compared to the Ba and Bc, which is consistent with the pressure dependence of the normalized lattice parameters (Fig. 2) and our previous ideal tensile strength.15 In addition, the universal elastic anisotropy index AU is defined by Ranganathan and Ostoja-Starzewski from the bulk modulus B and shear modulus G by Voight and Reuss approaches, as35 AU ¼ 5

GV BV þ  6; GR BR

(18)

where AU ¼ 0 is for isotropic materials. The calculated AU of FeB4 under high pressure is also summarized in Table IV. At 0 GPa, AU ¼ 0.451, it indicates that FeB4 is anisotropic materials, and AU increases with increasing pressure. D. Thermodynamic properties

Through the quasiharmonic Debye model, numerous thermodynamic quantities of FeB4 at various temperatures and pressures are obtained from the calculated energy-volume points at 0 K and 0 GPa. The calculated normalized volume V/V0 and isothermal bulk modulus B as a function of temperature up to 1200 K at various pressures are plotted in Fig. 6. It can be seen that the V/V0 curve becomes steeper with increasing temperature, which indicates that FeB4 is compressed more easily at higher temperatures. It is also found that the bulk modulus decreases slightly and linearly with temperature at given pressure and increases with pressure at a given temperature, which is consistent with the trend of volume. Meanwhile, pressure has a more significant effect on B than temperature.

183517-6

Zhang et al.

J. Appl. Phys. 114, 183517 (2013)

FIG. 6. The calculated normalized volume V/V0 and bulk modulus of FeB4 as a function of pressure at temperatures 0, 300, 600, 900, and 1200 K.

Debye temperature is a fundamental parameter of a material, which is linked to many physical properties, such as specific heat, elastic constants, and melting temperature. One of the standard methods to calculate the Debye temperature is from the elastic constants. The Debye temperature can be estimated from the averaged sound velocity,  m, by Eq. (1). At T ¼ 0 and P ¼ 0, we calculated the Debye temperature HD ¼ 1089 K, which is lower than 1336 K of MnB4,36 but higher than that of other transition metal borides (TcB4, HD ¼ 1050 K,36 ReB4, HD ¼ 824 K,36 ReB2, HD ¼ 744 K,37 OsB2, HD ¼ 780 K,38 OsB4, HD ¼ 781 K (Ref. 27)). The calculated Debye temperatures as well as the compressional velocity and the shear wave velocity obtained under different pressures are presented in Table V. It is shown that the Debye temperature increases with increasing pressure. As one of the most important thermodynamic parameters, the heat capacity CV of a substance not only provides essential information about its vibrational properties but also is fundamental to many applications. In Fig. 7, the temperature dependence of the heat capacity CV at various pressures is presented. It can be seen that as temperature increases, CV increases, while when pressure increases, CV decreases. As shown furthermore, at low temperatures (