Interplay of Structural and Bonding Characters in Thermal Conductivity

Subscriber access provided by UNIVERSITY OF LEEDS

Article

Interplay of Structural and Bonding Characters in Thermal Conductivity and Born-Effective Charge of Transition Metal Dichalcogenides George Yumnam, Tribhuwan Pandey, and Abhishek K. Singh J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.7b11160 • Publication Date (Web): 15 Jan 2018 Downloaded from http://pubs.acs.org on January 15, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Interplay of Structural and Bonding Characters in Thermal Conductivity and Born-Effective Charge of Transition Metal Dichalcogenides George Yumnam, Tribhuwan Pandey, and Abhishek Kumar Singh∗ Materials Research Centre, Indian Institute of Science, Bangalore 560012, India E-mail: [email protected]

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Thermal transport in a material is governed by anharmonicity of crystal potential, which depends on the type of inter-atomic interaction. Using first-principles calculations, we report that lattice thermal conductivity (κlatt ) and its anisotropy (κx,y − κz ) of transition metal dichalcogenides (TMDs) increase by orders of magnitude with the change of constituent metal atom from Zr/Hf to Mo/W. This unprecedented difference in κlatt is substantiated by lower phonon group velocity, and four times larger anharmonicity of Zr/Hf based TMDs compared to Mo/W based TMDs. The sign and the absolute value of the Born effective charges, which emerges from the ionicity of the bonds, are found to be different for these two classes of materials. This leads to a significant difference in their interlayer van der Waals (vdW) interaction strengths, which are shown to be inversely related to the anisotropy in κlatt .

Introduction Understanding the origin, effects and implications of thermal transport in materials plays an essential role in advancing the forefronts of sustainable technological applications, such as waste heat recovery (thermoelectrics), 1 solar-energy harvesting (photovoltaics), 2 exploring novel properties (superconductivity), 3 minimizing transmission loss of electrical power, channelling excess heat in daily gadgets (cellphones and laptops), 4 etc. Thermal transport is controlled by the material’s thermal conductivity, which is an intrinsic property that originates from scattering of phonons with: (a) phonons, (b) electrons, (c) mass difference impurity (isotopes), or (d) discontinuity in the lattice (defects, grain boundaries, and interfaces). 5 The primarily dominating contributor of impedance to thermal transport is the phonon-phonon scattering, 5 which is governed by the strength of interactions arising from ionic, metallic or covalent nature of bonds. In most of the materials, atoms are bonded by one type of chemical interaction, giving rise to fairly isotropic transport. 6 However, materials having a mixture of bonds along different crystallographic directions, exhibit anisotropic transport. Layered materials offer such anisotropy in both thermal and electronic transport, due to their 2


Page 2 of 20

Page 3 of 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in-plane strong covalent or ionic bonds, and out-plane weak van der Waals (vdW) interactions. 7 Popular layered materials like graphene and boron-nitride are often not suitable for semiconductor devices due to incompatible bandgaps. Meanwhile, transition metal dichalcogenides (TMDs) have been considered promising for applications in thermal, electronic, and optoelectronic devices due to their moderate bandgaps. 7–10 They show tunability of electronic and thermal transport properties via external perturbations such as application of strain, pressure, and doping. 8,9,11 Due to the layered structure of TMDs, the metal-chalcogen (M-X) bonds along in-plane, x-y, direction is stronger than the interlayer vdW interaction along cross-plane, z, direction, giving rise to anisotropy in lattice dynamics. Using first principles calculations in combination with lattice Boltzmann transport theory, we elucidate the origin of unprecedented variation of κlatt and its anisotropy in Gr-IV (Zr and Hf) and Gr-VI (Mo and W) TMDs. The variation of κlatt is maximum while going from ZrSe2 to WS2 , with a change of κx,y (κz ) from ∼5 (∼2.5) W/m-K to ∼290 (∼40) W/m-K, respectively. The lower κlatt in Gr-IV TMDs is substantiated by four times larger mode Grüneisen parameter (γ) and twice as large anharmonic scattering rate (Wanhar ) as compared to Gr-VI TMDs. The Born effective charges are found to be very different for these two classes of materials, which lead to a significant difference in their interlayer van der Waals (vdW) interaction strengths. We demonstrate an inverse correlation between the strength of vdW interaction and anisotropy in thermal conductivity of these TMDs. Our findings can be in principle generalized to other class of materials having mixedbonding characters.

Computational Details The calculations were performed using density functional theory (DFT) as implemented in Vienna ab initio simulation package (VASP). 12–15 The interactions of electrons and ions were treated using all-electron projector augmented-wave (PAW) pseudopotentials. 13 The Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional was used within the generalized gradient approximation

3



Page 4 of 20

(GGA), 15 with a wave function cut-off energy of 500 eV, and energy convergence criterion of 10−8 eV. The effects of dispersion interactions were incorporated by using optimized vdW correlation functionals: optB86b-vdW, and optB88-vdW 16,17 for Gr-IV and Gr-VI TMDs, respectively (See Figure S-1 of Supporting Information). The Born effective charge is calculated by taking account of long range electrostatic interactions within density functional perturbation theory. 18 κlatt is obtained by calculating the iterative solution to phonon Boltzmann transport equation (PBTE), 19 by considering the three phonon scattering processes within the single mode relaxation time approximation (SMRTA). 20 Within the framework of PBTE, κlatt at a temperature T can be expressed as

κlatt = καα =

1 X ∂f (~ωλ )vλα vλα τλ N V λ ∂T

(1)

where N , V , f and vλα are the number of uniformly spaced q points in the Brillouin zone, volume of the unit cell, the Bose-Einstein distribution function, which depends on the phonon frequency ωλ and the velocity along α direction. The numerical calculation of κlatt via this method requires the harmonic second-order interatomic force constants (IFCs) along with anharmonic third-order IFCs. The second (third) order IFCs were calculated using Hellmann-Feynman theorem, as implemented in VASP 12–15 within the finite-displacement method on a supercell of 5×5×3 (3×3×2), as implemented in Phonopy 21 (thirdorder.py 22,23 ). A converged k-mesh of 3×3×2 was used for the Brillouin zone integration. Finally, the linearized PBTE is numerically solved by using the ShengBTE code, 22–24 which has been shown to accurately predict κlatt of different materials. 25–27 For performing the pCOHP calculations, the required plane-wave calculations were performed using PBE-GGA within PAW as implemented in VASP with a fine-meshed Γ-centered k-point of 16 × 16 × 16 sampled with the tetrahedron method with Blöchl corrections. The pCOHP calculations were performed with the LOBSTER code. 28

4


Page 5 of 20

(a)

X M X

(b)

X M X

Frequency (cm-1)

400

Frequency (cm-1)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


AP

q-AOP

OP

PDOS

(c) ZrS2

Zr S

300

400

200

100

100 0 Γ M K ΓA L H Zr 300 (f) WSe2

LH A0 1 2 3

Se

200

200

100

100

0 Γ MK

(d) WS2

300

200

0 Γ M K ΓA 300 (e) ZrSe2

PDOS

ΓA L H A0 1 2 3 0Γ M K

ΓA L H

W S

A01 2 3

W Se

A0 1 2 3

Figure 1: Crystal structures of (a) Gr-IV TMDs (Zr/Hf)X2 with space-group P¯3m1 and (b) Gr-VI TMDs (Mo/W)X2 with space-group P63 /mmc. The phonon dispersion and corresponding phonon density of states (PDOS) of (c) ZrS2 , (d) WS2 , (e) ZrSe2 , and (f) WSe2 .

Results and Discussion In our investigation, we have chosen ZrSe2 and WS2 for representing Gr-IV and Gr-VI TMDs as they portray the opposite extremes of anisotropy in κlatt . Figure 1 (a) and (b) show the structures of Gr-IV and Gr-VI TMDs having space-groups P¯3m1 and P63 /mmc, respectively. The adjacent layers of Gr-IV and Gr-VI TMDs are A-A and A-B stacked, respectively. As shown in Figure 1 (c - f), there are no imaginary modes in the phonon dispersion for all computed structures which indicate their dynamical stability. ZrX2 has 3 acoustic phonon (AP) branches, whereas WX2 has 3 low-lying quasi-acoustic optical phonon (q-AOP) branches in addition to 3 AP branches. These low lying q-AOP branches in WX2 arise due to in-phase rigid-layer shear and vertical motion of the atoms. 29 The phonon density of states (PDOS) show that Gr-IV TMDs have no phonon band gap, whereas a finite gap exists between optical (OP) and AP/q-AOP branches in Gr-VI TMDs. Gr-IV and Gr-VI TMDs have respectively 6 and 12 normal modes of lattice vibration in optical 5



Page 6 of 20

branch at zone center, which can be described by the irreducible representations, Eqn. (2) and Eqn. (3), respectively: 29–31 ΓGr-IV =A1g ⊕ E1g ⊕ 2A2u ⊕ 2E1u

(2)

ΓGr-VI =A1g ⊕ 2A2u ⊕ B1u ⊕ 2B2g ⊕ E1g ⊕ 2E1u ⊕ E2u ⊕ 2E2g

(3)

Figure 2: The Grüneisen parameter (γ) of (a) ZrS2 , (b) ZrSe2 , (c) WS2 , and (d) WSe2 . The anharmonic scattering rate (Wanhar ) of (e) ZrS2 , (f) ZrSe2 , (g) WS2 , and (h) WSe2 . 1 Gr-IV TMDs have two Raman active (A1g , E1g ) and two infra-red (IR) active (A12u , E1u ) modes, 1 2 1 whereas Gr-VI TMDs have four Raman active (A1g , E1g , E2g , E2g ), two IR (A12u , E1u ), and four 1 2 1 inactive (B1u , B2g , B2g , E2u ) modes. As shown in Figure 1 (c - f), the frequencies of branches

representing similar modes of vibration decrease with increasing mass of chalcogen in both class of TMDs. Moreover, the frequency of TA1 branch sensitively decreases to the increase in mass of transition metal atom with the same chalcogen atom. Consequently, the anharmonicity of crystals shows a striking difference with respect to change in transition metal atom, as shown by a huge difference in mode Grüneisen parameter γ, which is defined as the negative logarithmic derivative of frequency, ωq with respect to volume (V) for qth mode of vibration: γ = −∂(ln ωq )/∂(ln V ). 32 γ of Gr-IV TMDs are nearly four times higher than Gr-VI TMDs as shown in Figure 2 (a-d). The larger anharmonicity in Gr-IV TMDs is substantiated by twice as large anharmonic scattering rate (Wanhar ) which describes higher three phonon scattering phase space volume in Gr-IV TMDs as 6


Page 7 of 20

compared to Gr-VI, as shown in Figure 2 (e-h). Larger anharmonicity signifies a larger non-linear restoring force of vibration leading to higher scattering of phonons, 32 which is generally observed in materials with weak interactions between structural components, such as van der Waals (vdW) materials. 33 Previous study has shown that vdW materials exhibit anisotropic thermal 34 transport, though the nature and origin of such anisotropy in κlatt is not fully explored.

15

10

(a)

200 ZrS2 (x,y) ZrS2 (z) 150 ZrSe2 (x,y) ZrSe2 (z) 100

(b)

WS2 (x,y) WS2 (z) WSe2 (x,y) WSe2 (z)

103

ζ = |κxy-κz| (W/m.K)

20

κii (W/m.K)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


102

(c)

ZrS2 ZrSe2

WS2 WSe2

101

5

50

0

0

400 600 800 1000 1200 1400

T (K)

100

400 600 800 1000 1200 1400

T (K)

400 600 800 1000 1200 1400

T (K)

Figure 3: The directional component of κlatt along in-plane (x, y) and cross-plane (z) directions of z (a) ZrS2 and ZrSe2 , (b) WS2 and WSe2 , and (c) the absolute anisotropy of κlatt (ζ = κx,y latt − κlatt ) of ZrS2 , ZrSe2 , WS2 , and WSe2 . As shown in Figure 3 (a, b), Gr-VI exhibits orders of magnitude larger κlatt along both inplane (x, y) and cross-plane (z) directions as compared to Gr-IV TMDs. For instance, κlatt of WS2 along in-plane (κx,y ) and cross-plane (κz ) directions are ∼60 and ∼20 times larger than that of ZrSe2 , respectively. Moreover, the absolute anisotropy defined as ζ = κx,y − κz of WS2 is 100 times larger than ZrSe2 , as shown in Figure 3 (c). It can be noted that the difference in crystal structure of these TMDs might also affect the κlatt . To investigate the effect of different crystal structures, we considered both crystal structures (P¯3m1, and P63 /mmc) for a given TMD. However, Gr-IV and Gr-VI TMDs exhibit imaginary modes showing dynamical instability in P63 /mmc and P¯3m1 structure, respectively; hence, completely ruling out the calculation of thermal conductivity in these configurations. Though, this limits the first-principle calculation of κlatt , previous multiscale modelling study 35 have shown that the room temperature thermal conductivity of MoS2 can change from 40±4 W/m-K to 32±3 W/m-K when its crystal structure changes from 2H (P63 /mmc) to 1T (P¯3m1) phase, respectively which is relatively very small. Further details of these differences in κlatt and its anisotropy are obtained from the mode resolved κilatt , where i = TA1, TA2, LA, 7



2 1.5 1 0.5 0 2 1.5 1 0.5 0 2 1.5 1 0.5 0 2 1.5 1 0.5 0

κlatt (W/m.K)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 20

ZrSe2

WS

T (K)

T (K)

2 40 κx = κy 30 (a) TA1 κz 20 10 0 40 (b) TA2 30 20 10 0 40 (c) LA 30 20 10 0 150 (d) OP 100 50 0 400 600 800 100012001400 400 600 800 100012001400

Figure 4: The contribution of (a) TA1, (b) TA2, (c) LA, and (d) OP branch (phonon modes) to the directional lattice thermal conductivity (κx,y , κz ) of ZrS2 and WS2 , respectively. and OP, which are the two transverse acoustic, longitudinal acoustic, and optical phonon modes, respectively. As shown in Figure 4, κz > κx,y in the TA1 mode of all Gr-IV TMDs, which minimizes the anisotropy arising from κz < κx,y in other modes. However, κx,y > κz for all modes in Gr-VI TMDs, leading to a larger anisotropy as compared to Gr-IV TMDs. Furthermore, the contribution of OP and q-AOP modes to κlatt are very large (60%) in WS2 , as compared to just 15% in ZrSe2 . 36 The coupling of AP and q-AOP modes in WS2 leads to a much larger contribution of optical modes (q-AOP + OP) to κlatt , as compared to ZrSe2 . Another major contribution to κlatt comes form phonon group velocity (vg ) as κlatt ∝ vg2 . As shown in Figure 5, the cross-plane vgz of ZrSe2 and all other Gr-IV TMDs are higher than in-plane vgx,y for TA1 mode, corroborating the trend in κT A1 . However, vgx,y of WS2 and other Gr-VI TMDs are greater than vgz for all AP modes. Anharmonicity and phonon group velocity in a crystal is controlled by the nature of interaction between constituent atoms, therefore, we investigated the electron localization function (ELF), which maps the electron pair probability in multi-electron systems. 38 A qualitative impression of

8


Page 9 of 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5: The directional component of group velocity (vxx , vyy , vzz ) for (a) TA1, (b) TA2, and (c) LA phonon modes of ZrSe2 and WS2 . Inset in (a) shows that vzz > vx,y for TA1 mode in ZrSe2 . the ionicity or covalency of bonds can be obtained from the form of ELF. 39 As shown in Figure 6 (a), the ELF of ZrSe2 shows perfect localization (ELF = 1.0) around Se atom, and electron depletion (ELF < 0.1) around Zr atom. This is a characteristic of ionic bond where electrons are spatially localized around the higher electronegative element. 32,37 However, as shown in Figure 6 (b), in the case of WS2 electrons are shared between W and S atoms with a gradient of ELF from W = 0.5 to S = 0.9, 37 which indicates polar covalent character of W-S bond. The difference in nature of bonds between these two compounds lies in the electronegativity (X el ) difference between constituent metal and chalcogen atom which controls the ionic character of a bond. For instance, 9



1.0

3

(a) ZrSe2 Zr

(c)

M

2

Se

M+

X--

Zzz*

1 ZrS2 ZrSe2 HfS2 HfSe2

0

X

0.1

0

*

-1 -4 -2

2

6

4

*

8

Zxx ≈ Zyy

(b) WS2 1.0

0.5

S W

Zzz*

0

(d)

MoS2 MoSe2 WS2 WSe2

X

M

-0.5 M--

0.1

-1.5 -1 -0.5

X+

0

*

*

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 20

0.5

Zxx ≈ Zyy Figure 6: The electron localization function (ELF) of (a) ZrSe2 and (b) WS2 along (110) plane. In the case of ZrSe2 , electrons are perfectly localized around the Se atoms (ELF=1.0) and surrounded by the electron gas (ELF=0.5), whereas Zr atoms are completly delocalized (ELF

Interplay of Structural and Bonding Characters in Thermal Conductivity

Recommend Documents