A first-principles study on the phonon transport in layered BiCuOSe

First-principles calculations are employed to investigate the phonon transport of BiCuOSe. Our calculations reproduce the lattice thermal conductivity of BiCuOSe. The calculated grüneisen parameter is 2.4 ~ 2.6 at room temperature, a fairly large value indicating a strong anharmonicity in BiCuOSe, which leads to its ultralow lattice thermal conductivity. The contribution to total thermal conductivity from high-frequency optical phonons, which are mostly contributed by the vibrations of O atoms, is larger than 1/3, remarkably different from the usual picture with very little contribution from high-frequency optical phonons. Our calculations show that both the high group velocities and low scattering processes involved make the high-frequency optical modes contribute considerably to the total lattice thermal conductivity. In addition, we show that the sound velocity and bulk modulus along a and c axes exhibit strong anisotropy, which results in the anisotropic thermal conductivity in BiCuOSe.

relatively rare. Recently, ab initio calculation of lattice thermal conductivity has been developed by combining the first-principles calculations of interatomic force constants with solving the phonon Boltzmann transport equation iteratively 27,28 . The first-principles calculations can accurately predict the phonon thermal conductivity without any assumption on phonon lifetimes, and capture the transport properties of each phonon, and have been used in the investigations of phonon transport properties for many materials [28][29][30][31][32][33][34][35][36][37][38][39][40] .
In this work, we present a detailed investigation of phonon transport property for BiCuOSe to understand the physical mechanism of ultralow thermal conductivity and strong anisotropy in phonon transport. Our calculations reproduce the lattice thermal conductivity of BiCuOSe. The grüneisen parameter, which measures the anharmonicity of bonding, is 2.4 ~ 2.6 at room temperature, indicating a strong anharmonicity in BiCuOSe. This leads to the low thermal conductivity. Moreover, we find that the high-frequency optical phonons of BiCuOSe contribute considerably to the total thermal conductivity, which is remarkably different from usual picture with little contribution from optical phonons. Additionally, our calculations show that the obvious anisotropy of phonon conduction originates from the anisotropic sound velocities in BiCuOSe.

Results and Discussions
Crystal structure. BiCuOSe has tetragonal structure with space group of P4/nmm 3  To determine the equilibrium lattice parameters, several volumes around the expected equilibrium volume are used to be relaxed for the ions and shapes. And a set of volumes V i and energies E i are obtained. As shown in Fig. 2(a), we make a fit of these ( , ) V E i i s to the Birch-Murnaghan 3rd-order equation of state 41 . Table 1 presents the calculated lattice constants, atomic positions, bulk modulus, and pressure derivative of bulk modulus for BiCuOSe. We also plot the relative change of lattice constants with respect to volume in Fig. 2(b). It shows that in BiCuOSe, the a axis is less sensitive to pressure or temperature than the c axis. The calculations are performed at zero temperature, and the experimental data are obtained at room temperature. Regardless of the thermal expansion, the HSE06 method gives a good prediction for the equilibrium volume within an error less than 0.4%, while the PBE method overestimates the equilibrium volume of BiCuOSe with an error more than 3%. However, now it is too expensive to perform the calculations of force constants by using the HSE06 method implemented in VASP, because the large amount of supercells with large number atoms are employed in the calculations. Then we use the PBE method to calculate the phonon property and thermal conductivity of BiCuOSe. Nevertheless, we could estimate roughly the errors of lattice thermal conductivity obtained by the PBE method. From Table 1, the bulk modulus by PBE method is about 14% less than that by the HSE06 method. Because in solids, the velocity of elastic wave, ṽ B. The lattice thermal conductivity κ~v 2 . And then κ~B. The PBE method may underestimate the lattice thermal conductivity by around 14%. Phonon transport properties. The calculated dielectric constants and Born effective charges are listed in Table 2. Our PBE calculations of dielectric constants and Born effective charges are slightly larger than those by PBEsol calculations of Saha 26 . Because BiCuOSe has layered tetragonal structure with space group of P4/nmm, the Born effective charges and dielectric constants exhibit anisotropy with different in-plane and out-of-plane values. The dynamical matrix for calculation of phonon spectra of BiCuOSe can be constructed from the harmonic force constants by adding a non-analytical correction with dielectric constants and Born effective charges, where Z and ε ∞ are the Born effective charge and dielectric constant, respectively. And the dynamical matrix where n denotes the nth primitive cell in supercell and R n a translation vector for the nth primitive cell, τ and τ M refer to the atom and its mass in the primitive cell, α is the Cartesian components of x, y, or z, and Φ τα τ α ′ ′ is the harmonic force constants. The parameter ρ (= 0.25) is set so that the non-analytical term becomes negligible at the zone boundaries. Figure 3 presents the phonon spectra along several high symmetry lines in Brillouin zone and phonon density of states. The primitive cell of BiCuOSe contains eight atoms, and there are 24 phonon branches in the phonon spectra. As shown in Fig. 3, the high-frequency optical modes are separated from the low-frequency modes by a gap of 46 cm −1 . Because of the different directions of Γ -M and Γ -Z, the non-analytical term describing the LO-TO splitting effect bring different corrections to the phonon modes near Γ point. Then the phonon spectrum seems to be not continuous at Γ point. From Fig. 3, the high-frequency optical modes above 213 cm −1 are mainly from O vibrations, and these modes exhibit obvious dispersion, indicating that they have relatively large group velocities. While the acoustic and low-frequency optical modes are mainly from the vibrations of Bi, Cu, and Se atoms.  The phonon lifetimes of BiCuOSe at room temperature are presented in Fig. 4. We calculate the three-phonon scattering rates and the contribution of isotopic disorder to scattering probabilities. Figure 4 shows that compared to the anharmonic phonon-phonon scattering, the isotopic scattering τ ( / ) 1 isotopic is negligible to the thermal resistance. There is a slight difference between the lifetimes by relaxation time approximation (RTA) and those of exact numerical solutions of the phonon Boltzmann equation, which means that the RTA describes the lattice thermal conductivity of BiCuOSe quite well. It implies that the umklapp scattering is dominant for the thermal resistance, and the redistribution of phonons resulting from the normal scattering processes has little influence on the phonon transport of BiCuOSe. The lifetimes of most acoustic and optical phonons with frequencies lower than 60 cm −1 are several ps' , and for high-frequency optical modes, they lie in 0.1 ps to 1 ps in BiCuOSe. Such frequency-dependent lifetimes have the same magnitude with those of PbTe 34 , while two orders of magnitude lower than those of Si 32 . Additionally, in the low-frequency range, the lifetimes show ω −2 dependence, which agrees with Klemens' prediction 42 . While the lifetimes of high-frequency optical modes (above 213 cm −1 ) behave one order of magnitude higher than those predicted by Klemens' relation. Figure 5 shows the average temperature-dependant lattice thermal conductivity of BiCuOSe from 100 to 800 K, and as well as the contributions from high-frequency optical modes to the total lattice thermal conductivity. The average lattice thermal conductivity is calculated by κ κ κ κ . As shown in Fig. 5, the thermal conductivities of BiCuOSe predicted by RTA are slightly lower than those of iterative solutions of the phonon Boltzmann equation. The calculated κ L at 300 K is 0.8 Wm −1 K −1 , which is consistent with experimental results 2,4-6,10 . We note here that the PBE calculations may underestimate the thermal conductivity by around 14%. On the other hand, we merely consider the intrinsic phonon-phonon scattering, and ignore other scattering source such as crystal interfaces, and may overestimate the thermal conductivity. Although BiCuOSe exhibits ultralow intrinsic thermal conductivity, the lattice thermal conductivities obey the ~1/T relation at high temperature (above Debye temperature ~240 K). In addition, we find that the contribution of high-frequency (above 213 cm −1 ) optical phonons, which are mostly from vibrations of O atom, to overall lattice thermal conductivity is larger than 1/3. This is remarkably different from the usual picture with very little contribution of thermal conductivity from such high-frequency phonons in most bulk materials 43 .
To explore the mechanism of ultralow thermal conductivity, we continue to discuss the strength of phonon-phonon scattering processes. At high temperatures, the phonon lifetime in crystal is determined by phonon-phonon scattering processes. A useful measurement of the strength of phonon-phonon scattering is the mode Grüneisen parameter γ ( ) q p , which can be calculated by  where V 0 is the equilibrium volume. Figure 6 presents the mode Grüneisen parameters with respect to frequencies calculated by Eq. (3) and (4) for BiCuOSe. Eq. (3) and (4) give consistent results for the mode Grüneisen parameters of acoustic and low-frequency optical phonons. While for high-frequency optical phonons, Eq. (3) predicts some larger mode Grüneisen parameters than those calculated by Eq. (4). And the discrepancy between them becomes less and less with increased frequencies.
As shown in Fig. 6, throughout the Brillouin zone, the mode Grüneisen parameters of BiCuOSe are mostly positive. Moreover the maxima of mode Grüneisen parameters are very high near Γ point and at low frequency range, up to 14. We further average γ ( ) q p to obtain the average Grüneisen parameter with temperatures respectively. The calculated γ ave of BiCuOSe is larger than that of PbTe, of which the γ ave is 1.96 ~ 2.18, which was also obtained by first-principles calculations 44 . It is known that PbTe has very strong anharmonic phonon scattering 45,46 . Therefore, BiCuOSe behaves even more anharmonic than PbTe, and this leads to its ultralow intrinsic phonon conductivity.
In common materials, the high-frequency optical phonons exhibit small group velocities and suffer large scattering, thus make little contribution to the thermal conductivity. In contrast to usual picture, the high-frequency optical modes of BiCuOSe make remarkable contribution to the total thermal conductivity. To compare the  The JDOS is a representative of the phase space available for scattering events. + Q corresponds to the absorption processes, and − Q is for the emission processes. We present the JDOS of BiCuOSe in Fig. 7, which shows that both absorption and emission processes dominate for low-frequency phonons, while the high-frequency phonons involve much less scattering. Additionally, the high-frequency phonons are mainly from the vibrations of O atom. Due to the light atomic mass and strong interaction with Bi, the vibrations modes related to O atom exhibit strong dispersion, and their frequency range expands from 210 cm −1 to 460 cm −1 . These high-frequency modes have relatively high group velocities. Therefore, the high-frequency optical phonons make considerable contribution to the total thermal conductivity in BiCuOSe. Figure 8 presents the temperature-dependent lattice thermal conductivity along a (κ xx ) and c (κ zz ) axes of BiCuOSe, which shows strong anisotropy. The ratio of κ xx to κ zz is 2.3 ~ 3. Recently, Saha has proposed that the different Grüneisen parameters between in-plane and out-of-plane directions resulted in the anisotropy of thermal conductivity 26 . However, in Slack's expression 47 , which serves as estimation of lattice thermal conductivity, κ = To clarify the mechanism of anisotropy of the thermal conductivities, we firstly discuss the anisotropy of group velocity (GV) in BiCuOSe. Figure 9 presents the group velocities along Γ -X and Γ -Z, which shows that the GVs along Γ -X are much larger than those along Γ -Z, especially for high-frequency optical modes. The GV of longitudinal acoustic (LA) branch along Γ -X near zone center is about 3650 m/s, and the GV of transverse acoustic (TA) branches are 2138 and 1576 m/s. Whereas the GVs of LA and TA along Γ -Z near zone center are about 3354 and 1634 m/s. For high-frequency optical modes, the maximum GV along Γ -X reaches 3284 m/s, while that along Γ -Z is only 442 m/s. Therefore, the GVs along in-plane and out-of-plane directions are very different.

Anisotropy of the thermal conductivities.
Because κ~B, to give a quantity estimation for the ratio of in-plane to out-of-plane GV, we calculate the bulk modulus B along a and c axes. The elastic constants are calculated and tabulated in Table 3. The elastic constants satisfy all the Born's mechanical stability criteria, indicating that the system of BiCuOSe is in a mechanical stable  state. Using calculated elastic constants, the bulk and Young's modulus can be obtained by Voigt-Reuss-Hill approximations 48 . The obtained bulk modulus from the elastic constants is 76.2 GPa, which is very close to that calculated from fitting the Birch-Murnaghan's 3rd-order equation of state. The Young's modulus is a measure of the stiffness of an elastic material, and reflects the strength of bonding in the covalent crystals. Our calculated Young's modulus of BiCuOSe is = .
E 79 6 GPa, which agrees well with the experimental results of = . E 76 5 GPa by Pei et al. 10 . And the relatively small Young's modulus indicates a totally weak covalent bonding in BiCuOSe. As is known, the elastic constants c 11 and c 33 reflect the stiffness-to-uniaxial strains along the crystallographic a and c axes, respectively, while the elastic constants c 12 , c 13 , c 44 , and c 66 are related to the elasticity in shape. From Table 3, the value of c 11 is much higher than c 33 , which means that BiCuOSe is stiffer for strains along the a axis than along the c axis. Furthermore, the ratio of bulk modulus along a axis to that along c axis can be calculated by 49 The obtained / = .
for BiCuOSe, which shows strong anisotropy in elastic property. And this value is close to the ratio of κ xx to κ zz . Therefore, our calculations show that the strong anisotropy of lattice thermal conductivity in BiCuOSe should be ascribed to its anisotropic velocities. The anisotropy of bulk modulus, which leads to anisotropic sound velocity in BiCuOSe, reflects the layered structure of BiCuOSe. And the relationship between such strong anisotropy and the layered structure and anisotropically chemical bonding in BiCuOSe remains to be uncovered.

Conclusion
In summary, we have employed first-principles calculations to investigate the phonon transport of BiCuOSe. The obtained phonon lifetime of BiCuOSe is mostly in the range from 0.1 ps to several ps, which is comparable to those of PbTe. Our calculations show that there is strong bonding anharmonicity in BiCuOSe, which leads to its ultralow lattice thermal conductivities. Due to the high group velocities and low scattering processes involved, the high-frequency optical modes make a considerable contribution to the total lattice thermal conductivity. The anisotropy of thermal conductivities is ascribed to the anisotropy of wave velocities along a and c axes, which may originate from the layered structure and anisotropic chemical bonding in BiCuOSe. Our work facilitates deep understanding of the phonon transport properties in such layered materials.

Methods
The calculations are based on density functional theory method in the generalized gradient approximation with the Perdew, Burke, and Ernzerhof functional (PBE) 50 , as implemented in the Vienna Ab initio Simulation Package (VASP), which employs a plane-wave basis 51,52 . The plane-wave energy cutoff is set to be 500.00 eV, and the electronic energy convergence is 10 −8 eV. During relaxations, the force convergence for ions is 10 −3 eV/Å. The PBE functional overestimates the equilibrium volume, and then underestimates the bulk modulus and group velocity of phonons. The hybrid functionals could be more reliable to obtain the equilibrium volume 53 and the bulk modulus. We then employ HSE06 54,55 to calculate the equilibrium lattice parameters and bulk modulus of BiCuOSe, then estimate the errors by the PBE calculations.
The phonon frequencies and velocities are calculated within the harmonic approximation. Parlinski-Li-Kawazoe method, which is based on the supercell approach with finite displacement method as implemented in the Phonopy package 56,57 , is employed to obtain the phonon dispersion, phonon density of states, and group velocities for BiCuOSe. To obtain convergent phonon property, in the calculations of harmonic interatomic force   constants, a × × 3 3 2 supercell of primitive cell containing 144 atoms is employed, and a Γ -centered × × 2 2 1 Monkhorst-Pack k-point mesh is used to sample the irreducible Brillouin zone.
The phonon Boltzmann transport equation solved iteratively as implemented in ShengBTE 28 is employed to study the phonon transport of BiCuOSe. The phonon lifetime can be solved numerically with an iterative method. The zeroth-order solution corresponds to the relaxation time approximation. The RTA can well describe the phonon conduction of some materials such as Si and Ge 31 , while it fails to predict the thermal conductivity of carbon-based materials such as graphene 40 . We also inspect how accurate the prediction by RTA is for the thermal conductivity of BiCuOSe.
To get the three-phonon scattering matrix elements for calculating the phonon lifetimes, the third-order force constants should be calculated firstly. In ShengBTE, the third-order force constants are calculated by a finite-difference supercell approach. For the force calculations, a × × 3 3 2 supercell of primitive cell containing 144 atoms is employed so that there is only negligible interaction between atoms in the center and at the boundary. To get convergent lattice thermal conductivity, we adopt the eighth nearest neighbors as the cutoff for third-order force interactions. And there is 716 supercells with displaced atoms needed to be calculated in VASP. During the calculations of force constants, only the Γ point is set for the k-point grid. And in the calculation of lattice thermal conductivity, a 15 × 15 × 9 q-point grid (180 inequivalent q points) is employed. Apart from three-phonon processes, the contribution to scattering probabilities from isotopic disorder can also be calculated in ShengBTE. After calculating the phonon lifetime, the lattice thermal conductivity can be obtained by where Ω is the volume of primitive cell, α and β are the Cartesian components of x, y, or z, k B is Boltzmann constant, and ω ( ) λ f 0 is the Bose-Einstein distribution function.