Lattice vibration modes of the layered material BiCuSeO and first principles study of its thermoelectric properties

BiCuSeO has recently been shown to be one of the best oxide-based thermoelectric materials. The electrical properties of this material have been widely studied; however, the reasons for its intrinsically low thermal conductivity have only been briefly discussed. In this paper, we calculated the band structure and the electrical properties of BiCuSeO. The phonon spectrum, mode Grüneisen parameters and the thermal properties were also investigated. Additionally, we proposed a new method for illustrating the interlayer interactions in this material. For the first time, using first principles calculations, we provide direct evidence of the structural in-layer and interlayer off-phase vibration modes, which contribute to the anharmonic vibrations and structural scattering of phonons and result in an intrinsic low lattice thermal conductivity for BiCuSeO.


Introduction
Thermoelectric materials can achieve direct mutual conversion of thermal and electrical energy. This property provides these materials with potential applications in energy recovery and temperature control. The thermoelectric conversion efficiency can be characterized by the dimensionless figure of merit ZT S 2 σ κ = , where S, σ, κ, and T are the Seebeck coefficient, electrical conductivity, thermal conductivity, and absolute temperature, respectively. This simple equation provides researchers with several clear objectives, such as increasing the power factor and reducing the thermal conductivity. However, these properties are structure dependent, and they are difficult to improve simultaneously. Thus, to obtain promising thermoelectric materials, one can either search for materials with high electrical properties, find approaches to reduce the thermal conductivity, or optimize the power factor of materials with intrinsically low thermal conductivities.
Recently, the layered compound BiCuSeO has been found to exhibit fascinating thermoelectric properties. It was reported that the figure of merit of fine pristine BiCuSeO can reach 0.70 at 773 K [1]. Moreover, as a p-type semiconductor, the electrical properties of BiCuSeO can be improved by doping with divalent elements such as Mg, Ca, Sr and Ba [2][3][4][5][6]. By doping with the heavy metal Ba, the ZT of BiCuSeO can reach 1.1 at 923 K. By introducing Cu deficiencies, the electrical conductivity can be improved while maintaining the Seebeck coefficient and low thermal conductivity, resulting in a ZT of 0.81 at 923 K [7]. Liu et al investigated the influence of Pb substitution at the Bi sites and achieved a ZT of 1.14 at 823 K [8]. Sui et al synthesized textured BiCuSeO and obtained a ZT of 1.4 at 923 K, which is one of the best values achieved for this system [9].
From a theoretical perspective, the electronic structure of BiCuSeO was investigated by Stampler et al following the synthesis of pristine BiCuSeO samples. The density of states (DOS) and the band structure of BiCuSeO were compared with those of BiCuSO, and the lone pair-like ns electrons at the valence band maximum (VBM) were discussed [10]. The enhancement of the power factor was thoroughly explained. Barreteau et al reported an experimental investigation of the structural and electronic transport properties of Srdoped BiCuSeO. The electronic band structures and the projected DOS were investigated using DFT calculations, which revealed that a mixture of heavy and light hole pockets in the upper valence band is favorable for thermoelectric performance [5]. Zou et al recently reported the electronic structures and transport properties of BiCuChO (Ch = S, Se, Te) materials using first principles calculations and semi-classical Boltzmann transport theory, and they explored the relationship between their band structures and transport behavior; moreover, a prediction of the ZT was also provided [11].
Based on all of the aforementioned studies, the electrical properties have been thoroughly investigated, and the enhancement of the power factor has been well explained. However, the most important feature of BiCuSeO, its intrinsically low lattice thermal conductivity, was only briefly discussed, and the thermal conductivity used in predicting the figure of merit was fitted from experimental data. As it is known, phonons play the most important role in heat transfer and shed light on tuning the thermal conductivity. For instance, the leading thermoelectric material PbTe exhibits giant anharmonic phonon scattering [12], and resonant bonding leads to low thermal conductivity in rock salt IV-VI compounds such as SnTe and Bi Te 2 3 [13]. Moreover, a single crystal of the newly emerged thermoelectric material SnSe achieved an unprecedented ZT of 2.6, which also originated from its structural dissipation of phonon transport [14].
In this paper, we report the electronic transport properties, the phonon spectrum, and mode Grüneisen parameters of BiCuSeO and estimate its thermoelectric properties using first principles calculations, semiclassical Boltzmann transport theory and the quasi-harmonic approximation. Then, we explore the relationship between phonon behavior and thermal conductivity. In particular, the structural off-phase vibration modes of in-layer and interlayer interactions in BiCuSeO are revealed, which play a crucial role in scattering phonons and eventually lead to the intrinsically low thermal conductivity of this material.

Computational details
Density functional theory (DFT) based on first principles was used as the calculation method in the Vienna ab initio Simulation Package (VASP), BoltzTraP and Phonopy [15][16][17][18]. Previous calculations using the GGA+U method reported a band gap of 0.38 eV [11], which is an underestimation of 0.42 eV compared to the experimentally measured value of 0.8 eV [19] because of insufficient consideration of the exchange-correlation effects of the localized Cu 3d electrons [20]. Therefore, the B3LYP hybrid functional was used in our electronic structure calculations. The B3LYP hybrid functional is a hybrid of exact Hartree-Fock exchange with local gradient-corrected exchanges and correlation terms [21,22]. E xc is given by Δ is Beche's gradient correction to the exchange functional and E c PW91 Δ is the Perdew-Wang gradient correction to the correlation functional [23]. The choices of coefficients are a 0.2 0 = , a 0.72 x = , and a 0.81 c = , which means that the exchange energy consists of 80% LDA exchange and 20% non-local HF exchange. 72% of the gradient correction of the Becke88 exchange function is included [24]. The correlation energy consists of 81% LYP [25] and 19% of the VWN correlation functional [26], which is fitted to the correlation energy in the random phase approximation of the homogeneous electron gas.
The transport properties were obtained using semi-classical Boltzmann theory in conjunction with the rigid band approximation. The relaxation time was considered to be an energy-independent constant, and it was determined by comparing the calculated σ τ value with the measured σ, as reported in previous studies [11]. The phonon spectrum was calculated by adding finite displacements to atoms in BiCuSeO supercells. The force constants and eigenvectors were then calculated to obtain the phonon vibration frequencies. The atomic motions were subsequently calculated using Phonopy and were visualized through analysis software called 'v_sim' [27]. The bulk thermal properties were determined using the quasi-harmonic approximation. The volume dependence of the phonon frequencies was introduced and was attributed to the anharmonic effect. At constant pressure, the Gibbs free energy is defined as where, min V means finding a unique minimum value of the function in the brackets by changing the volume. The volume dependencies of the energies in the electronic and phonon structures are different. With varying temperature, the volume dependence of the phonon free energy changes, and then the equilibrium volume changes with the temperature [18]. Once equilibrium volumes are known, a number of thermodynamic properties can be derived, namely, the Helmholtz free energy (F), the equilibrium entropy (S), the internal energy (U), and the constant-volume heat capacity (C v ) [28]: The Grüneisen parameter can be computed from the volume derivative of TS − : The mode Grüneisen parameters are also calculated from first principles calculations at three different volumes to explore the anharmonic contributions. For a phonon with wave vector q and band index ν, the mode Grüneisen parameter q ( ) γ ν is given by the following approximation [18]: were V is the volume, q ( ) ω ν are the phonon frequency, D(q) is the dynamic matrix, and e q ( ) ν are the eigenvectors at the equilibrium volume. The remaining two volumes are for D q ( ) Δ with slightly larger and smaller volume than V. In our calculation, we used a volume difference of approximately 2 percent.
The velocity autocorrelation function (VAF) was obtained by ab initio molecular dynamics (AIMD) [29] using the GGA of the PBE form as implemented in the VASP code with the canonical ensemble. The PAW method is adopted, and a 3 × 3 × 2 superlattice with 216 atoms was used. The calculation ran for more than 2000 ps, and the simulation temperature was set to 20 K to estimate the intrinsic behavior.
In all of these calculations, we used a plane-wave energy cutoff of 600 eV and an energy convergence criterion of 10 −8 eV. The B3LYP calculations were conducted using a 4 × 4 × 4 k-mesh for the complex computational process. A 39 × 39 × 19 k-mesh was used to accurately determine the transport properties. For phonon spectra and other thermal property calculations, 6 × 6 × 4 k-meshes were used because of the use of the superlattice. For AIMD, a 1 × 1 × 1 k-mesh was employed.

Structural information
BiCuSeO materials crystallize in the P nmm 4 space group with a ZrSiCuAs-type structure [30,31]. Figure 1 shows the crystal structure of BiCuSeO.
The antifluorite (Cu Ch ) 2 2 2− chalcogenide layers, which are composed of slightly distorted edge-sharing CuCh 4 tetrahedra, alternate with fluorite-like (Bi O ) 2 2 2+ oxide layers and are stacked perpendicularly along the c axis of the tetragonal unit cell [10,30,32]. The Cu-Se bonds are covalent, Bi-O bonds are ionic and the bonding between layers is ionic because the layers are not electroneutral. The calculations were performed on the primitive cell containing eight atoms (two atoms for each element). To determine the optimized structure, we manually selected the lattice parameter, and fixed the shape while relaxing the ions inside the structure. For each value of a, the value of c increased evenly, forming a series of structures, and a mesh of lattice parameters could be obtained. Calculations were conducted for each structure. GGA+U [33] calculations were employed to find the optimized structure because the B3LYP hybrid functional required a complex computational process; however, it achieved similar performance in estimating the lattice parameters [34]. The Coulomb parameter U and exchange terms (screened) J were chosen according to a previous study [35]. The total energy was the criterion for optimizing the structure, and the structure with the lowest energy was considered to be the most stable structure. In figure 2, the optimized structure is highlighted in red, and the optimized lattice constants were a = 3.955 Å, c = 9.123 Å.

Electrical properties
The calculated electronic band structure of BiCuSeO and the density of states (DOS) are shown in figure 3. BiCuSeO is a direct-band gap semiconductor with both the conduction band minimum (CBM) and the valence band maximum (VBM) at point Z. The band gap calculated using B3LYP is 0.62 eV, which is a quite large correction to the GGA+U value of 0.38 eV but still smaller than the experimentally estimated value of 0.8 eV, indicating that the exchange-correlation effects of the localized Cu 3d electrons were still not sufficient. The deviation is not surprising because the value of the band gap of all bulk materials cannot be simply estimated using only one set of the fraction of exchange and correlation in the hybrid functional; the band gaps for real bulk materials are complicated and deviation occurs when approximations are adopted.
BiCuSeO tends to form p-type semiconductors [10,19]; thus, the upper part of the valence band has the greatest effect on the electrical properties. As shown in figure 3, the dispersion along the M Γ − line is greater than that along the Z R Γ − − line, which represents the mixture of heavy and light bands in the upper valence band of BiCuSeO. The light band provides good electrical conduction, and the heavy band is favorable for a high Seebeck coefficient, which has been proven by previous studies [11]. The DOS near the band gap is shown in figure 3; due to the complex computational process, only a 4 × 4 × 4 k-mesh is used for the band structure and DOS calculations. Thus, the precision of the DOS is low, and only a qualitative analysis is meaningful. As shown  in figure 3, although rapid changes in the DOS represent a large Seebeck coefficient [11,36], the carrier density is rather small around the Fermi level, which results in a moderate power factor. Doping divalent elements at the Bi sites is promising for optimizing the electrical conductivity and obtaining good thermoelectric properties.
In general, the underestimation of the band gap needs to be corrected to obtain more realistic electrical properties. As has been previously reported, in the calculation of the electronic transport properties, the discrepancy is corrected through the use of a scissors operator, which moves the conduction band upwards to match the experimental value of the band gap [11]. The electrical properties of BiCuSeO as a function of temperature are shown in figure 4. The calculated Seebeck coefficient increases when the temperature increases, and decreases as the doping level increases. Because the underestimation of the band gap has been corrected, consistency between the calculated results and the experimental data [2,8] is expected. The agreement between the calculated and experimental electrical conductivity is good in the high temperature regime, whereas deviation occurs in the low temperature regime because the relaxation time is fitted using high temperature results. The electrical conductivity deceases with increasing temperature, and the electrical conductivity increases as the level of doping increases at the same temperature, as shown in figure 4(b).

Vibration modes analysis and thermal properties
To investigate the thermal properties of BiCuSeO, we first calculated the phonon spectrum using a 3 × 3 × 2 superlattice, as shown in figure 5. The primitive cell of BiCuSeO contains eight atoms; thus, the calculated phonon spectrum has 24 dispersion relations. The phonon spectrum reveals the lattice vibration modes of BiCuSeO. By analyzing the atomic vibrations provided by Phonopy, the motion of each atom can be distinguished, and the interactions between atoms can be described in detail. The upper six bands ranging from 7 to 14 THz are the spectrum of the oxygen atoms, including four transverse branches that consist of two acoustic branches and two optical branches each. Additionally, two longitudinal branches consist of one acoustic branch and one optical branch. The lower eighteen bands ranging from 0 to 5 THz show strong interactions between the acoustic and optical branches, which indicates strong anharmonicity. Moreover, the phonon spectrum possesses relatively flat curves, indicating a small phonon group velocity based on the formula , and is beneficial to the low thermal conductivity. The vibration modes were analyzed by carefully investigating the atomic motion for each band of the phonon spectrum, especially at the high symmetry points. The index of each band was labeled in figure 5. For instance, at the G point, the lowest three bands exhibited pure normal modes, as shown in figure 6(a), whereas the other bands at G points exhibited different motions, such as Cu and Se atom layers moving associatively or Cu atoms moving while the other atoms remained fixed. Additionally, for the fourth to sixth bands at the R point, the Cu-Se atom layer moved as a unit, and the Bi-O atom layer moved with a different phase, which was believed to be direct evidence for the interlayer interaction, as shown in figure 6(b). Interestingly, Bi atoms and Se atoms with the same x and y coordinates tended to move associatively, forming a vibration pair and which occurred in almost every band, as shown in figure 6(c). The same motion also occurred with Cu atoms and O atoms. Moreover, bands with degenerate energy exhibited the same atom motions with different phases or directions.
At the same time, several peculiar vibration modes were detected, which can be attributed to the presence of Cu atoms [37]. For example, at the ninth and tenth bands of the R point or the seventh and eighth bands of the A point, Cu atoms exhibited uncommon circular motion, as shown in figure 6(d). The arrows in this figure are not indicators of the initial vibration direction but rather to show an deviation from normal vibration. Moreover, a large vibration amplitude of Cu was also observed at the 13th band of the Z point and also occurs at other points. The presence of these peculiarities in the vibration modes may be attributed to the low thermal conductivity.
As it is known, strong anharmonicity in bonding contributes to low lattice thermal conductivity in materials [38,39]. The Grüneisen parameter is often considered to be an anharmonicity parameter that reflects how much the phonon vibrations in a crystal lattice deviate from harmonic oscillations. Therefore, mode Grüneisen parameters were investigated using the quasi-harmonic approximation to clarify the origin of the intrinsically low thermal conductivity. As shown in figure 7, the red squares correspond to the mode Grüneisen parameters of acoustic branches, black circles correspond to the frequencies between 2.5 and 6 THz, and blue diamonds represent the phonon modes of O above 6 THz. The mode Grüneisen parameters vary between 0.08 and 6.74. The average Grüneisen parameter is approximately 2.5, which is considerably larger than the previous value of 1.5 estimated from the average sound velocity [4]. Comparably, the single crystal of the newly emerged thermoelectric material SnSe has a Grüneisen parameter of 2.1 in the b axis, and its thermal conductivity is 0.35 W mK [14], which is similar to that of BiCuSeO. As has been previously reported for BiCuSeO, the anharmonicity of the chemical bond drives the phonon-phonon Umklapp and normal processes that limit the lattice thermal conductivity. The presence of the heavy element Bi also contributes to the large Grüneisen parameter due to its large valence shell formed by lone pair electrons [4]. As can be seen from figure 7, the mode Grüneisen parameters of the low-frequency modes are very large, which indicates the strong interaction between  acoustic phonons and optical phonons. These phonon-phonon Umklapp processes can significantly scatter acoustic phonons carrying heat and reduce the lattice thermal conductivity. The mode Grüneisen parameters between 2.5 to 6 THz are around 2 to 2.5 for most branches with higher contribution from heavy element Bi around the G and Z points. The mode Grüneisen parameters for frequencies larger than 6 THz come from O, and should be small because O is light. However, its bonding with Bi increases its contribution to the mode Grüneisen parameters. In contrast, the O-O bonds have little contribution to the anharmonicity, resulting in mode Grüneisen parameters lower than 1.
The anomalously high Grüneisen parameters of BiCuSeO also reflect its crystal structure, which contains distorted edge-sharing CuCh 4 tetrahedra and a layered structure with zigzag geometry. This implies a soft lattice, which can dissipate phonons easily.
Using previously obtained information, the thermal conductivity is calculated. In general, the total thermal conductivity is the sum of the lattice thermal conductivity l κ and the electronic thermal conductivity e κ . Additionally, the lattice thermal conductivity can be estimated using the following relationship [40].
where n is the number of atoms in the primitive unit cell, 3 δ is the volume per atom, θ is the Debye temperature, M is the average mass of atoms in the crystal, and A is a collection of physical constants (A 3.1 10 6 ≈ × − if κ is in W mK, M in atomic mass unit, and δ in Å). The γ here is the temperature-dependent average Grüneisen parameter and it can be obtained from the quasi-harmonic approximation. The temperature dependence of the Grüneisen parameter is shown in figure 8. One can observe that the Grüneisen parameter tends to increase with increasing temperature, which causes an increase in the degree of thermal expansion, as shown in the inset of figure 8.
With the obtained Grüneisen parameter, the lattice thermal conductivity can be calculated using equation (8), as mentioned above, whereas the electronic thermal conductivity is calculated using the Wiedemann-Franz relation L T , where σ is the electrical conductivity and L is the Lorenz number which can be determined according to previous studies [2,41]. Figure 8 presents the results from the calculation of thermal conductivity. At high temperature, the calculation results are comparable with the experimental data, and the deviation arises from point defects or grain boundary scattering in real crystals, which are not included in the calculations. At low temperature, the experimental data are far below the calculated results, even after the addition of scattering contributions from the defects, which indicates intrinsic scattering effects that are not included in the quasi-harmonic approximation. This effect is favorable for the intrinsically low thermal conductivity of BiCuSeO, which will be discussed in detail in the next section.

Inlayer and interlayer off-phase dissipation
To explain the origin of the discrepancy in the previous section, the velocity autocorrelation function (VAF) as a function of the correlation time for each element was calculated via AIMD simulations. The VAF Z ( ) τ is defined as follows: where N is the number of atoms, t 0 is the simulation starting time, t max is the selected ending time and τ is the time difference, which ranges from 0 to t t max 0 − . The total simulation time used in the calculations is t max + range of Z ( ) τ . Z ( ) τ is an average over time and over the trajectories of all atoms. The self-correlation functions were normalized by v t ( ) 2 0 . The self-velocity correlation functions of Bi, Cu and Se are shown in figure 9. In general, the correlation function and the variation in the velocities of atoms have the same period of oscillations; thus, the decay of the correlation function can exactly reflect the decay in the correlations in atomic motion along the trajectories of the atoms. Moreover, a damped oscillating VAF generally represents strong interactions in solids, such as Se elements, indicating strong interactions between them, and their vibrations are very much localized [37]. However, the Bi and Cu elements behave differently; typical solid-like features were not observed, but the observed features were far from those of liquids [42]. We therefore analyzed the power spectrum (PS) of the VAF, which is defined as    In previous studies, the power spectra of VAF were discussed to distinguish the states of matter. The PS for the solid-like phase has distinct peaks, and the lack of translational motion [intensity I ( 0) 0] ω = = is a wellaccepted attribute of solid-like behavior. In contrast, the PS for the liquid-like phase has a peak that is less pronounced than that for the solid-like phase [42,43]. Here, we used the power spectrum in an innovative way, that is, to elucidate the vibration phase difference between atoms in a layered structure, and consequently, the interaction between them. The power spectra of the self-VAF of Bi, Cu, and Se and the cross-VAF of Bi/Se and Cu/Se are shown in figure 10. From the PS, we can determine the frequencies of different atomic motions for both the self-VAF and cross-VAF. Additionally, the relative intensities of the self-VAF and cross-VAF represent the frequency distributions of atomic motion and the phase difference of the motion between the atoms in the BiCuSeO layered structure, respectively. As shown in figure 10, the PS of Bi consists of four main peaks, namely, 0.28, 0.47, 0.75, and 0.93 THz. The PS of Cu contains three positive main peaks at 0.37, 0.75 and 0.93 THz. Moreover, two positive peaks were observed at 0.37 and 0.93 THz for Se. The PS of the cross-VAF for Bi/Se and Cu/Se are also shown in figure 10. If both elements vibrate with the same frequency, then their peaks in the PS of self-VAF will appear at the same frequencies in the cross-VAF spectrum, which can be easily understood as the superposition of the 'wave functions' of different atoms. However, the PS of Bi/Se has negative peaks at 0.37 and 0.84, which do not exist for the PS of single types of atoms, indicating different phases for the Bi atoms and Se atoms. Similarly, negative peaks in the PS of Cu/Se are located at 0.47 and 0.75 THz. Considering the relative positions of the Bi, Cu and Se atoms, the negative signals of the cross-VAF of Cu/Se can be attributed to the in-layer vibration modes with a reverse phase, and those of Bi and Se can be considered to be proof for the interlayer off-phase interaction. The results are consistent with the visualized atomic motions mentioned above.
The aforementioned in-layer and interlayer off-phase atomic motions can significantly increase the possibility of scattering the propagation of the in-plane phonons, which could eventually reduce the lattice thermal conductivity. More importantly, such behavior originates from the intrinsic structural characteristics. Thus, at low temperature, the dissipation significantly reduces the lattice thermal conductivity. As the temperature increases, the thermal fluctuations begin to dominate the dissipations of the phonon propagation, which leads to small dissipations because the negative off-phase can be concealed under thermal noise. Consequently, the scattering effect would decrease at high temperature. As indicated by the previous calculation, at low temperature, the underestimation of such scattering could result in a difference of approximately 1. − − with deviation due to other scattering effects of phonons originating from the grain boundaries, defects and other factors in real polycrystalline BiCuSeO; for our calculations, the material was considered to be a perfect single crystal.

Conclusion
In summary, the band structures, the phonon spectrum, the mode Grüneisen parameters, the thermoelectric properties and the structural scattering effects of BiCuSeO were investigated using first principles calculations. The electronic structure indicates that the thermoelectric properties can be enhanced by p-type doping at the Bi sites. Strong anharmonic bonding of the heavy element Bi and abnormal atomic motion originating from the layered structure caused a large anharmonicity and resulted in a Grüneisen parameter of 2.5 at room temperature. The structural off-phase in-layer and interlayer interactions significantly scattered phonons carrying heat, which was hypothesized to be the origin of the overestimation of the lattice thermal conductivity using the quasi-harmonic approximation. The significance of this paper is that for the first time we highlighted in detail the vibration modes concerning the layer interactions in a layered structure and established the connection between these modes and the phonon scattering. The results revealed that an intrinsic scattering pattern contributes to the low lattice thermal conductivity in layered materials, which could be of considerable interest to recently developed layer-structured thermoelectric materials.