First-principles study on the lattice dynamics of FeSb2

Solid State Communications 152 (2012) 231–234

Contents lists available at SciVerse ScienceDirect

Solid State Communications journal homepage: www.elsevier.com/locate/ssc

First-principles study on the lattice dynamics of FeSb2 Rende Miao a,∗ , Guiqin Huang b , Chunhui Fan c , Zhong Bai a , Yanbiao Li a , Liang Wang a , Li_an Chen d , Wenguang Song e , Qiangui Xu e a

Institute of Science, PLA University of Science and Technology, Nanjing 211101, China

b

Department of Physics and Institute of Theoretical Physics, Nanjing Normal University, Nanjing 210046, China

c

Institute of Meteorology, PLA University of Science and Technology, Nanjing 211101, China

d

School of Physics and Electronic Engineering, Nanjing Xiao Zhuang College, Nanjing 211171, China

e

Department of Basic Disciplines, International Studies University of PLA, Nanjing 210039, China

article

info

Article history: Received 7 January 2011 Received in revised form 14 October 2011 Accepted 19 October 2011 by F. Peeters Available online 25 October 2011 Keywords: A. FeSb2 D. Phonons D. Electron–phonon interactions

abstract The lattice dynamics of FeSb2 are investigated by first-principles calculations based on a plane wave pseudopotential method. Phonon spectrum and the electron–phonon linewidths at the Γ -point are obtained using the density-functional perturbation theory within the linear response method. Nine phonon modes are in good agreement with the experimental data, but the B3g mode shows unusually large disagreement. In order to investigate the possibility of anharmonicity of phonon modes, frozen phonon calculations have been performed for the B11g , B21g and B3g modes. But the results are all equal to the density-functional perturbation theory calculations, indicating that the phonons are all harmonic. Our calculated electron–phonon linewidth of the B11g mode is consistent with experimental data, and we also confirm the existence of electron–phonon interactions for the A1g , B11g , and B1u modes in FeSb2 . Published by Elsevier Ltd

1. Introduction Iron diantimonide FeSb2 is a narrow-gap semiconductor, which recently attracted much attention both for fundamental physics and potential applications. FeSb2 represents a number of interesting and extraordinary behaviors, such as unusual magnetic properties [1], anisotropic transport properties [1–4], and the noticeably thermoelectric properties (colossal Seebeck coefficient about −45 mV K−1 at 10 K which is the largest power factor ever reported), which suggest their potential application as a future solidstate thermoelectric cooling device at cryogenic temperatures [5]. Doping into FeSb2 with Te or Sn (electrons or holes doping) metalizes the compounds [6–8], and reduces the Seebeck coefficient [6, 8]. Ab initio local density approximation plus U electronic calculations found that the ground state of FeSb2 could be close to a ferromagnetic instability [9]. On the other hand, unpolarized and polarized Raman scattering measurements of FeSb2 have been reported [10–13]. Early, Lutz and Müller observed two Raman active modes by unpolarized Raman scattering measurements [10]. Latterly, using polarized Raman scattering measurements, Racu et al.

∗

Corresponding author. E-mail address: [email protected] (R. Miao).

0038-1098/$ – see front matter. Published by Elsevier Ltd doi:10.1016/j.ssc.2011.10.022

found three Raman active modes, and suggested that the temperature dependence of linewidths A1g of and the two B1g Raman active modes come from phonon–phonon (ph–ph) interactions, while the infrared (IR) spectroscopy showed an electron–phonon (e–ph) interaction for B1u phonon mode [11]. All six Raman active modes have not been observed until the recent studies by Lazarević et al. [12]. Contrary to Racu et al., the temperature dependence of energy and linewidth of the highest energy B11g mode was believed to come from a strong e–ph interaction for T < 40 K and ph–ph interaction for T > 40 K, respectively [13]. Very recently, Tomczak et al. found that, using GW approximation, the insulating ground state of FeSb2 could be well captured, but the low-temperature Seebeck coefficient of FeSb2 is incompatible with a local electronic picture, suggesting the importance of vertex corrections and nonlocal self-energy effects, or the presence of a substantial phonondrag effect [14]. However, to our best knowledge, there is still a lack of theory study on the lattice dynamics of FeSb2 . In this paper, we report a first-principles investigation of the lattice dynamics of FeSb2 . Density-functional perturbation theory (DFPT) was employed to evaluate the phonon spectrum, and e–ph coupling strengths at the Γ -point. The possibility of anharmonicity of phonon modes was probed by a frozen phonon approach. A careful comparison with available experimental data shows that our density-functional theory (DFT) calculations of phonon frequencies are in reasonable agreement with experimental data

R. Miao et al. / Solid State Communications 152 (2012) 231–234

except for B3g mode. We also provide support for the existence of e–ph interactions for the A1g , B11g , and B1u modes. The rest of paper has been organized as follows. In Section 2, we briefly describe details of the methods used in the present study. Section 3 presents the results, including the structural relaxation, bulk modulus, electronic band structure, phonon spectra, the zone-center phonon frequencies, as well as e–ph linewidths and coupling strengths. Section 4 concludes the paper. 2. Computational details

4 3 2 Energy [eV]

232

1 EF

0 -1

3. Results and discussion FeSb2 crystallizes in the orthorhombic crystal structure with two formula units per unit cell. The space group is Pnnm (D12 2h ) (No. 58), with Fe atoms occupying the 2a Wyckoff positions at (0, 0, 0), and Sb atoms occupying 4g Wyckoff positions at (u, v, 0), where u and v are internal parameters. The optimized structural parameters a = 5.86 Å, b = 6.60 Å, c = 3.17 Å with (u = 0.188, v = 0.355) are in good agreement with the experimental data a = 5.82117 Å, b = 6.50987 Å, c = 3.16707 Å (T = 20 K) with (u = 0.1875, v = 0.3554) [2], the discrepancies are about 0.7%, 1.4%, and 0.09% for a, b and c, respectively, and the equilibrium volume V0 = 73.75 cm3 /mol is about 2.2% larger than experimental data. We also optimized structural parameters under pressure to obtain volume–pressure relation. Then the bulk modulus was calculated from the expression as B0 = (−V0 dP /dV )P =0 , we obtain B0 = 100.1 GPa which is somewhat larger than the experimental data of 84 GPa [2], and in agreement with the theoretical value of 94 GPa [22]. With the equilibrium lattice parameters, we calculated the electronic band structure of FeSb2 , and the results are shown in Fig. 1. As can be seen from Fig. 1, there are doubly degenerate bands crossing the Fermi surface with a small hole pocket at R-point, indicating the weak metallicity of FeSb2 in our DFT—GGA calculation. FeSb2 represents anisotropic transport properties. The electrical resistivity along the a and b axes of FeSb2 shows semiconducting behavior with a rapid increase for T < 100 K, while along the c axis, the resistivity exhibits a metal to semiconductor transition at about 40 K [1,3], so FeSb2 is not strictly semiconducting. Since our DFT

-2 -3

Γ

X

S

Y

Γ

Z

T

R

U

Z

Fig. 1. Band structure of FeSb2 . The Fermi energy is set to zero.

300 250 Frequency (cm–1)

We use DFT calculations [15] as implemented in ABINIT package [16]. The all-electron potentials were replaced by normconserving pseudopotentials as generated in the scheme of Troullier–Martins [17] with Fe (3d6 , 4s2 ), and Sb (5s2 , 5p3 ) levels treated as valence states. The electronic exchange and correlation effects were treated through the generalized gradient approximation (GGA) in the formalism of Perdew–Burke–Ernzerhof (PBE) [18]. The plane wave cutoff was set at 50 Hartree and a 8 × 8 × 8 Monkhorst–Pack grid [19] was used to sample the Brillouin zone (BZ) in the self-consistent calculations. Structural optimization has been performed using the Broyden–Fletcher–Goldfarb–Shanno minimization [20], modified to take into account the total energy as well as the gradients. The symmetry was maintained during optimization, and the relaxation was terminated when the residual force on each atom was less than 5.0d−5 Hartree/Bohr. Dynamical matrices and e–ph linewidths have been computed within DFPT [21] in the harmonic approximation. In order to obtain full phonon spectrum, first the dynamical matrix was calculated on a 4 × 4 × 4 mesh of q-points, then a Fourier interpolation was performed. Phonon density of states (DOS) and partial density of states (PDOS) were calculated using the Gaussian method. The e–ph linewidths at the Γ -point were calculated with a denser (16 × 16 × 16) mesh. All these parameters have been chosen so as to ensure that the values of lattice parameters and phonon frequencies are well converged.

200 150 100 50 0

Γ

X

S

Y

Γ

Z

T

R

U

Z

Fig. 2. Calculated phonon dispersion relations of FeSb2 .

calculations were performed at zero temperature, the weak metallicity seems to be justifiable. The overall band profiles are found to be in reasonable accordance with previous DFT calculations [9,14], showing the validity of our model. Accurately capturing the band gap of the semiconducting FeSb2 requires the inclusion of many-body effects, which can be found in Ref. [14]. Since this paper is focused on lattice-dynamical properties, which are electronic ground-state properties, rather than electronic excited-state properties, we do not discuss the band structure in further detail. In the following, we will present a comprehensive study on the latticedynamical properties of FeSb2 . The phonon-dispersion curves along high symmetry lines of FeSb2 are presented in Fig. 2. There are six atoms in the unit cell which give rise to 18 phonon bands. The absence of imaginary frequencies in the full phonon spectra indicates the dynamical stability of the orthorhombic crystal structure of FeSb2 . The phonon DOS and PDOS are shown in Fig. 3. As can be seen from the phonon PDOS, there is a pseudogap at ∼190 cm−1 , above which is mainly characterized by vibrations of modes of Fe, while below which is predominantly composed of vibrational modes of Sb. The theoretical group analysis predicts the following irreducible representations of vibrational zone-center modes:

Γ = 2Ag + 2Au + 2B1g + 2B1u + B2g + 4B2u + B3g + 4B3u . The 2Ag + 2B1g + B2g + B3g modes are Raman active, while the 2B1u + 4B2u + 4B3u modes are IR-active, and the two Au modes are silent. The calculated frequencies as well as experimental measurements of FeSb2 [10–12,23] are summarized in Table 1. After subtracting the acoustic modes, the seven IR-active modes are B1u + 3B2u + 3B3u , however, in Ref. [23], the amount of B2u modes is four, and the total number of B1u and B3u modes is three. As for the B2u modes, the resonance frequencies of the modes at

R. Miao et al. / Solid State Communications 152 (2012) 231–234

233

Table 1 The calculated phonon frequencies (in cm−1 ), e–ph linewidths γ e–ph (in cm−1 ), and e–ph coupling strengths λ at the Γ -point. Mode

Symmetry

Raman

Exp. [10]

Exp. [12]

Exp. [23]

Calculations

5K

300 K

Room temperature

10 K

Frequencies

γ e–ph

λ

166

157

153.6

162.8

0.15

0.015

150.7

157.1

0.082

0.0093

B11g

190

180

173.9

187.6

0.41

0.030

B21g

160

150

154.3

185.0

0.081

0.0066

B2g

90.4

91.7

0.0035

0.0013

B3g

151.7

96.6

0.0067

0.0023

A1g

175

A2g

154

195a

B1u

IR

Silent a

Exp. [11]

198a

191.3

0.14

0.011

B12u

271

298.7

0.012

0.00045

B22u

231

235.3

0.089

0.0045

B32u

106.4

116.3

0.0019

0.00048

B13u

261.4

255.2

0.0060

0.00025

B23u

216

240.7

0.11

0.0053

B33u

121

126.3

0.0072

0.0013

A1u

203.1

0.013

0.00097

A2u

87.6

0.0011

0.00043

Extracted from Fig. 4 of Ref. [11].

0.30

g(E)(states/cm-1)

0.25 0.20 0.15 0.10 0.05 0.00 0

50

100

150

200

250

300

E (cm-1) Fig. 3. Calculated phonon DOS and PDOS of FeSb2 . The black solid line represents total phonon DOS, while the red and blue dashed lines correspond to Fe and Sb PDOS, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

257 and 271 cm−1 are quite close to each other (see Ref. [23]), so we suppose they are only one B2u mode for the present. For the B1u and B3u modes, since B1u ≃ 195 cm−1 have been assigned in Ref. [11], we suppose the observed 121, 216, 261.4 cm−1 modes are all B3u modes. Following the above presumption, our calculated phonon frequencies of B1u , B22u , B32u , B13u and B33u modes are in good agreement with the experimental data of Ref. [23], while B12u and B23u modes are overestimated by about 10%. For the Raman active modes, the calculated frequency of A1g mode is closer to Refs. [11,12] than Ref. [10], so our result supports the A1g mode assignment of Refs. [11,12]. Calculated

frequencies of A2g , B11g , and B2g modes are in good agreement with the experimental data, while B21g modes are overestimated by about 10%. However, the B3g mode is underestimated by more than 30%. It is quite unusual to give such a large discrepancy between calculated values and the experiment. In order to ensure that our conclusions are independent of the basis wave expansion, we have used different exchange–correlation functionals (LDA and GGA), and also used the ultrasoft pseudopotential plane

wave method [24] as implemented in QUANTUM-ESPRESSO package [25] for cross-checking of phonon calculations, but have not found any significant differences. The discrepancy of the calculated same phonon frequencies, including the B3g mode, are only within several cm−1 . We notice that, experimentally, temperature dependence of energy and linewidth of the B11g mode was suggested to come from a strong e–ph interaction for T < 40 K and ph–ph interaction (anharmonicity) for T > 40 K (Ref. [13]), while our above phonon calculations were within DFPT in the harmonic approximation, beyond some temperature, anharmonic effects might give a sizeable contribution. The large discrepancy is likely due to the possibility of vibrational anharmonicity. As a check on this assumption, first, we explore the possibility of anharmonicity of the B11g , B21g and B3g modes using the frozen phonon method. Then, the calculated e–ph linewidths and coupling strengths will be shown below. According to the eigendisplacements of the three phonon modes, Sb atoms were given eight different displacements u, with u/a, u/b, or u/c ranging from 0.005 to 0.040, and the eight calculated distorted energies Eu were fitted to a quartic polynomial Eu = k2 u2 + k4 u4 . The first-order term k2 u2 corresponds to the harmonic contribution, the fourth-order term k4 u4 relates to anharmonic effects, and the value of k4 /k22 represents the anharmonicity of phonon mode. Using Hartree–Fock decoupling, one gets Eu = (k2 + 3k4 ⟨u2 ⟩)u2 , where ⟨u2 ⟩ = h¯ /2M ωsch , h¯ = me = 1, MSb = 121.75 ∗ 1836 in atom unit, and ωsch = [(k2 + 3k4 ⟨u2 ⟩)/M ]1/2 is the self-consistent harmonic (sch) frequency [26,27]. The results are shown in Table 2. The fact that the values of k4 /k22 of B11g , B21g and B3g modes are all small [28] indicates that the phonons are all harmonic. As a result, the self-consistent harmonic frequencies of B11g , B21g and B3g modes are all equal to the DFPT calculations within the linear response method. Consequently, our DFT calculation of B3g mode of FeSb2 yields a relatively large difference with experimental data, and we hope that this result will motivate further studies to understand the discrepancy. The intrinsic linewidth γ in a perfect crystal is γ = γ e–ph + γ ph–ph , where γ e–ph and γ ph–ph represent the e–ph and ph–ph interactions, respectively [29,30]. Here, our calculated e–ph linewidths γ e–ph and coupling strengths λ at the Γ -point are also

234

R. Miao et al. / Solid State Communications 152 (2012) 231–234

Table 2 The fit coefficients of distorted energies Eu to a quartic polynomial according to the eigendisplacements of B11g , B21g and B3g modes, and the corresponding selfconsistent harmonic frequencies ωsch (in cm−1 ) using the frozen phonon method. Mode

k2 (Hartree/Bohr2 )

k4 (Hartree/Bohr4 )

ωsch (cm−1 )

B11g

0.163

−0.00581

187.1

B21g B3g

0.159 0.0430

−0.00729 8.27 ∗ 10−5

185.0 96.3

data, but the B3g mode shows a relatively large discrepancy. The calculated frozen phonon frequencies of B11g , B21g and B3g modes are all equal to the DFPT results within the linear response method, indicating that the phonons are all harmonic. Our calculated e–ph linewidth of the B11g mode is consistent with experimental data, and we also confirm the existence of e–ph interactions for the A1g , B11g , and B1u modes in FeSb2 . Acknowledgments

listed in Table 1. It can be seen from Table 1 that the e–ph linewidth and coupling strength of the B11g mode are both larger than that of all other modes. Considering the fact that phonon linewidth of the B11g mode drastically decreases with temperature decrease for

T < 40 K (see Fig. 2 of Ref. [13]), our calculated value of 0.41 cm−1 (T = 0 K) is consistent with experimental data 1.4 cm−1 (T = 15 K). Racu et al. suggested that the e–ph interaction is present in the B1u mode, while the temperature dependence of linewidths of A1g and the two B1g modes only come from ph–ph interactions, with no additional e–ph interactions [11]. On the contrary, Lazarević et al. argued that the temperature dependence of linewidth of the B11g mode comes from a strong e–ph interaction for T < 40 K [13]. We obtain λ ∼ 0.011 of the B1u mode, confirming the existence of e–ph interaction for the B1u mode [11]. The obtained λ ∼ 0.030 of the B11g mode is about three times that of the B1u mode, and larger than that of other modes. This result seems to support the existence of a strong e–ph interaction for the B11g mode of Ref. [13], or at least, the e–ph interaction is present in the B11g mode. The

calculated λ ∼ 0.015 of A1g mode is larger than that of B1u mode, indicating that an e–ph interaction is also present in A1g mode at low temperature. For the A2g mode, although the obtained λ ∼ 0.0093

of A2g mode is a little smaller than that of B1u mode, it seems to be non-negligible. It should to be pointed out that our electronic band structure calculations get a small hole pocket at R-point, finding a metal instead of a narrow-gap semiconductor for FeSb2 . This might influence the e–ph coupling as well as the phonon frequencies. It can be improved by performing a hybrid functional calculation, an LDA + U calculation or by applying the scissor operator. For example, Zhang et al. showed that inclusion of an effective Coulomb interaction within the LDA + U method can greatly influence the e–ph coupling strengths of phonon modes in cuprates [31]. Our calculated frequencies of B11g and B21g modes are quite close to each other, while the corresponding experimental data is different. It may be caused by the renormalization of the phonon frequency due to the presence of the strong e–ph interaction in B11g mode [32]. 4. Conclusions The lattice dynamics of FeSb2 have been computed using DFPT. The phonon frequencies together with e–ph linewidths at the Γ point are compared with the available experimental data. Our results support the A1g mode assignment of Refs. [11,12]. Nine phonon modes are in good agreement with the experimental

The author Rende Miao is very grateful to the Reviewer for reading the paper and helpful advice. This research was supported by the Pre-Research Foundation of PLA University of Science and Technology (Grant no. 20110520 and Grant no. 20110508). References [1] C. Petrovic, J.W. Kim, S.L. Bud’ko, A.I. Goldman, P.C. Canfield, W. Choe, G.J. Miller, Phys. Rev. B 67 (2003) 155205. [2] C. Petrovic, Y. Lee, T. Vogt, N.D. Lazarov, S.L. Bud’ko, P.C. Canfield, Phys. Rev. B 72 (2005) 045103. [3] R. Hu, V.F. Mitrovic, C. Petrovic, Phys. Rev. B 74 (2006) 195130. [4] R. Hu, V.F. Mitrovic, C. Petrovic, Appl. Phys. Lett. 92 (2008) 182108. [5] A. Bentien, S. Johnsen, G.K.H. Madsen, B.B. Iversen, F. Steglich, Europhys. Lett. 80 (2007) 17008. [6] A. Bentien, G.K.H. Madsen, S. Johnsen, B.B. Iversen, Phys. Rev. B 74 (2006) 205105. [7] R. Hu, V.F. Mitrović, C. Petrovic, Phys. Rev. B 79 (2009) 064510. [8] Y. Sun, S. Johnsen, P. Eklund, M. Sillassen, J. Bottiger, N. Oeschler, P. Sun, F. Steglich, B.B. Iversen, J. Appl. Phys. 106 (2009) 033710. [9] A. Lukoyanov, V. Mazurenko, V. Anisimov, M. Sigrist, T. Rice, Eur. Phys. J. B 53 (2006) 205. [10] H.D. Lutz, B. Müller, Phys. Chem. Miner. 18 (1991) 265. [11] A.-M. Racu, D. Menzel, J. Schoenes, M. Marutzky, S. Johnsen, B.B. Iversen, J. Appl. Phys. 103 (2008) 07C912. [12] N. Lazarević, Z.V. Popović, R. Hu, C. Petrovic, Phys. Rev. B 80 (2009) 014302. [13] N. Lazarević, Z.V. Popović, R. Hu, C. Petrovic, Phys. Rev. B 81 (2010) 144302. [14] J.M. Tomczak, K. Haule, T. Miyake, A. Georges, G. Kotliar, Phys. Rev. B 82 (2010) 085104. [15] W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) A1133. [16] X. Gonze, M.J. Verstraete, G. Zerah, J.W. Zwanziger, et al., Comput. Phys. Commun. 180 (2009) 2582. [17] N. Troullier, J.L. Martins, Phys. Rev. B 43 (1991) 1993. [18] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. [19] H.J. Monkhorst, J.D. Pack, Phys. Rev. B 13 (1976) 5188. [20] H.B. Schlegel, J. Comput. Chem. 3 (1982) 214. [21] S. Baroni, S. de Gironcoli, A. Dal Corso, P. Giannozzi, Rev. Modern Phys. 73 (2001) 515. [22] X. Wu, G. Steinle-Neumann, S. Qin, M. Kanzaki, L. Dubrovinsky, J. Phys.: Condens. Matter 21 (2009) 185403. [23] A. Perucchi, L. Degiorgi, R. Hu, C. Petrovic, V.F. Mitrović, Eur. Phys. J. B 54 (2006) 175. [24] D. Vanderbilt, Phys. Rev. B 41 (1990) 7892. [25] P. Giannozzi, et al., J. Phys.: Condens. Matter 21 (2009) 395502. [26] R.W.H. Stevenson (Ed.), Phonons in Perfect Lattices and in Lattices with Point Imperfections, Plenum Press, New York, 1966. [27] G.Q. Huang, Z.W. Xing, D.Y. Xing, Phys. Rev. B 82 (2010) 014511. [28] We changed the maximum values of u/a, u/b, or u/c both from 0.040 to a larger value 0.095 for B11g , B21g and B3g modes, respectively, but have not found significant differences in the results [29] A. Shukla, M. Calandra, M. d’Astuto, M. Lazzeri, F. Mauri, C. Bellin, M. Krisch, J. Karpinski, S.M. Kazakov, J. Jun, D. Daghero, K. Parlinski, Phys. Rev. Lett. 90 (2003) 095506. [30] A. Debernardi, S. Baroni, E. Molinari, Phys. Rev. Lett. 75 (1995) 1819. [31] P. Zhang, S.G. Louie, M.L. Cohen, Phys. Rev. Lett. 98 (2007) 067005. [32] P. Zhang, Y. Xue, P. Dev, Solid State Commun. 148 (2008) 151.