Anharmonic inter-layer bonding leads to intrinsically low thermal conductivity of bismuth oxychalcogenides

TL;DR

This study reveals that electrostatic inter-layer coupling and strong anharmonic chemical bonds in bismuth oxychalcogenides cause intrinsically low thermal conductivity, providing insights into phonon anharmonicity linked to chemical bonding.

Contribution

The paper establishes a direct connection between inter-layer chemical bonding and phonon anharmonicity using density functional theory calculations.

Findings

Electrostatic inter-layer coupling reduces thermal conductivity.

Strong anharmonic Bi-chalcogen bonds are responsible for phonon scattering.

Low thermal conductivity is intrinsic to Bi₂O₂X materials.

Abstract

The anharmonicity of phonons in solid is ultimately rooted in the chemical bonding. However, the direct connection between phonon anharmoncity and chemical bonding is difficult to make experimentally or theoretically, due mainly to their complicated lattice structures. Here, with the help of density functional theory based calculations, we discovery that electrostatic inter-layer coupling in Bi$_2$O$_2$X (X=S,Se,Te) leads to intrinsically low lattice thermal conductivity. We explain our discovery by the strong anharmonic chemical bonding between Bi and chalcogen atoms. Our results shed light on the connection between inter-layer chemical bonding and phonon anharmonicity, which could be explored in a wide range of layered materials.

Tables

Table 1. Table 1: Projected Grüneisen parameters γ 𝛾 \gamma on different atoms and directions in Bi 2 O 2 X.

		Bi(x)	Bi(y)	Bi(z)	X(x)	X(y)	X(z)	O(x)	O(y)	O(z)
Bi2O2S		4.79	4.79	4.02	4.70	2.40	4.93	1.30	1.43	1.48
Bi2O2Se		4.27	4.27	3.42	4.05	4.05	4.34	1.96	1.96	1.93
Bi2O2Te		3.60	3.60	3.07	3.76	3.76	5.20	1.71	1.71	1.50

Table 2. Table 2: Calculated lattice parameters, static dielectric tensor ( ε 𝜀 \varepsilon ) of BOX.

	$a$ (Å)	$b$ (Å)	$c$ (Å)	$\epsilon$x/y	$\epsilon$z
Bi2O2S	3.98	3.89	12.08	10.02/11.36	9.42
Bi2O2Se	3.93	3.93	12.40	13.86/13.86	10.34
Bi2O2Te	4.02	4.02	12.88	17.21/17.21	12.00

Table 3. Table 3: Calculated Born effective charge ( z 𝑧 z ) for ions in the unit of e 𝑒 e for BOX.

	$z_{\mathrm{Bi},x/y}$	$z_{\mathrm{Bi},z}$	$z_{X,x/y}$	$z_{X,z}$	$z_{O,x/y}$	$z_{O,z}$
Bi2O2S	5.17/5.91	5.17	-3.30/-4.32	-3.26	-3.50/-3.75	-3.54
Bi2O2Se	6.12/ 6.12	5.39	-4.33/-4.33	-3.12	-3.96/-3.96	-3.82
Bi2O2Te	6.54/6.54	5.67	-4.35 /-4.35	-3.21	-4.33/-4.33	-4.05

Table 4. Table 4: Projected Grüneisen parameters γ 𝛾 \gamma on different atoms and directions in BOX and BiCuOX.

		Bi(x/y)	Bi(z)	X(x/y)	X(z)	O(x/y)	O(z)	Cu(x/y)	Cu(z)
Bi2O2S		4.79/4.79	4.02	4.70/2.40	4.93	1.30/1.43	1.48	-	-
Bi2O2Se		4.27/4.27	3.42	4.05/4.05	4.34	1.96/1.96	1.93	-	-
Bi2O2Te		3.60/3.60	3.07	3.76/3.76	5.20	1.71/1.71	1.50	-	-
BiCuOS		2.80/2.80	2.37	3.16/3.16	2.53	1.99/1.99	2.0	2.74/2.74	3.50
BiCuOSe		3.50/3.50	3.00	3.30/3.30	2.81	1.98/1.98	1.93	3.95/3.95	3.26
BiCuOTe		3.47/3.47	2.90	3.05/3.05	2.97	2.26/2.26	2.06	3.56/3.56	2.50

Equations

$${\vec{\kappa }}_{\lambda }=C_{\lambda }{\vec{v}}_{\lambda }\otimes {\vec{v}}_{\lambda }{\tau }_{\lambda },$$

$$\vec{\kappa }(\omega )=\sum _{\lambda }\delta (\omega -{\omega }_{\lambda }){\vec{\kappa }}_{\lambda }.$$

$${\vec{\kappa }}_s(\omega )=\sum _{{\omega }_{\lambda }\le \omega }{\vec{\kappa }}_{\lambda },$$

$$\vec{\kappa }=\frac{1}{V}\sum _{\lambda }{\kappa }_{\lambda },$$

$$D(\omega )$$

$${\gamma }_{i\alpha }=\frac{{\sum }_{\lambda }{\gamma }_{\lambda }{e}_{\lambda ,i\alpha }}{{\sum }_{\lambda }{e}_{\lambda ,i\alpha }}.$$

Full text

I Introduction

Materials with low thermal conductivity may find applications in many disciplines including thermoelectricsHe and Tritt (2017); Zeier et al. (2016), heat insulation and phononic devicesLi et al. (2012). Phonons are main heat carriers in semiconductors and insulators. Phonon engineering of thermal conductivity has witnessed tremendous progress in recent yearsBera et al. (2010); Vineis et al. (2010); Luckyanova et al. (2012). Different extrinsic and intrinsic approaches have been developed to reduce the phonon thermal conductivity ($\kappa$). Among them are nanostructuring, defect engineering, and enhancing phonon anharmonicity. Anharmonic phonon scattering is an intrinsic mechanism that leads to finite $\kappa$He and Tritt (2017); Li and Mingo (2015); Carrete et al. (2014a); Mukhopadhyay et al. (2018). Thus, utilizing strong phonon anharmonicity is an attractive way to reduce $\kappa$, while keeping other properties intactJana and Biswas (2018). Generally, weak chemical bonds lead to large anharmonicity due to large atomic displacement involved. As a rule of thumb, complicate lattice structure with heavy elements is believed to lead to strong anharmonicity. But their microscopic mechanism is difficult to pinpoint.

Layered materials have weak inter-layer bonding and many show low $\kappa$ and good thermoelectric performanceZhou and Zhao (2017). They normally have good electrical transport property in the plane and low $\kappa$ across the plane. For thermoelectric applications, it is highly desirable to reduce the in-plane $\kappa$ ($\kappa_{\parallel}$) through phonon engineering. As a kind of thermoelectric material, layered bismuth oxychalcogenides Bi2O2X (X=S, Se, Te) (BOX) have been studied experimentallyRuleova et al. (2010); Luu and Vaqueiro (2015); Tan et al. (2018); Luu and Vaqueiro (2016); Zhang et al. (2013). But their thermoelectric performance is mainly hindered by the poor electrical transport property. Very recently, thin layer of single crystal layered oxychalcogenide, Bi2O2Se, has been successfully synthesizedWu et al. (2017a, b). Its high electron mobility, strong spin-orbit interaction and ultrafast infrared response lead to potential applications in nano-electronics, opto-electronics, topological devices and ferroelectricityYin et al. (2018); Meng et al. (2018); Wu and Zeng (2017). The electronic band structure of Bi2O2Se has been mapped out combining angle-resolved photoemission spectroscopy and density functional theory (DFT) calculationsChen et al. (2018). The experimental progress makes it possible to enhance the thermoelectric performance by, i.e., tuning of carrier concentrationWu et al. (2018).

Despite the above mentioned progress in characterizing the electric properties, the phonon transport properties of BOX are still poorly understood, which hinders the in-depth understanding in its electronic, optoelectronic and thermal properties. Here, using DFT based calculations, we show that, BOX has low intrinsic in-plane $\kappa$ ($\kappa_{\parallel}$). This means high electrical conductivity and low phonon thermal conductivity can in principle be realised simultaneously in the same direction. Based on the analysis of its phonon spectrum, scattering lifetime, Grüneisen parameters and real-space electron distribution, we are able to show unambiguously that, the low $\kappa_{\parallel}$ originates from strong anharmonic inter-layer bonding between Bi and chalcogen atoms. Our results shed light on the connection between chemical bonding and phonon anharmonicity, and showed that inter-layer coupling can be used to tune $\kappa_{\parallel}$ in layered materials.

II Results and discussions

II.1 Structure and phonon dispersion

Bulk Bi2O2Se and Bi2O2Te crystallize in the (Na0.25Bi0.75)2O2Cl type structure, and belong to tetragonal space group I4/mmm (139) with 10 atoms in one unit cell as shown in Figure 1 (b). Meanwhile, Bi2O2S has a distorted structure, where the Bi atoms slide slightly apart, and belongs to the Pnnm (58) group with lower symmetry (Figure 1 (a)). This structure distortion changes the chemical bonding environment of S and brings a small anisotropy between $x$ and $y$ direction. We will show below that it has important effect on $\kappa$. There are two electron transfer from the Bi2O2 layer to chalcogen layer, and the two layers are bonded through electrostatic force. The optimized lattice parameters, the dielectric constant and Born effective charge of each atom are listed in Table 1 and 2 of the Supporting Information (SI). We have also calculated the electronic band structures from the optimized lattice parameters. They show good agreement with previous works (SI, Figure 1).

Figure 1 (d-f) show the atomically resolved phonon dispersion and the corresponding density of states (DOS) projected onto different atoms. We observe a phonon gap at $\sim 6$ THz, separating the spectrum, albeit small in Bi2O2S. The high frequency part of the spectrum is contributed dominantly by O atoms. The corresponding dispersion and DOS are very similar for all three materials. Although O and Bi atoms are strongly bonded in the same layer, their motions are decoupled due to large mass mis-match. On the other hand, although Bi and chalcogen atoms are at different layers, they couple together and form the low frequency phonon band. This clear separation of O vibration from others makes our analysis easier.

II.2 Thermal conductivity

We calculate $\kappa$ using the Boltzmann transport equation (BTE) within the single mode relaxation time approximation (RTA) as implemented in PhonopyTogo and Tanaka (2015) and Phono3pyTogo et al. (2015) package. The results are shown in Figure 2. More details of the results are shown in Figures 2-4 of SI. Consistent with their layered structure, all materials show anisotropic $\kappa$, with $\kappa_{\parallel}>\kappa_{z}$. Bi2O2S and Bi2O2Te show similar low $\kappa_{x/y}$, which is below 1 W/m-K at T=300 K, while Bi2O2Se has around two-fold larger value. We note that, $\kappa$ we obtained for BOX is comparable to values of other thermoelectric materials, including SnSeCarrete et al. (2014b); Xiao et al. (2016); Li et al. (2015), Bi2Te3, PbX (X=S, Se, Te) and BiCuOX (X=S, Se, Te)Lee et al. (2014); Ji et al. (2016); Shao et al. (2016). Due to their similar lattice structure, we have performed detailed comparison to BiCuOX in the SI (Figures 5-9 and Table 3). Additionally, we notice that Bi2O2S shows lower $\kappa$ than Bi2O2Se. Normally, from S to Te, we expect $\kappa$ to decrease with increasing atomic mass for similar lattice structures. We may attribute this abnormal behavior of Bi2O2S to its structure distortion. The purpose of the rest analysis is to give a microscopic explanation of the above mentioned features in $\kappa$.

As shown in Eq. (2) of Methods, within all the parameters that $\kappa$ depends, the relaxation time $\tau_{\lambda}$ is determined by anharmonic phonon scattering, while all the rest terms are determined by the harmonic phonon spectrum. To find the origin of low $\kappa$ in Bi2O2X, firstly we set $\tau_{\lambda}=1$, and show the accumulative sum of frequency dependent $\kappa(\omega)$ (Eq. 3 of Methods) in Figure 3 (a). Bi2O2S shows the largest value, opposite to the results of $\kappa$. This means lower $\kappa$ of Bi2O2S does not come from the harmonic phonon spectrum, but from their shorter $\tau_{\lambda}$. We have also shown a comparison to BiCuOX (X=S, Se, Te) in Figure 7 of the SI. We find that, although they have similar $\kappa$, the band structure contribution to $\kappa$ in BOX is much larger than BiCuOX. This suggests stronger anharmonic scattering and hence smaller $\tau_{\lambda}$ in BOX, which we focus on in the following.

II.3 Anharmonic scattering

Figure 4 shows the distribution of scattering lifetime $\tau$ as a function of frequency for the three materials considered. The high frequency part shows little frequency dependence and is of similar magnitude for all materials. The low frequency part has a much wider distribution. But it is still clear that, the data of Bi2O2Se locate in higher $\tau$ region compared to Bi2O2Te and Bi2O2S. This confirms the expectation that, lower $\kappa$ of Bi2O2S and Bi2O2Te originates from its shorter $\tau$ and hence stronger anharmonic scattering. Considering their similar lattice structure, shorter $\tau$ of Bi2O2Te compared to Bi2O2Se can be understood due to larger mass of Te. However, this can not explain why Bi2O2S shows the lowest $\tau$ distribution in the low frequency regime. Thus, we need further analysis of $\tau$.

The magnitude of $\tau$ depends both on the scattering phase space and on the strength of anharmonic potential. We can separate the two contributions. Figure 3(b) compares the joint DOS (JDOS) of 3-phonon scattering in all materials. We find that with increasing atomic number from S, Se to Te, JDOS in the low frequency regime grows up. The reason is that, heavier mass leads to narrower phonon spectrum and stronger overlap between different modes. This tends to increase the phase space for anharmonic phonon scattering. Although the change of JDOS may partly explain different $\kappa$ of Bi2O2Se and Bi2O2Te, it still can not explain the lower $\kappa$ of Bi2O2S.

We now turn to the anharmonic potential, the strength of which can be characterized by the Grüneisen parameters ($\gamma$). To further correlate with the atomic chemical bonding, in Table 4, we show the projected $\gamma$ on different atoms in three directions. As a common feature, we get larger $\gamma$ for Bi and chalcogen atoms compared to O atoms. These results suggest that inter-layer Bi-X bonding in BOX is strongly anharmonic and generates stronger phonon scattering. This is the common feature of all three kinds of materials. The low $\kappa$ of Bi2O2S can also be understood from the Grüneisen parameter. The lattice distortion in Bi2O2S makes the $x$ and $y$ direction anisotropic. This is reflected in the Grüneisen parameters. $\gamma$ of S in $y$ direction is reduced, while that in $x$ and $z$ gets much larger, promoting the inter-layer anharmonicity. This information suggests that, we can attribute the smaller scattering lifetime and lower $\kappa$ of Bi2O2S to the enhancement of anharmonic inter-layer coupling generated by lattice distortion.

Further comparison to BiCuOX supports our above arguments. In BiCuOX, chalcogen atom and Cu form stronger bonds. This reduces the anharmonic inter-layer coupling between BiO and CuX layers. This is reflected in the projected Grüneisen parameters. First, in BiCuOX, $\gamma$ in $z$ direction is consistently smaller than that of in-plane ($x/y$) direction for all atoms. Second, Cu atoms show the largest $\gamma$ instead of chalcogen atoms in BOX. The reduction of $\gamma$, together with their smaller JDOS (SI, Figure 8), leads to much longer scattering lifetime (SI, Figure 9).

II.4 Origin of the strong anharmonicity

From the above analysis, we can attribute low $\kappa$ of BOX to large anharmonic coupling between Bi and chalcogen atoms. The lattice distortion of Bi2O2S further enhances the anharmonicity. The next question we ask is what is the microscopic mechanism of the strong anharmonicity. We will now make connection to its atomic structure and chemical bonding. In BOX, each chalcogen atom is surrounded by a Bi cage with eight Bi atoms. Due to the inter-layer charge transfer, the chalcogen atom interacts mainly through electrostatic force with surrounding Bi atoms, forming weak ionic bonds. Their bonding environment can be characterized through the electron localization function (ELF)Becke and Edgecombe (1990) (Figure 5), which shows clear ionic bonding between atoms. The small caps near Bi represents its lone-pair electrons (LPE), which also contribute to the inter-layer coupling, i.e., they interact through electrostatic force with the chalcogen atoms. As for Bi2O2S, after lattice distortion, four Bi atoms move even closer S. This further enhances their mutual interaction, which can be seen from the plane cut in the lower panel of Figure 5 (a). Comparing to (b) and (c), the strongest bonding between S and Bi can be deduced. This is the microscopic reason why Bi2O2S has $\kappa$ as low as Bi2O2Te and much lower than Bi2O2Se.

III Conclusions

To conclude, we have predicted low phonon thermal conductivity of bismuth oxychalcogenides Bi2O2X (X=S, Se, Te). Through careful analysis of their phonon properties, we can ultimately correlate the strong anharmonicity and hence low thermal conductivity to the inter-layer bonding between Bi and chalcogen atoms. The strong correlation between bond anharmonicity and low thermal conductivity gives atomic insights of the thermal properties of materials. The same principle can be applied to a broad range of layered materials with electrostatic inter-layer coupling.

IV Methods

IV.1 DFT calculations

For the first-principles calculations, we use density functional theory (DFT) with the projected augmented wave (PAW) method as implemented within the Vienna $ab$ $initio$ Simulation Package (VASP)Kresse and Furthmüller (1996a, b). We choose the Perdew-Bueke-Ernzerhof (PBE)Perdew et al. (1997) version of generalized gradient approximation (GGA) to treat the exchange-correlation interaction. The plane wave cut-off energy is set as $550$ eV. The Brillouin zone is sampled by using the Monkhorst-Pack schemeMonkhorst and Pack (1976) with $9\times 9\times 3$ mesh k-points to optimize the structure until the forces on the atoms are less than $0.01$ eV/Å.

To calculate the phonon dispersion and phonon conductivity, we use phonopyTogo and Tanaka (2015) and phono3pyTogo et al. (2015) codes together with VASPKresse and Furthmüller (1996a, b). The second and third order force constants are calculated by finite-difference method. We use $5\times 5\times 2$ supercell for the second order force constant, and $3\times 3\times 2$ supercell for the third order force constant. We have also performed the $4\times 4\times 2$ supercell calculation to confirm the convergence of the third order forces. We set the convergence criteria to $10^{-8}$ eV for self-consistent loop and $0.01$ eV/Å for the force. The $\Gamma$-only scheme is used to sample the reciprocal cell of the supercell.

IV.2 Thermal conductivity

The linearized phonon BTE with the single mode RTA is used to calculate the phonon thermal conductivity with $21\times 21\times 21$ sampling mesh using phono3pyTogo et al. (2015). We first define

[TABLE]

as contribution of each mode $\lambda$ to the thermal conductivity. Here, $C_{\lambda}$ is the heat capacity of mode $\lambda$, $\vec{v}_{\lambda}$ is the group velocity, and $\tau_{\lambda}$ is the relaxation time. The frequency-resolve version is then written as

[TABLE]

The accumulated sum plotted in Fig. 4 is calculated as

[TABLE]

with $\tau_{\lambda}$ set to 1. Finally, the thermal conductivity $\kappa$ is expressed as

[TABLE]

where $V$ is the volume of the system.

IV.3 The joint density of

states

The joint density states can be used to quantify the phase space for phonon anharmonic scattering, determined by the phonon dispersion relation. The results shown in Fig. 4 (b) is for 3-phonon scattering and calculated from

[TABLE]

IV.4 The projected Grüneisen parameter

The mode-resolved Grüneisen parameters are calculated within the quasi-harmonic approximation using Phonopy. To characterize the anharmonicity of each atom, we define the projected Grüneisen parameter by projecting all the modes to given atom in a given direction as following

[TABLE]

Here, $\gamma_{\lambda}$ is the Grüneisen parameter of mode $\lambda$, $i$ is the atom index, $\alpha$ represents the direction, and $e_{\lambda,i\alpha}$ is the element of the eigen vector $e_{\lambda}$ corresponding to atom $i$ in direction $\alpha$.

Acknowledgements.

The authors are supported by the National Key Research and Development Program of China (Grant No. 2017YFA0403501), the National Natural Science Foundation of China (Grant No. 21873033) and the program for HUST academic frontier youth team. They thank the National Supercomputing Center in Shanghai for providing computational resources.

Appendix A Electronic structure

Based on the relaxed lattice structure, the electronic band structure is calculated in the modified Becke-Johnson meta-GGA potentialBecke and Johnson (2006); Tran and Blaha (2009) including the spin-orbit interaction. The results are shown in Figure 6. We get a band gap of 1.38 eV, 0.91 eV and 0.23 eV for Bi2O2S, Bi2O2Se and Bi2O2Te, respectively. They are comparable to experimental observations. The DOS plots show that at the top of the valence band, the contribution of S/Se/Te atom is dominant, while at the bottom of the conduction band the contribution of Bi atoms is dominant. Further analysis shows that, for all cases, the contribution mainly comes from the $p$ orbitals.

Appendix B Phonon properties of Bi2O2X

Tables 2-3 give the calculated lattice parameters, the static dielectric constants and the Born effective charges of Bi2O2X (X=S, Se, Te) (BOX). Figures 7-9 give the details of the phonon thermal conductivity ($\kappa$) calculation.

Appendix C Comparison to BiCuOX

BiCuOX (X=S, Se, Te) has similar structure to BOX. Figure 10 shows their lattice structure, corresponding phonon spectrum, and projected density of states onto different atoms. $\kappa$ of BiCuOX is also comparable to BOX. Thus, we here give some comparison between these two kinds of materials. Figure 11 shows that, $\kappa$ of BOX and BiCuOX is of similar magnitude. Figure 12 shows that, the phonon band structure contribution to $\kappa$ is different for these two kinds of materials. BOX would have larger $\kappa$ than BiCuOX if they have the same relaxation time. Thus, the scattering lifetime $\tau_{\lambda}$ in BOX is relatively smaller than that in BiCuOX.

There are two possible factors that influence $\tau_{\lambda}$. First, Figure 13 shows that, the joint density of states (JDOS) for 3-phonon scattering in BOX are larger than that in BiCuOX. Larger JDOS gives rise to smaller $\tau_{\lambda}$. Second, comparison of projected Grüneisen parameters (Table 4) shows that, the anharmonicity in BOX is larger than that in BiCuOX. Especially, the anharmonicity of inter-layer bonding is stronger in BOX. This can be seen from the relative magnitude of in-plane ($x/y$) and out-of-plane ($z$) $\gamma$ of chalcogen atoms in these two materials. Summing together, both JDOS and the Grüneisen parameters suggest that BOX has a shorter scattering lifetime than BiCuOX. This is actually confirmed in our numerical calculations (Figure 4 in main text and Figure 14).