Weiliang Ma

Marie-Christine Record

Jing Tian

Pascal Boulet

Influence of the stacking sequence on layered-chalcogenides properties: First principle investigation of Pb 2 Bi 2 Te 5

The Pb 2 Bi 2 Te 5 compound has been reported in literature with two stacking sequences -Te-Pb-Te-Bi-Te-Bi-Te-Pb-Te-and -Te-Bi-Te-Pb-Te-Pb-Te-Bi-Te-labelled in this work as A and B, respectively. The electronic and the thermoelectric properties of the Pb 2 Bi 2 Te 5 compound with the 2 different stacking sequences have been determined from a series of first principles calculations using density functional theory (DFT). The related compounds PbTe and Bi 2 Te 3 have also been investigated for comparison. Different exchange-correlation functionals have been tested, w/o spin-orbit coupling, which has been found to have important effects. The elastic moduli, dielectric constants, Born effective charges, and phonon dispersion within the quasi-harmonic approximation have also been calculated and based on these calculations results, the thermal conductivity has been determined by solving the Boltzmann transport equation. Additionally, the QTAIM theory was employed to explain the differences in the properties of the 2 stackings. The most interesting compound for thermoelectric applications has been found to be Pb 2 Bi 2 Te 5 with the stacking B sequence. The highest zT values have been found to be 4.02 in the a-axis direction and 2.26 in the c-axis one.

Introduction

Due to their low lattice heat conductivity, complex layered chalcogenides such as ternary compounds in the quasibinary A IV C V I -B V 2 C V I 3 systems (A IV = Ge, Sn, Pb; B V = Bi, Sb ; C V I = Te, Se) are potential candidates for thermoelectric applications. [START_REF] Kanatzidis | The Role of Solid State Chemistry in the Discovery of New Thermoelectric Materials[END_REF] The existence of homologous series of layered compounds of the n(A IV C V I )-m(B V 2 C V I 3 ) type in this system [START_REF] Shelimova | [END_REF] offers large possibilities for optimizing the thermoelectric properties. Indeed by using different types of stacks as well as various combinations, these homologous series behave as compounds generating devices. Many n(A IV C V I )m(B V 2 C V I 3 ) compounds have been already investigated as potential thermoelectric materials [3][4][5][6] , phase change memories for optical data storage 6,7 and topological insulators (TIs) [8][9][10][11][12][13][14][15][16] .

In the PbTe-Bi 2 Te 3 quasi-binary system, six ternary compounds belonging to the homologous series n(A IV C V I )m(B V 2 C V I 3 ) have been evidenced. Five of them, namely PbBi 2 Te 4 , Pb 2 Bi 6 Te 11 , PbBi 4 Te 7 , PbBi 6 Te 10 , and PbBi 8 Te 13 with a trigonal symmetry, are located in the Bi 2 Te 3 rich side [17][18][19] , whereas Pb 2 Bi 2 Te 5 with a nine-layers hexagonal unit cell 18,20 is located in the PbTe rich side. Most of the studies reported in literature concerns the compounds belonging to the Bi 2 Te 3 rich side, e.g. thermodynamic properties of PbBi 4 Te 7 and PbBi 2 Te 4 21 , vibrational spectra of Pb 2 Bi 2 Te 3 , PbBi 2 Te 4 and PbBi 4 Te 7 22 , thermoelectric properties of PbBi 2 Te 4 4 and PbBi 4 Te 7 23 . The Pb 2 Bi 2 Te 5 compound with a stability that was even questioned 17 has been much less investigated. In 1969, Petrov et al. 18 synthesized Pb 2 Bi 2 Te 5 from PbBi 4 Te 7 using a long time annealing to evaporate Bi 2 Te 3 and characterized it by electronography. The unit cell of Pb 2 Bi 2 Te 5 consists in nine-atomic layers, which are made by inserting two PbTe into Bi 2 Te 3 according to the following sequence: -Te-Pb-Te-Bi-Te-Bi-Te-Pb-Te-, whereas the unit cell of PbBi 4 Te 7 consists in 12-atomic layers made of a 5-layers slab (-Te-Bi-Te-Bi-Te-) and a 7-layers slab (-Te-Bi-Te-Pb-Te-Bi-Te-). The nine atomic layers in Pb 2 Bi 2 Te 5 are linked by van der Waals interactions. Recently Chartterjee et al. 20 synthesized by a simple solution-phase method at low temperatures a two-dimensional Pb 2 Bi 2 Te 5 nano-sheet and characterized it by several methods including powder X-ray diffraction, inductively coupled plasma atomic emission spectroscopy (ICP-AES) and energy dispersive X-ray analysis (EDAX). In agreement with the previous results 18 , the unit cell of Pb 2 Bi 2 Te 5 is constituted of nine-atomic layers but with a different atomic sequence, namely -Te-Bi-Te-Pb-Te-Pb-Te-Bi-Te-.

Therefore the stability and the atomic layer stacking sequence for this compound are still questioned.

In the present work, by using first principles calculations and density functional perturbation theory (DFPT), we investigated the stability, the electronic structure, the thermoelectric properties and the thermal conductivity of Pb 2 Bi 2 Te 5 for each stacking sequence.

Computational procedure

The electronic structure calculations have been performed within the frame of density functional theory using allelectron FP-LAPW approach with local orbital method as implemented in WIEN2K 24 . Several exchange-correlation functionals have been used the detail of which will be mentioned when appropriate in the results section. For structural optimizations the Brillouin zone has been sampled with the k-meshes 8 × 8 × 8, 8 × 8 × 8 and 12 × 12 × 2, for PbTe, 0a Aix-Marseille University, University of Toulon, CNRS, IM2NP, Marseille, France. E-mail: weiliang.ma@etu.univ-amu.fr, Corresponding author: m-c.record@univ-amu.fr 0b Aix-Marseille University, CNRS, MADIREL, Marseille, France. E-mail: jing.tian@etu.univ-amu.fr, pascal.boulet@univ-amu.fr

Bi 2 Te 3 and Pb 2 Bi 2 Te 5 , respectively. For the subsequent convergence of the self-consistent energy the k-meshes have been set as 16 × 16×16, 16 × 16 × 16 and 18 × 18 × 4. The total energy and atomic forces convergence threshold have been defined as 0.136 meV and 0.257 meV/Å, respectively. The cut off energy for the core-valence separation has been fixed at -5.0 Ry for PbTe and Pb 2 Bi 2 Te 5 , and -6.0 Ry for Bi 2 Te 3 and the R mt K max value as 9.0. The structure and charge density calculated above are used to analyze the topological properties within QTAIM theory 25,26 . The thermoelectric properties have been computed by solving the Boltzmann transport equation within the constant relaxation time τ and the rigid band structure approximations as implemented in BoltzTraP2 package 27 . The implementation of BoltzTraP2 is based on the use of full bands structure in the Brillouin zone. The sampling, which is important in transport calculation, has been done using a very dense k-mesh of 48 × 48 × 48 for PbTe, 36 × 36 × 36 for Bi 2 Te 3 , and 64 × 64 × 14 for Pb 2 Bi 2 Te 5 . Within the limit of validity of the semi-classical linearized Boltzmann transport equation, the strongly band structure dependent Seebeck coefficient (S) can be directly obtained. By contrast, for the calculation of the electronic part κ e of the thermal conductivity and the electronic conductivity σ that depend on τ, the τ value has to be evaluated independently. For this, we have used the effective mass approximation with deformation potential (DP) theory to calculate. 28 Second and third order anharmonic force constants have been calculated by means of DFPT method 29 by using the QUANTUM-ESPRESSO package 30 together with the Phono3py 31 program. A plane-wave energy cutoff of 70 Ry (952 eV) has been employed. Total energies have been minimized with a convergence criterium of 10 -7 Ry and a total force threshold of 10 -4 Ry/Bohr. A supercell of 2 × 2 × 2 has been considered, which consists in a total of 64, 40, and 72 atoms for PbTe, Bi 2 Te 3 and Pb 2 Bi 2 Te 5 , respectively, with 4 × 4 × 4, 6 × 6 × 6 and 4 × 4 × 1 q-mesh sampling. In subsequent post-processing, phonon lifetimes have been sampled with a finer 19 × 19 × 5 mesh.

Structural information

The bulk crystal structure of PbTe is a face-centered cubic one (space group Fm 3m) with 2 atoms in the primitive cell (Fig. 1a). Bi 2 Te 3 crystallizes in the rhombohedral lattice system (R 3m) with 5 atoms in the primitive cell stacked along the c-axis and can also be described with a hexagonal unit cell (Fig. 1b). There are three equivalent atomic positions, which locate at Bi (µ, µ, µ), Te1 (ν, ν, ν) and Te2 (0, 0, 0). The structure of Bi 2 Te 3 consisting in three slabs containing five atom layers each can be easily visualized in the conventional hexagonal cell. The weak Te-Te bonds (373.3 pm 32 ) show that adjacent slabs are linked together through van der Waals forces that contribute to the crystal cohesion. The studied compound Pb 2 Bi 2 Te 5 , which can also be identified as the quasi-binary compound 2(PbTe)-Bi 2 Te 3 , crystallizes in the P3m1 space group with nine atoms per unit cell stacked along the c-axis. To get a direct view of the atoms stacking, a 1 × 1 × 2 supercell along the c-axis can be considered (see Fig. 1c,d). Among the nine atoms, five inequivalent atoms can be seen. As with Bi 2 Te 3 a weak Te-Te interaction is observed between slabs. As mentioned above, two possible atoms sequences have been reported in the literature 18,20 : -Te3-Pb-Te2-Bi-Te1-Bi-Te2-Pb-Te3-(stacking A) and -Te3-Bi-Te2-Pb-Te1-Pb-Te2-Bi-Te3-(stacking B); from one structure to the next the bismuth atoms and the lead ones are interchanged.

Results

Structural properties

The equilibrium lattice constants of PbTe, Bi 2 Te 3 and Pb 2 Bi 2 Te 5 have been calculated with three exchange-correlation functionals, namely LDA 33,34 , PBE 35 , and rev-vdW-DF2 36 with and without spin-orbit coupling (SOC). As there are weak Te-Te bonds in Bi 2 Te 3 and Pb 2 Bi 2 Te 5 , van der Waals forces have been considered in these compounds. The rev-vdW-DF2 36 has been shown to be accurate in weakly bound solid systems, such as rare gases and graphitelike structure, and also reliable in strongly bound systems. 37 The calculated lattice constants have been determined by fitting the Birch-Murnaghan equation and the results are listed in Table1 together with experimental data for comparison. Although the calculations on Bi 2 Te 3 have been performed with the rhombohedral primitive cell, in order to be consistent with the hexagonal structure of Pb 2 Bi 2 Te 5 we have given in Table 1 the hexagonal lattice parameters of Bi 2 Te 3 . The effect of SOC on the lattice constants is negligible as we can see from Fig. S1-4 (supplementary materials). For PbTe, the lattice constants are in better agreement with the experimental ones when calculated with PBE, whereas the agreement is better when calculated with rev-vdW-DF2 for Bi 2 Te 3 . The agreement is also better when using rev-vdW-DF2 for Pb 2 Bi 2 Te 5 stacking B. For Pb 2 Bi 2 Te 5 stacking A from PBE to rev-vdW-DF2 the c lattice parameter agreement is worsen whereas the a one is improved. This could suggest that the contribution of the van der Waals interactions is larger in stacking B. To allow for comparisons between forthcoming results on stacking A and B, the calculations have been performed with the rev-vdW-DF2 functional.

Electronic and thermal properties 4.2.1 Electronic band structure

Based on the rev-vdW-DF2 optimized structure, the band gaps have been calculated with six different functionals with and without SOC. The Engel-Vosko functional (EV-GGA) 40 allows for reproducing the real exchange potential and the Tran-Blaha-modified Becke-Johnson exchange potential (TB-mBJ) 41 for getting band gaps with a high accuracy similar to that obtained from GW calculations. The HSE hybrid functional, with 25% Hartree-Fock exchange, obviates the band-gap underestimation commonly observed with conventional GGA functionals. The results are presented together with literature values in Table 2. As can be seen, the effect of SOC is significant for all the compounds. The non-SOC energy gaps are overestimated compared with the experimental values. The inclusion of SOC yields closer values to experimental ones, in particular with the PBE functional. Hence this functional has been chosen for all subsequent calculations. Even if gaps are enlarged with HSE, it is still closed for Pb 2 Bi 2 Te 5 with stacking A, which reflects its metallic character. The band gap of Pb 2 Bi 2 Te 5 stacking B is about 0.113 eV, thus this structure can be considered as a narrow gap semiconductor. The gaps of Bi 2 Te 3 and Pb 2 Bi 2 Te 5 calculated with and without vdW corrections are close, which indicates that the effect of vdW on their electronic structures is negligible.

The calculated bands structures with and without SOC along high symmetry points of the first Brillouin zone are shown in Figure 2 for PbTe and Bi 2 Te 3 , and in Figure 3 for Pb 2 Bi 2 Te 5 . All the compounds are found to be semiconductors, except Pb 2 Bi 2 Te 5 in stacking A, which shows typical metallic character with overlapping between conduction and valence bands. The valence band maximum (VBM) and conduction band minimum (CBM) for PbTe, Bi 2 Te 3 and Pb 2 Bi 2 Te 5 stacking B locate at the Γ point with non-SOC calculations. For Bi 2 Te 3 and stacking B the location of both the VBM and CBM differs with and without SOC. Both structures bear a direct and indirect gap without and with SOC, respectively. For Bi 2 Te 3 with SOC the VBM locates near the Z point and the CBM is along the Z-F direction with relatively flat bands. For all the compounds, one can observe that SOC splits both high-lying conduction band degenerated orbitals and low-lying valence band ones at high symmetry points. For Pb 2 Bi 2 Te 5 , SOC affects more markedly the bottom conduction orbitals than the top valence ones. Therefore the changes in thermoelectric properties should be more noticeable for the n-type doped compound. In addition without SOC the bottom conduction orbital exhibits a typical parabolic behaviour around the point whereas with SOC the band is flat and extends itself between and A. In the latter case, the band effective mass and electron mobility should then be lighter and higher, respectively, both leading to lower Seebeck coefficient and to higher electrical conductivity. From these effects, coupled with the fact that the gap is significantly lowered by SOC, and assuming that the thermal conductivity is not affected, the zT value for n-type Pb 2 Bi 2 Te 5 should be degraded when considering the SOC effect in the calculations and all the more that the temperature increases. These assumptions have been verified by additional calculations (see Fig. S5 in supplementary materials).

Since the effect of SOC on the band structures is significant for all the compounds, except in very few cases where notified, we have included this effect in the following electronic calculations. Both the valence and conduction bands of the two Pb 2 Bi 2 Te 5 stacking sequences have mixed contributions from Pb-Te and Bi-Te states belonging from different electrons states as shown in Figure 4. The DOS of the stacking B is more closely related to that of PbTe than that of Bi 2 Te 3 with a strong Pb 6s state below the Fermi level, while the DOS of stacking A is more similar to that of Bi 2 Te 3 . The strong peak near -0.4 eV comes from the hybridization between Pb-6s, Pb-6p, Bi-6s and Te-5p for stacking A, and from a strong hybridization of Pb-6s, Bi-6s, Bi-6p and Te-5p for stacking B. In the stacking A where Pb atoms locate near the slab edge, an additional peak occurs at about -0.23 eV from a strong hybridization of Bi-6p, Pb-6p and Te-5p, especially with Te3 located at the edge of the ennead atoms slab. These results show that the inner atoms of the slab are the main contributors to the electronic band structure, the contribution of the atoms located at the edge of the slab being larger around the Fermi level. For stacking B, the tellurium atom at the slab edge (Te3) contributes less to the states near the top of the valence band due to the large electronegativity of the Bi atom to which it is bonded.

Electronic transport properties

The electronic conductivity (σ ) and the electronic thermal conductivity (κ e ) calculated by BolzTrap2 are given as the ratios σ /τ and κ/τ, respectively. The relaxation time τ depends on the charge carrier concentration, the temperature and the electrons energy. However, treating τ as a constant is not a satisfactory approximation to get thermoelectric properties. According to Ref 28 , τ can be evaluated from effective mass approximation with deformation potential (DP) theory. The mean free path , the mobility µ and the relaxation time τ can be defined as:

1 = m 2 e E i d 2 k 0 T π h4 c ii µ = 2(2π) 1/2 eh 4 c ii 3m 5/2 e (k 0 T ) 3/2 E i d 2 τ = µm e
where k 0 is the Boltzmann constant, c ii is the lattice elastic constant, and E i d is the deformation potential expressed as:

E i d = ∂ E ∂ a/a 0
where a 0 is the equilibrium lattice constant.

The effective mass m * is derived from the second derivative of E along the a-axis and the c-axis direction:

m * k,l = h2 /(∂ 2 E/∂ k k k l ).
For isotropic PbTe, only the c-axis has been considered while for the anisotropic Bi 2 Te 3 and Pb 2 Bi 2 Te 5 , both in-layers a-axis and cross-layers c-axis have been considered. The DP has been calculated by using seven values of lattice constant, namely the equilibrium one and six around this value chosen by step of 0.01 Å. A deep core state has been used as a reference to evaluate the absolute DP. The results are listed in Table 3.

The elastic constants of bulk PbTe, Bi 2 Te 3 and Pb 2 Bi 2 Te 5 have been calculated by fitting the energy-strain curve to quadratic or sextuple polynomial, as detailed in Ref 45 . We have implemented a series of distortions on the equilibrium structure. The bulk modulus has been obtained by fitting the results with the Birch-Murnaghan equation. Based on the space group, the elastic constants of PbTe (C 4.

The temperature and doping dependence of the electronic properties S, σ and κ e are shown in Fig. 5 and in Fig. 6. Since the best thermoelectric properties are obtained for semiconductors with carriers concentrations around 10 19 carriers/cm 3 , 50 the doping dependence was investigated in the domain of carriers concentration ranging from 10 17 /cm 3 to 10 22 /cm 3 . The calculations have been performed for both electron-doping and hole-doping in in-layers a-axis (xx) and cross-layers c-axis (zz). As shown in Fig. 5, the maximum | S | of stacking A appears in the area of 20 carriers per cubic centimeter) in a wide temperature range. Irrespective of the doping type the Seebeck coefficient absolute values are of the same order of magnitude and the anisotropy is small. By contrast to S, the in-layer σ is more than an order of magnitude larger than that in cross-layer, which is expected from the structure since the separated layer slabs along the c-axis hold together by vdW interactions. As to the stacking B in Fig. 6, the optimum doping level is relatively lower than for stacking A (around 10 18 -10 19 e/cm 3 ). Generally, stacking B has a higher S and a lower σ than stacking A, due to the metallic character of stacking A.

Lattice thermal conductivity

After having comprehensively determined the electronic part of TE properties we have accurately calculated the lattice thermal conductivity κ l of Pb 2 Bi 2 Te In a previous work [START_REF] Hellman | [END_REF] it has been reported that, van der Waals functionals yield underestimated bulk modulus by about 20-40% compared to LDA, resulting in too low phonon frequencies that do not agree with experimental data, hence these functionals have been discarded for our study. Furthermore, the harmonic phonon spectrum mostly depends on the choice of the lattice parameters and little on the choice of the functional 52 . Hence, the LDA functionals have been used to determine the harmonic and anharmonic force constants.

The longitudinal optical (LO) -transverse optical (TO) splitting at the Gamma point, which has a significant effect on the phonon dispersion of PbTe caused by large Born effective charges, has been evaluated by adding the nonanalytical term to the dynamical matrix. The spectra of PbTe and Bi 2 Te 3 are shown with and without LO-TO splitting in Fig. S6. A good agreement is obtained between the phonon dispersion curves of PbTe and Bi 2 Te 3 along a high symmetry path and experimental data from Ref 53,54 (Fig. 7). The Born effective charges calculated by the Berry phase method 55 and dielectric constants are listed in Table 5. Because of the rock salt structure, the Born effective charge tensor of PbTe has identical diagonal elements. Furthermore, PbTe being a binary compound we have the relation

Z * Pb = -Z * Te .
The TO modes are very soft, which is well in agreement with both experimental data 54 and previous DFT 8). Indeed according to Mal et al. 22 , the peak at 131 cm -1 is related to the interactions at the boundary of the layers and in the present work the atomic contributions to the 4 THz frequencies for Pb 2 Bi 2 Te 5 stacking B solely belong to the Bi and Te3 atoms.

To evaluate the lattice thermal conductivity κ l , the Boltzmann transport equation (BTE) [START_REF] Kittel | Introduction to Solid State Physics[END_REF] has to be solved self consistently including all the scattering processes. The compounds of interest are considered as ideal ones without scattering impurities. The κ l can be calculated using a full solution of the linearized phonon Boltzmann equation (LBTE) following the method introduced in Ref [START_REF] Chaput | [END_REF] . The LBTE gives a rigorous way to evaluate lattice thermal conductivity by considering phonon-phonon interactions but necessitates huge calculations. In previous researches on metal and semiconductor 64,65 the single-mode relaxation time approximation (RTA) 66 has been proved to be a fast and reliable choice to simplify the calculations. Within the RTA theory, the lattice thermal conductivity is given by:

κ l = 1 NV 0 ∑ λ C λ υ λ ⊗ υ λ τ λ
where V 0 is the unit cell volume, υ λ is the group velocity and τ λ is the phonon scattering time in the specific phonon mode (λ ). The summation runs over all the band in the Brillouin zone. The specific mode heat capacity of the (λ ) phonon with frequency ω λ is C λ = k B (hω λ β ) 2 n 0 λ (n 0 λ + 1), with β = 1/(k B T ). n 0 λ is the Bose-Einstein distribution function.

We have calculated the lattice thermal conductivity using both LBTE and RTA at 300 K, 500 K, 700 K and 900 K for the two stacking sequences. The results of in-layer a-axis, cross-layer c-axis, and average thermal conductivity κ ave = (2 × κ xx + κ zz )/3 67 are summarized in Table 6. With LBTE, the thermal conductivity at 300 K amounts to κ xx = 0.98 Wm -1 K -1 in the in-layer direction and κ zz = 0.35 Wm -1 K -1 in the cross-layer direction for stacking B. Pb 2 Bi 2 Te 5 with ordered sequence A has a similar κ xx value to that of stacking B, while the cross-layer one κ zz (0.19 Wm -1 K -1 ) is noticeably smaller. Both structures are found to have very low and anisotropic lattice thermal conductivity. Comparing LBTE and RTA results, irrespective of the stacking, one can note that, the cross-layer κ zz values are similar while the in-layer ones are different. When increasing temperature, κ l decreases in the two directions.

The cumulative lattice thermal conductivity, phonon lifetime τ, group velocity in a-axis and c-axis directions for Pb 2 Bi 2 Te 5 with the two stacking sequences, calculated at 300 K by using RTA, are shown in Fig. 9. The temperature dependence of the thermal conductivity of the two stacking sequences in a-axis and c-axis directions as well as that of the average thermal conductivity are presented in Fig. 10. It can be clearly seen that the in-layers thermal conductivities of the two stacking sequences are close, while the cross-layers thermal conductivities of stacking B is almost twice as large as that of stacking A. The interaction between the acoustic branches and the optical ones plays an important role in the lattice thermal conductivities. Based on the cumulative κ l curves, we can divide the frequency domain into three parts: (a) 0-1.8 THz, (b) 1.8-3.4 THz and (c) 3.4-4.1 THz. In the first part (a), the cumulative κ l curves of both structures show a S-shape like behavior. In this frequency domain the cumulative lattice thermal conductivities of the two stacking sequences are similar. In the second part (b), correlated with the phonon band gap observed between 2.3 to 2.5 THz in stacking A (Fig. 8), the cumulative κ l of this stacking sequence reaches a plateau. In the third part (c) the optical phonons come into play and participate, though weakly, to the heat conduction. The contribution provided by the first, second and third parts to heat conduction amounts to 78%, 10%, 12% for stacking A and 66%, 28%, 6% for stacking B, respectively. Meanwhile, the gap in the optical branches in stacking A leads to fewer phonon scattering, producing much longer phonon lifetimes (Fig. 9b).

Thermoelectric figure of merit

From the Seebeck coefficients and the electronic and thermal conductivities the figure of merit zT can be evaluated as a function of carriers concentration and temperature. The results are depicted in Fig. 11,12. The zT values of Pb 2 Bi 2 Te 5 with stacking A depend little on the considered direction and carriers types. By contrast, for stacking B the zT values are much higher for p-type carriers that for n-type ones. The conditions to obtain the maximum zT are 200-300 K and 1 × 10 18 -5 × 10 19 h/cm 3 along the a-axis direction, and 600-1000 K and 1 × 10 20 -5 × 10 21 h/cm 3 along the c-axis one. The highest value in the a-axis direction is 4.02 and that in the c-axis direction is 2.26, which are higher than those already reported in literature for other layered chalcogenide compounds (e.g. 1.3 for the Sb 2 Te 3 -GeTe compounds 68 , 1.4 for the Bi 2 Te 3 -GeTe ones 69 and 1.86 for Bi 0.5 Sb 1.5 Te 3 70 ). Despite the lower lattice thermal conductivity of stacking A, its electronic thermal conductivity is much higher than that of stacking B as a result of the metallicity (Figure 5 and6).

Summary of the electronic properties for both stackings of Pb 2 Bi 2 Te 5

We have gathered in Table 7 the electronic properties of stacking A and stacking B of Pb 2 Bi 2 Te 5 at low and high temperatures of 300 K and 800 K with p-and n-type doping of 10 19 cm -3 and 10 20 cm -3 , respectively. All the properties are reported for the xx and zz directions to highlight the anisotropy of the compounds.

The best zT values are obtained for the stacking B of Pb 2 Bi 2 Te 5 for p-type doping. This is correlated with a correspondingly high Seebeck coefficient and moderate electronic conductivity. Although the moderate electronic conductivity may be detrimental to the power factor S 2 σ appearing in the zT expression, this implies a low electronic thermal conductivity. Combined with a fairly low lattice thermal conductivity, especially in the zz direction, the figure of merit of stacking B Pb 2 Bi 2 Te 5 is improved compared to that of stacking A. Regarding stacking A, its metallic character leads to a too low Seebeck coefficient and too high electronic thermal conductivity that are not compensated by the indeed very low lattice thermal conductivity.

Due to both a thermal conductivity lying in the range of what is generally observed for glassy compounds and an electronic conductivity and a Seebeck coefficient typical of moderately doped semi-conductors, we can say that Pb 2 Bi 2 Te 5 stacking B meets the "phonon-glass, electron-crystal" (PGEC) concept.

Mechanical, dynamical and relative stabilities of stacking A and stacking B Pb 2 Bi 2 Te 5

The elastic constants that have been calculated in section 4.2.2 for permitting the determination of τ are also useful results (Table 4) to evaluate the mechanical stability of structures. The necessary and sufficient conditions of stability for the rhombohedral I system are given as 71 : Table 7 Electronic properties for Pb 2 Bi 2 Te 5 , stacking A and stacking B at 300 K and at a carrier concentration of 10 19 cm -3 (for both n and p). Electronic conductivity σ × 10 -5 in Ω -1 m -1 , electronic thermal conductivity in WK -1 m -1 , Seebeck coefficient in

C 11 >| C 12 |; C 44 > 0
µVK -1 , lattice thermal conductivity in WK -1 m -1 .
Electronic properties at 300 K and at a carrier concentration of 10

C 2 13 < 1 2 C 33 (C 11 +C 12 ) C 2 14 < 1 2 C 44 (C 11 -C 12 ) = C 44 C 66
whereas for cubic one they are:

C 11 -C 12 > 0; C 11 +C 12 > 0; C 44 > 0
In view of these mechanical stability conditions, both the two Pb 2 Bi 2 Te 5 stacking sequences are mechanically stable. This result is supported by the negative value of the cohesive energy at zero Kelvin for both stacking. The cohesion energies have been determined with the general formula:

E coh = E T -∑ i E i
where E T is the total energy of the compound of interest, and E i is the energy of the constitutive atom. The cohesive energies are -2.49 eV/at. for stacking A and -2.54 eV/at. for stacking B.

The phonon dispersion curves (section 4.2.3) provide a hint towards the dynamical stability of structures. The positive phonon values (Figure 8) in the whole Brillouin zone show that there is no dynamical instability for both stacking A and stacking B.

In agreement with literature results (see Introduction), we have found that both stacking A and stacking B can exist, however under given temperature and pressure only one of these should be stable, the other one being metastable. To investigate the relative stability of the two stacking, the Helmoltz free energy has been calculated. According to Wang et al. 72 , the free energy of a crystalline material F(T ), when neglecting the configurational contribution, can be expressed as follows :

F(T ) = E 0K + F vib + F el
Where E 0K is the energy at zero Kelvin, F vib is the vibrational free energy and F el is the electronic free energy. The vibrational free energy reads [START_REF] Dove | Introduction to lattice dynamics[END_REF] :

F vib = 1 2 ∑ k,ν hω(k, ν) + k B T ∑ k,ν ln 1 -exp hω(k, ν) k B T
and the electronic free energy is defined as 74

F el = E el -T S el with E el = n(ε) f εdε + n(ε)εdε and S el = -k B n(ε) [ f ln( f ) + (1 -f ) ln(1 -f )]
T , k B , h, f and ε are the temperature, Boltzmann constant, reduced Planck constant, Fermi-Dirac function and electrons energies, respectively. ω(k, ν) is the phonon frequency of the k th band and ν th vector. The Helmoltz free energy of both Pb 2 Bi 2 Te 5 stacking is depicted in Figure 13 as a function of the temperature. Negative free energy is found for both stacking from 0 K to 1000 K and the free energy of stacking B is found to be lower than that of stacking A, suggesting that stacking B is the stable form of Pb 2 Bi 2 Te 5 whereas stacking A is the metastable one.

Electron charge density and chemical bonding

Electron charge density value and distribution play an important role in bonding trend, which will affect the thermoelectric properties. Therefore, the analysis of electron charge density and chemical bonding developed by Bader 25 and following QTAIM methods 26,[START_REF] Otero-De-La-Roza | [END_REF] 8). The difference in scalar distribution between the two structures is illustrated in Fig. 14c. Based on the lattice parameters in Table 1, there is a small difference between atomic positions of the two stacking structures in c-axis direction. In order to compare the difference in charge density between the two stackings, we meshed the (110)-planes with the same grid in the direction [-1, 1, 0], and the same grid between each pair of adjacent atoms in the direction [0, 0, 1]. Zhang 76 found a strong positive correlation between the anisotropy in lattice thermal conductivity and the anisotropy in electron charge density. As shown in Fig. 14c, the charge difference mainly locates in the c-axis direction, which accords with the anisotropy in lattice thermal conductivity reported in Table 7.

As mentioned above and reported in Ref 77 , bonds properties play an important role in thermoelectric properties. Therefore, the topological properties of bond critical points (BCPs) in Fig. 14a,b are discussed, the results of which being tabulated in Table 8, on the basis of bond path length and angle, the charge density ρ and its Laplacian, the Table 8 Topological properties of Pb 2 Bi 2 Te 5 corresponding to the BCPs in Fig. 14 a,b; r1 is the distance between BCP and Bi or Pb atoms, or Te atoms for b5; r2 is the distance between BCP and Te; r1 and r2 are in unit of pm; the angle is between BCPs and their two ends; ρ and ∇ 2 ρ are the electron charge density and its Laplacian at the BCPs in unit of 10 -2 e/bohr 3 and 10 -2 e/bohr 5 , respectively; G, V and H represent the kinetic, potential and total energy densities at the BCPs in unit of 10 -2 a.u./bohr bond degree (H/ρ), the kinetic energy per electron (G/ρ) and the dimensionless | V | /G ratio. The electron charge density in the expected interlayer van der Waals Te-Te (b5) bondings of stacking A is much higher than that of stacking B. This lower charge density in stacking B results in a longer interatomic distance that can be explained by the higher electronegativity of Bi, which weakens the Te-Te interaction, compared to that of Pb. While within the ennead atoms slab of stacking A two poor electron density areas are found at A-b3 and A-b5 positions in Fig. 14c leading to smaller lattice thermal conductivity in c-axis direction, only one such area is found in stacking B at B-b5.

To explore the reason for the previous observations, we should catch a glimpse of bonds properties by identifying the BCPs. Usually, closed shell interactions (ionic, H-bonds and vdW) have a large positive value of ∇ 2 ρ over the entire interaction region, G/ρ > 1, | V | /G < 1 and a small ρ. Conversely, ∇ 2 ρ < 0, G/ρ < 1, | V | /G > 2 and a large ρ are expected for shared interactions (covalent or polar bonds). 78 In our study, all BCPs show small value of ρ and positive ∇ 2 ρ with 1 <| V | /G < 2, which cannot be considered as pure covalent nor pure closed-shell bonds. This region has been identified by Espinosa and called a transit region 79 . In order to clearly show the tendency of BCPs, as proposed in Ref 80 According to Espinoza's bonding classification, these weak Te-Te bonds belong to the transit region and not to the closed shell one as would be expected for van der Waals interactions. This could be explained by the fact that the transit region limits could be more or less displaced towards the shared shell or the closed shell regions, according to the nature of the interacting atomic pair 77 . The electronegativity values of Pb, Bi, and Te in the Pauling scale are 1.87, 2.02 and 2.1, respectively. Hence one can assume that the ionic-like interaction and the polarity of Pb-Te are much stronger than those of Bi-Te. The BCPs properties for the binary compounds PbTe and Bi 2 Te 3 are gathered in Table S1. The results show one kind of Pb-Te bond and 3 kinds of Bi-Te ones. In stacking B, the Pb-Te bonds are at the center of the slab with an approximate periodic sequence (-Te-Pb-Te-Pb-Te-) (see Fig. 1d) leading to two kinds of Pb-Te bonds and as shown in Table 8 and in Fig. 15, the properties of the two central Pb-Te BCPs (B-b1 and B-b2) are close to each other. By contrast, the properties of the two Bi-Te BCPs (B-b3 and B-b4) are different. Symmetrically in stacking A, the Bi-Te bonds are in the center of the slab but they show different central BCPs (A-b1 and A-b2) properties. The Pb-Te BCPs (A-b3 and A-b4) also differ from each other. These differences could be related to the difference in the Te-Te interaction strength, which is larger when Pb is located at the slab edge (stacking A) due to a lower electronegativity for Pb in comparison to Bi as already mentioned above. The Te atom being less electron-rich when bonded to Bi in stacking B, this induces a weaker Te-Te bond across the slabs gap, which is correlated by the larger bond length and lower absolute bond degree (Table 8). The bond distance between the Te atoms in stacking B structure (370 pm) is much smaller than the sum of the Te van der Waals radii (412 pm) suggesting that these bonds, although being weak, are not strictly of van der Waals type. This is in agreement with Gatti's 78 and Yang's 77 works that assign positive bond degree to van der Waals interactions. This difference in the Te-Te interaction strength can also explain the differences in bulk modulus values calculated above: 47.3 GPa for stacking A and 45.5 GPa for stacking B (see §4.2.2 Electronic transport properties). A lower symmetry in BCPs properties combined with stronger Te-Te interactions in stacking A should lead to a stronger anisotropy and a larger electronic conductivity in this structure. The calculated electronic properties reported in Table 7 confirm this assumption.

As mentioned above, the non-equivalent Pb-Te bonds are dissimilar in stacking A as the electron density differs from one another. One can assume that this difference originates from the position of the Pb-Te bond with respect to that of the Te-Te one. In order to check this assumption, we have built a new stacking-A-sequence structure with an increased number of layers in the slabs (21 instead of 9) and calculated, after optimisation, the topological properties of their bonds. The values are gathered in table S2. Since the structures are symmetrical, for sake of clarity we only give the properties of the non-equivalent bonds. The charge densities are plotted in Fig. 16. As expected the differences between the charge density at the critical points decrease when one move away from the Te-Te bond. In stacking B, the non-equivalent Pb-Te bonds are similar and this is also true when the number of layers in the slab is increased. However if we replace the Bi atom at the edge of the slab by a Pb one the Pb-Te bonds differentiate from one another and the difference decreases when one move away from the Te-Te bond. Therefore it seems that the Bi atom at the edge of the slab in stacking B isolates the Pb-Te bonds from the influence of the Te-Te ones between the slabs.

Conclusion

In this work, we have investigated the stability, the thermoelectric properties and the bonding ones of Pb 2 Bi 2 Te 5 for two stacking sequences. In agreement with experimental findings, for both stacking A and B no instability has been evidenced. In the range of temperatures 0 K-1000 K stacking B seems to be the stable form of Pb 2 Bi 2 Te 5 , whereas stacking A is the metastable one. The two stacking sequences show remarkable differences in electronic properties. Pb 2 Bi 2 Te 5 with stacking B is a narrow energy gap semiconductor, with calculated indirect gap of 0.116 eV, while stacking A shows metallic character without energy gap. The thermoelectric properties of the two stacking sequences, namely the Seebeck coefficient, electronic conductivity and thermal conductivity have been calculated. To go beyond the relaxation time approximation, τ has been evaluated from effective mass approximation with deformation potential (DP) theory. The calculated results show that Pb 2 Bi 2 Te 5 with stacking B has higher Seebeck coefficient, lower electronic conductivity and lower electronic thermal conductivity. Based on the calculation of phonons spectra, stacking A has stronger anisotropy and lower lattice thermal conductivity than stacking B. As a consequence, the zT values of stacking B are much higher than those of stacking A. The highest value in the a-axis direction is 4.02 and that in the c-axis direction is 2.26 with p-type carriers. The conditions of high zT are 200-300 K and 1 × 10 18 -5 × 10 19 h/cm 3 in a-axis direction, and 600-1000 K and 1 × 10 20 -5 × 10 21 h/cm 3 in the c-axis one.

Quantum theory of atoms in molecules (QTAIM) calculations have been performed to explore the reasons for the properties differences. The stacking sequence does not change the number or type of CPs, however, the bond Te-Te between two slabs and the Pb-Te bond near the slab edge (stacking A) are strengthened when changing the Pb atom location from the center of the slab (stacking B) to the edge (stacking A). The changing of sequence from stacking B to stacking A breaks the symmetry of the bonds strength and leads to various BCPs properties in the structure, which yield to both higher anisotropy and electronic conductivity.

Figure 1

 1 Figure 1 (a) Primitive cell of PbTe; (b) primitive and conventional cells of Bi 2 Te 3 ; (c) 1 × 1 × 2 supercell of Pb 2 Bi 2 Te 5 in stacking A; (d) 1 × 1 × 2 supercell of Pb 2 Bi 2 Te 5 in stacking B.

Figure 2

 2 Figure 2 Calculated electronic band structures of PbTe (a), Bi 2 Te 3 (b) with PBE-GGA functional. The blue lines and the red lines represent the band structures without and with SOC, respectively.

Figure 3

 3 Figure 3 Calculated electronic band structures of Pb 2 Bi 2 Te 5 stacking A (a) and Pb 2 Bi 2 Te 5 stacking B (b) with PBE-GGA functional. The blue lines and the red lines represent the band structures without and with SOC, respectively.

Figure 4

 4 Figure 4 Total DOS (a), projected DOS of PbTe and Bi 2 Te 3 (b) and projected DOS of Pb 2 Bi 2 Te 5 stacking A and B (c,d) calculated with PBE-GGA functional.

5 .

 5 In order to use QUANTUM ESPRESSO for this purpose we have first re-optimized the structures of PbTe, Bi 2 Te 3 and Pb 2 Bi 2 Te 5 stacking A and B with this package using the PBE associated with vdW-DF2 exchange-correlation functional without SOC. The obtained equilibrium lattice constants are a = 657 pm for PbTe, a = 435 pm, c = 2985 pm for Bi 2 Te 3 , a = 441 pm, c = 1752 pm for stacking A and a = 442 pm, c = 1732 pm for stacking B, which are quite close to the constants optimized by WIEN2k and the experimental ones.

Figure 5

 5 Figure 5 Calculated Seebeck coefficient, electronic conductivity and electronic thermal conductivity versus temperature and n-and p-type doping levels in in-layer a-axis (xx) and cross-layer c-axis (zz) direction for Pb 2 Bi 2 Te 5 stacking A.

Figure 6

 6 Figure 6 Calculated Seebeck coefficient, electronic conductivity and electronic thermal conductivity versus temperature and n-and p-type doping levels in in-layer a-axis (xx) and cross-layer c-axis (zz) direction for Pb 2 Bi 2 Te 5 stacking B.

Figure 7

 7 Figure 7 Phonon dispersion relations and total and projected phonon density of states for PbTe (a) and Bi 2 Te 3 (b) in a unit cell at equilibrium. The solid red lines are from this work, and the blue dots are experimental results extracted from Ref 54 for PbTe and Ref 53 for Bi 2 Te 3 .

Figure 8

 8 Figure 8 Phonon dispersion relations of Pb 2 Bi 2 Te 5 for stacking A (a) and stacking B (b), as well as total and projected phonon density of states.

Figure 9

 9 Figure 9 Calculated cumulative lattice thermal conductivity for Pb 2 Bi 2 Te 5 stacking A (blue) and stacking B (orange) with RTA (a). Phonon lifetimes for A and B (b). Group velocity for A (c) and B (d) in the in-layer a-axis and cross-layers c-axis. Results obtained at 300 K.

Figure 10

 10 Figure 10 Calculated lattice thermal conductivity κ of Pb 2 Bi 2 Te 5 stacking A (a) and stacking B (b) vs. temperature. The evolution is plotted for the in-layer a-axis, cross-layers c-axis and average lattice thermal conductivity.

Figure 11

 11 Figure 11 Calculated figure of merit zT of Pb 2 Bi 2 Te 5 stacking A in a-axis(xx) and c-axis(zz) directions vs. temperature and carriers concentration.

Figure 12

 12 Figure 12 Calculated figure of merit zT of Pb 2 Bi 2 Te 5 stacking B in a-axis(xx) and c-axis(zz) directions vs. temperature and carriers concentration.

Figure 13

 13 Figure 13 Helmholtz free energy per unit cell versus temperature for stacking A and stacking B of Pb 2 Bi 2 Te 5 .

  will help us understand the reason of remarkable properties difference between the two stacking sequences. The electron charge density together with the flux lines generated by ∇ρ of the two stacking sequences in (110)-plane are shown in Fig. 14a,b with all critical points (CPs) namely nuclei, bonds (b), rings (r), and cages (c). In the studied structures, all the bond critical points (BCP) of the Pb-Te and Bi-Te bonds lie closer to the lead and bismuth atoms than to the Te ones. The change of stacking sequence does not change the numbers and types of critical points. Moreover, there is no charge accumulation along the bonds in the whole plane as evidenced by the positive Laplacian values at the BCP (Table

Figure 14

 14 Figure 14 CPs and flux lines derived from ∇ρ(r) (light green lines) of stacking A (a) and B (b) in (110)-plane, as well as the electron charge density difference ρ Aρ B (c) between the A and B stackings.

  the | V | /G versus H/ρ plot is illustrated in Fig. 15. The Te-Te bonds of Bi 2 Te 3 , Pb 2 Bi 2 Te 5 stacking A and Pb 2 Bi 2 Te 5 stacking B locate in the top left area of Fig. 15, with a small magnitude of H/ρ and | V | /G, the latter being larger than though close to 1. The bond degree and dimensionless | V | /G decrease along the sequence of Pb 2 Bi 2 Te 5 stacking A, Bi 2 Te 3 , Pb 2 Bi 2 Te 5 stacking B.

Figure 15

 15 Figure 15 Dimensionless | V | /G vs. bond degree H/ρ.

Figure 16

 16 Figure 16 Electron density at the bond critical points along the slab from the edge (Te-Te interaction on the left) to the middle of the slab (on the right) in stacking-A-sequence structures with 9 (S1-1) and 21 (S1-2) layers, in stacking-B-sequence structures with 9 (S2-1) and 21 (S2-2) layers, and in stacking-B-sequence structures with the Bi atom replaced by a Pb one with 9 (S3-1) and 21 (S3-2) layers. The dashed lines are just meant to guide the eye.

Table 1

 1 Calculated lattice constants in pm and atomic position of fractional coordinates compared with experimental values. The first atomic position is 0 for Te2 in Bi 2 Te 3 , Te1 in Pb 2 Bi 2 Te 5 for both stacking A and stacking B.

			LDA	PBE	rev-vdW-DF2 Other works
	PbTe	a	638.4	657.2	-	645.4 38
		a	436.0	445.7	435.0	437 32 , 438 39
	Bi 2 Te 3	c ν	2999.9 3074.7 3020.8 0.209 0.211 0.208	3034 32 , 3049 39 0.212
		µ	0.400	0.399	0.400	0.400
		a	443.0	454.4	451.6	446 18
	Pb 2 Bi 2 Te 5 stacking A	c Bi Te2 0.2110 0.2066 0.2088 1687.8 1733.4 1721.1 0.1128 0.1117 0.1125 Pb 0.3391 0.3394 0.3396	1752 18 0.109 0.218 0.329
		Te3 0.4315 0.4293 0.4301	0.440
		a	442.8	454.6	450.8	442 20 *
	Pb 2 Bi 2 Te 5 stacking B	c Pb Te2 0.2192 0.2194 0.2189 1732.0 1740.8 1768.7 0.1099 0.1097 0.1100 Bi 0.3326 0.3349 0.3323	1786 20 *
		Te3 0.4299 0.4294 0.4261
	* Lattice parameter from thin film.

Table 2

 2 Band gaps (eV) for PbTe, Bi 2 Te 3 , Pb 2 Bi 2 Te 5 stacking A and Pb 2 Bi 2 Te 5 stacking B calculated with various functionals compared to the experimental ones. Both the non-SOC and SOC results (in parenthesis) are presented.

		LDA	PBE	EV	TB-mBJ PBE/vdw HSE (this work) HSE (from lit.) Exp.
	PbTe	0.574	0.827	1.149	1.196	-	1.26	1.13 42	0.1941 43
		(0.084) (0.172) (0.358) (0.441)	-	(0.379)	(0.29) 42
	Bi 2 Te 3	0.168 (0.088) (0.113) (0.116) (0.126) 0.256 0.521 0.554	0.316 (0.114)	0.633 (0.306)	-(0.311) 44	0.1658 16
	Pb 2 Bi 2 Te 5 stacking A 0	0	0	0	0	0	-	-
		(0)	(0)	(0)	(0)	(0)	(0)
	Pb 2 Bi 2 Te 5 stacking B	0.380	0.504	0.766	0.785	0.489	0.887	-	-
		(0.145) (0.113) (0.072) (0.07)	(0.116)	(0.110)

  11 , C 12 , C 14 ), Bi 2 Te 3 (C 11 , C 12 , C 13 , C 14 , C 33 , C 44 ) and Pb 2 Bi 2 Te 5 (same as with Bi 2 Te 3 ) are listed in Table

Table 3

 3 Deformation potential constants Ed and carriers effective mass m * at zero Kelvin in the c direction and in the ab plane calculated with PBE-GGA functional.

	Structures Charge	c	m *	a,b	c	E d (eV) a,b
	PbTe	holes electrons -0.368 0.259			10.64 9.62
	Bi 2 Te 3	holes electrons -0.131 -0.019 0.123 0.018	10.73 14.05 9.74 13.92
	Pb 2 Bi 2 Te 5 holes	0.077		0.023	10.84 13.22
	stacking A electrons -0.305 -0.06	8.86	12.04
	Pb 2 Bi 2 Te 5 holes	1.12		0.106	9.75	12.51
	stacking B electrons -0.343 -0.132	9.89	12.69

intermediate doping level (range 10 19 -3×10

Table 4

 4 Bulk modulus B (GPa) and elastic constants C i j (GPa) of PbTe, Bi 2 Te 3 and Pb 2 Bi 2 Te 5 calculated with PBE-GGA together with calculated results and experimental data from the literature.

		B	C 11	C 12	C 13	C 14	C 33	C 44
	PbTe	37.5 99.1	6.7	-	-	-	13.4
	Other Calc. 46 40.4 111.0 6.7	-	-	-	14.4
	Exp. 47	39.8 105.3 7.0	-	-	-	13.2
	47.0 78.0 Other Calc. 48 36.4 82.8 Bi 2 Te 3	22.7 38.1 23.2 20.8 40.1 -	69.4 46.7 60.4 40.1
	Exp. 49	40.9 74.4	29.2 29.2 15.4	51.6 29.2
	Pb 2 Bi 2 Te 5 stacking A	47.3 77.7	27.1 38.5 -21.1 62.4 44.1
	Pb 2 Bi 2 Te 5 stacking B	45.6 77.0	21.0 37.3 -22.4 65.4 45.4
	works 46,54 .

Table 5

 5 Calculated Born effective charges Z * (e) and dielectric constants ε ∞ in in-layer a-axis(xx) and in cross-layer c-axis(zz).

			This work [xx, zz] Other works [xx, zz]
	PbTe	Z * Pb ε ∞xx 37.61 [5.75, 5.75]	[6.5, 6.5] 56 [5.8, 5.8] 46 36.9(77 K) 32.8(300 K) 57
	Bi 2 Te 3	Z * Bi Z * Te1 Z * Te2 ε ∞xx 43.96 [6.83, 6.83] [-6.79, -6.79] [-3.46, -3.46]	[7.29, 3.58] 58 [-6.03, -5.25] 58 [-4.28, -0.96] 58 45-50 59-61
		ε ∞zz	42.36	38 60
	Pb 2 Bi 2 Te 5 stacking B	Z * Bi Z * Pb Z * Te1 Z * Te2 Z * Te3	[7.86, 3.93] [5.84, 5.98] [-4.41, -1.27] [-6.29, -5.71] [6.69, -5.93]
		ε ∞xx 54.57
		ε ∞zz	42.77

Table 6

 6 Lattice thermal conductivity (W/(m•K)) for Pb 2 Bi 2 Te 5 A and B stacking sequences obtained from LBTE and RTA at the temperatures of 300 K, 500 K, 700 K and 900 K. The average thermal conductivity values as well as that in the in-layer xx and cross-layer zz direction are given

				LBTE			RTA
			κ xx	κ zz	κ ave	κ xx	κ zz	κ ave
		300 K	0.98	0.19	0.70	0.65	0.18	0.50
	Stacking A	500 K 700 K	0.59 0.42	0.11 0.08	0.43 0.31	0.39 0.28	0.11 0.08	0.30 0.21
		900 K	0.33	0.06	0.24	0.22	0.06	0.17
		300 K	0.98	0.35	0.75	0.91	0.34	0.72
	Stacking B	500 K 700 K	0.59 0.42	0.21 0.15	0.46 0.33	0.55 0.39	0.20 0.14	0.43 0.31
		900 K	0.32	0.12	0.25	0.30	0.11	0.24

  19 cm -3

			σ xx	σ zz	κ e xx	κ e z	S xx	S zz	κ l xx	κ l zz	zT
	stacking A	n-type p-type	106 393	5.73 40.1	216 715	12.2 56.6	-31.7 49.9	54.2 108	0.65	0.65	0.0001 0.08
	stacking B	n-type p-type	6.71 11.4	0.202 0.154	4.98 7.46	0.28 0.10	-161 296	-171 312	0.91	0.34	0.88 3.63
			Electronic properties at 800 K and at a carrier concentration of 10 20 cm -3
			σ xx	σ zz	κ e xx	κ e z	S xx	S zz	κ l xx	κ l zz	zT
	stacking A	n-type p-type	86.2 315	5.37 32.9	414 1429	30.2 146	-24.6 56.8	2.44 87.5	0.24	0.07	0.004 0.08
	stacking B	n-type p-type	9.69 13.8	0.729 0.274	34.75 53.9	5.68 0.90	-138 206	-175 298	0.91	0.34	0.48 1.15

  3 .

	BCP	r1	r2	r1/r2	angle	ρ	∇ 2 ρ	G	V	H	| V | /G	H/ρ
	A-b1	158.0	165.8	0.953	179.54	3.064	3.068	1.37	-1.98	-0.61	1.44	-1.98
	A-b2	150.7	158.3	0.952	179.95	4.024	2.896	1.84	-2.95	-1.12	1.61	-2.77
	A-b3	168.3	176.6	0.953	179.06	1.985	2.730	0.87	-1.06	-0.19	1.21	-0.96
	A-b4	146.9	157.1	0.935	179.53	4.151	3.175	1.96	-3.12	-1.16	1.59	-2.80
	A-b5	177.7	177.7	1.000	180.00	1.773	3.152	0.75	-0.90	-0.14	1.18	-0.80
	B-b1	156.4	167.8	0.932	179.81	2.861	3.293	1.32	-1.81	-0.49	1.37	-1.72
	B-b2	156.1	167.3	0.933	179.89	2.898	3.316	1.34	-1.85	-0.51	1.38	-1.76
	B-b3	159.9	168.6	0.948	179.81	2.834	2.965	1.25	-1.76	-0.51	1.41	-1.80
	B-b4	150.6	158.6	0.949	179.89	4.075	2.592	1.82	-2.99	-1.17	1.64	-2.87
	B-b5	185.2	185.2	1.000	180.00	1.379	2.075	0.57	-0.63	-0.055	1.09	-0.40

Acknowledgements

The PhD thesis of Mr. W. Ma is financially supported by the China Scholarship Council (CSC). This work was granted access to the HPC resources of the Centre Informatique National de l'Enseignement Supérieur (CINES), Montpellier, France under allocation A0070806881 made by the Grand Equipement National de Calcul Intensif (GENCI). It was also granted access to the HPC resources of Aix-Marseille University financed by the project Equip@Meso (ANR-10-EQPX-29-01) of the program "Investissements d'Avenir" supervised by the Agence Nationale de la Recherche.

Notes and references


Hailong Yang

Marie-Christine Record
email: m-c.record@univ-amu.fr

P Boulet

Thermoelectric properties of Cu-Sb system compounds from DFT calculations

Keywords: DFT, QTAIM, transport properties, chalcogenides, thermoelectrics, structure-properties relationships

In this work the thermoelectric properties of Cu-Sb system compounds have been investigated to understand their relation with the electronic structures. The approach used is a theoretical one that combines the density functional theory, the Boltzmann transport theory and the quantum theory of atoms in molecules. The results show that the Seebeck effect in metallic systems originates from the existence of various types of chemical bonds. However this feature is not sufficient enough to generate effective thermoelectric effect. As in semi-conducting compounds, the intrinsic driving force for thermoelectric effect remains the variation of the density of states (DOS) near the Fermi level.

INTRODUCTION

As the thermoelectric (TE) properties are electronic ones, the analysis of the electronic interactions and as a consequence of bonds in materials should allow one to better understand the origin of the TE properties. For example, it has been shown in the compounds of the Cu-Sb-S system that weak chemical bondings in materials are utterly important to get a low lattice thermal conductivity. 1 In the late 1960s several models of the chemical bond, namely the ionicity, 2 partially ionic binding, 3,4 covalent bonding 5 and dielectric constants 6 ones, have been developed to investigate the properties of materials. However, due to the quantum nature of the chemical bonds, working out a general model that correlates the physical properties to the bond ones is made arduous. Nonetheless, it is widely admitted that the charge density is amenable to provide information on the bonding interaction. As it turns out Bader's quantum theory of atoms in molecules (QTAIM), [7][8][9][10] a topological tool applied to the electron density, is among the most advanced theories that allows for revealing bonding interaction nature in molecules and crystals. As with many other techniques, the crucial starting point is the partitioning of the whole system into subsystems. In QTAIM the molecule or the crystal is divided into the so-called atomic basins that are volumes of space containing a single atom delimited by a boundary where the flux of the electron density vanishes. From then on, the analysis of the electron density ρ(r) can be achieved through that of its gradient vector field ∇ρ(r) and the topological analysis of the electron density boils down to the search for critical points (CP) that are points in ρ(r) where ∇ρ(r) zeroes out. The electron density of stable crystals exhibits four types of CP among which the bond critical point (BCP) that characterizes, together with the bond path, the bonding interaction between two atoms. The bond path is the line joining two atoms where the electron density is maximal with respect to any neighbouring path. The BCP is characterized by the signature (3,-1), where 3 corresponds to the dimensionality and -1 is the sum of the signs of the electron density Laplacian ∇ 2 ρ. Within this tool, the chemical interactions have been classified into two types, namely the shared shell interactions for which ∇ 2 ρ is negative (incoming electron density) and the closed shell ones for which ∇ 2 ρ is positive (outgoing electron density). 11 Nonetheless, in many situations, in particular in semiconductors, this basic classification is insufficient and the electron density Laplacian information can be advantageously supplemented by both the local virial theorem that reads:

ℏ ! "# ∇ $ 𝜌(𝑟) = 2𝐺(𝑟) + 𝑉(𝑟) (1)
where G(r) and V(r) are the kinetic and potential energy densities, respectively, and the total energy density H(r) relation:

𝐺(𝑟) + 𝑉(𝑟) = 𝐻(𝑟) (2)
G(r) can be estimated using the approximation derived by Kirzhnits 12,[START_REF] Kirzhnits | Field theoretical methods in many-body systems[END_REF] and subsequently refined by Abramov 14 that relates ρ(r) to G(r):

𝐺(𝑟) = % &' (3𝜋 $ ) ! " 𝜌(𝑟) # " + & ($ [∇+(-)] ! 0(1) + & 2 ∇ $ 𝜌(𝑟) (3)
From the knowledge of G(r), V(r) can be determined from the virial equation (1).

According to equations ( 1), ( 2) and (3) it is possible to exploit the relation:

3 $ 0 $ = 4 $ 0 $ (1 - |6 $ | 4 $ ) (4)
where 𝐻(𝑟)/𝜌(𝑟) is the bond degree (BD) defined by Espinosa 15 that stands for the total energy per electron. The bond degree corresponds to the total electronic pressure and denotes the equilibrium between the pressure exerted by an electron on the others (𝐺(𝑟)/𝜌(𝑟)) and the pressure exerted by the other electrons on this one (𝑉(𝑟)/𝜌(𝑟)). |𝑉 7 |/𝐺 7 expresses the balance between the potential energy and kinetic energy densities for the formation of bonds.

On the basis of the electron density Laplacian and total energy density distributions, the local energy properties at the BCP allows one to analyse the interatomic interactions in crystals encompassing all types of bonding from the pure closed-shell interactions to the shared-shell ones. 1,16 𝐺 7 /𝜌 7 , which is an important parameter for the bonding classification, 11,15,17 is related to the evolution of 𝐻 7 /𝜌 7 vs. |𝑉 7 |/𝐺 7 (eq. 4). In the present work we use the latter two indicators to investigate the bonding features in Cu-Sb system intermediate compounds.

The Cu-Sb alloys are important for various applications, especially as lead-free solder in electronic products. The Cu-Sb system contains six intermetallic phases including Cu4Sb δ and Cu2Sb (η ones for which the compositional homogeneity ranges are significant in some of the Cu-Sb-based higher order systems. 18,19 Since the intermediate compounds in the Cu-Sb system are metallic conductors, their thermoelectric properties have never been investigated. The Seebeck coefficient has only been measured for Cu-rich two-phase Cu2Sb samples. 20 However the investigation of the thermoelectric properties in this system could be valuable to determine the specific influences of electronic transport on thermoelectric effect and how electronic motion generates Seebeck effect. In this work, three among the six aforementioned intermetallic phases, namely Cu2Sb (P4/nmm), 21 Cu3Sb (Fm3 2 m) 22 and Cu4Sb (P63/mmc), 23 have been computationally investigated. The choice for these three compounds has been dictated by the fact that they all present a large homogeneity range of compositions and temperatures, thus allowing for the formation of large ternary solid solutions. Moreover Cu4Sb δ and Cu2Sb (η are the only stable intermediate phases at room temperature. 18,19 The Cu2Sb, Cu3Sb and Cu4Sb compounds have been studied by analysing their electron density properties using the QTAIM approach, calculating their band structures and density of states by density functional theory (DFT) and predicting their thermoelectric properties using the Boltzmann transport theory.

COMPUTATIONAL DETAILS

The calculations have been performed at the generalized gradient approximation of the density functional theory using the Perdew-Burke-Ernzerhof (PBE) 24 exchange-correlation functionals. The Quantum Espresso (QE) package 25,26 has been used. Scalar-relativistic projector augmented plane wave pseudopotential have been employed and the cutoffs for the plane waves and electron density have been set to 48 Ry and 480 Ry, respectively. The Monkhorst-Pack procedure has been used to generate the k-points mesh for the Brillouin-zone sampling. The grid adopted in our calculations corresponds to a 4x4x4 k-points. The quality of the parameterization of the calculations can be assessed from the comparison between the optimized structure parameters and the experimental ones (Table 1). For Cu3Sb and Cu4Sb the calculated parameters agree nicely with the experimental ones with a relative deviation amounting to at most 2%. By contrast, for Cu2Sb a much larger relative deviation is observed with 6% for a and b parameters and 15% for c parameter. These discrepancies are hardly attributed to convergence insufficiency of the technical parameters (energy cutoff and k-point grid) since convergence tests have been conducted. Rather, since this type of discrepancies are often observed for tetragonal Bravais lattices with too large a and b parameters and too small c ones or the other way around, the discrepancies should be more likely due to DFT or exchange-correlation functionals weaknesses. The transport properties calculations have been performed with the BoltzTraP program. [START_REF] Madsen | [END_REF] For this purpose the number of k-points has been extremely increased to reach between six to ten thousand k-points in the Brillouin zone.

The thermoelectric properties have been determined as follows. The DFT-calculated bands are extrapolated using a smoothed Fourier expansion. Using the constant relaxation time approximation, the electrical conductivity and Seebeck coefficient can then be calculated according the following set of equations:

the Seebeck coefficient tensor reads:

𝑆 DE (𝑇, 𝜇) = F &' (=,G) H &( (=,G) (8)
The power factor has been calculated as the product PF=S 2 σ. The electron density has been analyzed using the program Critic2 29 based on QTAIM. The total, kinetic and potential energy densities have been determined as follows. The kinetic energy density has been calculated from the Kirzhnits expression eq. ( 3). From the knowledge of the electron density and its Laplacian, the potential energy density can be determined from the virial equation (1). Hence, the total energy density is obtained as the sum of the kinetic and potential ones eq. ( 2).

All the crystal structures presented in the following have been generated with VESTA. 30

RESULTS AND DISCUSSION

The Cu2Sb compound

In the tetragonal structure of Cu2Sb (P4/nmm), the copper atoms are located on the 2a (Cu2) and 2c (Cu1) sites, and the antimony atoms are located only on the 2c site as it can be seen in Figure 1. Along the z direction the crystal structure can be regarded as a "sandwich" layered structure. The central Cu2 plane is sandwiched between two atomic mono-layers, each being constituted of alternate Cu1 and Sb atoms, corrugated due to the strong bonding between the Cu1 atoms and the central Cu2 atoms as assessed from the high absolute value of the bond degree (pink point in Figure 2). The Cu1-Sb bonds in z direction connect these "sandwich" layers. This "sandwich" layered structure exhibits anisotropy with 2-D properties; as shown in Figure 3, besides a very high electrical conductivity, the electronic transport is more efficient in the x-y plane than along the z axis. Although electron transport displays some anisotropy ascribed to the distribution of Cu-Sb bonds, the energy density features exhibit only small differences (Figure 2), leading to a more or less equivalent electronic transport in the whole cell. This statement is in agreement with the Cu1 and Cu2 atoms having the same contributions around the Fermi level as shown in the atomic PDOS (Figure 4), which are dominated by similar 3d and 3p orbitals of both Cu1 and Cu2 atoms. This slightly anisotropic electronic transport could be not large enough to generate electronic entropy to drive the thermoelectric effect; indeed the Seebeck effect is negligible (around -2 μV/K at 300 K) near the Fermi level (Figure 5), which compares well with the experimental data (around 1μV/K at 300 K) reported in Cu-rich biphased samples. 20 The difference in sign of the Seebeck coefficient could be related to the presence of the secondary phase in the sample. The decent value of the Seebeck coefficient closest to the Fermi level is located at -0.18 Ry and originates from the large contribution of Cu1 and Cu2 atoms to the DOS. This negative value is hardly achievable by doping or alloying due to the impossibility to overcome such a large potential barrier. However, Cu2Sb can be regarded as an additional constituent of great importance for the improvement of thermoelectric materials. 31 Due to their layered structure, nanosized and thin films Cu2Sb are also promising as anode candidates for Na-ion batteries, showing high capacity and very high rate performance. [32][33][34] Moreover, near the absolute zero temperature, Cu2Sb can be a superconductor. 35

The Cu3Sb compound

The cubic Cu3Sb (𝐹𝑚3 2 𝑚) compound has also two positions for the Cu atoms, 8c for Cu1 and 4b for Cu2. The Cu1 atoms are bonded with four Sb atoms, whereas the Cu2 ones interact with six Sb atoms forming a regular octahedral configuration (Figure 6). The combined analysis of Figure 6 and Figure 7 shows that Cu and Sb atoms form two discriminative Cu-Sb bonds; the strong Cu1-Sb bonding along the < 110 > direction, so are the Cu1-Cu2 bonds, and the relatively weak Cu2-Sb bonding along the three axial directions. The same distribution of chemical bonding along the x-, y-, and z-direction implies the same electronic transport path and then the same electrical properties, in agreement with the results depicted in Figure 8. However, the different distributions of Cu1-Sb and Cu2-Sb bonds in the whole cell generate an electronic entropy which leads to a negative Seebeck coefficient slightly above the Fermi level and a corresponding peak of thermoelectric power factor (Figure 9). Indeed the Cu2 atoms have more prominent contribution to DOS than Cu1 atoms around the Fermi level (Figure 10a). This difference in DOS could be the intrinsic factor for driving thermoelectric effect. Indeed, the octahedral configuration of Cu2 with Sb atoms allows for a longer bond length, which is more ionic according to the smaller value of |𝑉 7 |/𝐺 7 , the charge density deformation (Figure 7), and the leading contribution of the Cu-2s orbital in the DOS (Figure 10b). By contrast, the Cu1 atoms have more obvious contribution of 3d orbitals in the DOS, and are more covalently bonded to Sb atoms with shorter bond length (Figure 10b and Figure 7). In fact, the Seebeck coefficient is related to the variation of the DOS near the Fermi level, corresponding to the small valley around +0.2 eV relative to the Fermi level in Figure 10. The negative Seebeck coefficient and the peak of the thermoelectric power factor in Figure 9 are a consequence of the DOS variation of the Cu2-s orbital from zero to maximum around the Fermi level (Figure 10b), having a close relationship with the weaker Cu2-Sb bonds.

In Figure 9, the peak of thermoelectric power factor almost locates around +0.015 Ry away from the Fermi level, corresponding to an approximate doping of one electron per unit cell (u.c.). This result indicates that the thermoelectric performance is more dramatically influenced by electron doping rather than hole doping. In order to further investigate the thermoelectric performance of Cu3Sb, the thermoelectric power factor is calculated at various electron doping levels (Figure 11).

When increasing the carrier concentration up to 1×10 21 e/cm 3 , the thermoelectric power factor is still very low and the thermoelectric effect is negligible. When the carrier concentration reaches 5×10 21 e/cm 3 , i.e. 0.6 e/u.c., the thermoelectric power factor has a tremendous increase in the whole temperature range. With a further increase in the electronic carrier concentration, the thermoelectric performance is severely degraded. By considering that the peak of thermoelectric power factor is maximal for one electron per unit cell, the optimum concentration should be around 8.4×10 21 e/cm 3 .

The Cu4Sb compound

The structurally optimized Cu4Sb has a large number of different interatomic interactions as shown in Figure 12. This crystal displays a four-layered structure along the z axis (Figure 13). In the first layer constituted by Sb1, Sb3 and Cu4 atoms, the strong Sb1-Sb3 bonding forms linear chains in the xy plane. The Sb2 atoms located in the third layer generate four weaker Cu-Sb bonds, among which two are forming the Cu-Sb linear chains in the xy plane involving Cu6 atoms, the remaining two are involving Cu3 atoms of adjacent copper layers. Obviously, the electronic transport differs depending on the direction (Figure 14), the conductivity around the Fermi level being much higher in the y direction than in the x and z ones, with the lowest value in the z direction. This electrical anisotropy originates from the distribution of the various interatomic interactions.

As mentioned above, the weaker Cu-Sb bondings solely involve the Sb2 atoms, being all located in the third layer. Therefore the Seebeck effect generated by these bonds cannot be extended to the whole cell and the Seebeck coefficient keeps rather small and fluctuant with the variation of chemical potential (Figure 15). In addition, the lowest electrical conductivity in the z-direction can be explained by the weak interlayer bonding; indeed only the Sb2 atom of the third layer is bonded with the Cu3 atom of the adjacent layer. By contrast, the electrical conductivity is larger in the xy plane, which can be explained by the presence of two layers containing only copper atoms, leading to a good conduction. The fact that the conduction in the y direction surpasses that in the x direction could be due to the presence of Cu-Sb-Cu chains oriented in the y-direction in the layers containing both Cu and Sb atoms.

CONCLUSIONS

The QTAIM approach coupled with the calculation of the electronic band structure, density of states and thermoelectric properties has allowed one to explain the origin and the anisotropy of the transport properties. It has been shown that the Seebeck effect in metallic systems, such as in Cu-Sb system compounds, is due to the presence of different distributions of chemical bonding in the whole cell. However a local difference of bonding distributions is not sufficient enough to generate effective thermoelectric effect.

  the unit cell volume, fμ(T,ϵ) the Fermi-Dirac distribution, μ the electron chemical potential, T the temperature, ϵ the electronic energy states, and σαβ(i,k) the electrical conductivity tensor for a given state i at kpoint k.𝜎 89 (𝑖, 𝑘) = 𝑒 $ 𝜏 A,B 𝑣 8 (𝑖, 𝑘)𝑣 9 (𝑖, 𝑘)(6)From the definition of the transport function: 𝜈 89 (𝑇, 𝜇) = & C=: ∫ 𝜎 89 (𝑖, 𝑘) :-;< % (=,?) ;?;

Fig. 1 .

 1 Fig. 1. Supercell of Cu2Sb (P4/𝑛𝑚𝑚). The a, b and c directions correspond to the x, y and z ones respectively. The small color balls are BCPs of different interactions.

Fig. 2 .Fig. 4 .

 24 Fig. 2. 𝐻 ! /𝜌 ! vs |𝑉 ! |/𝐺 ! for interactions in Cu2Sb. The dash line passes through the point (1,0) and that of the shortest Cu-Sb bond. The color words indicate the corresponding colors of the BCP in Fig. 1.

Fig. 5 .

 5 Fig. 5. Thermoelectric properties (Seebeck coefficient, electrical conductivity and power factor) as a function of the chemical potential at 300 K for Cu2Sb.

Fig. 6 .

 6 Fig. 6. Crystal structure of Cu3Sb (𝐹𝑚3 / 𝑚). The small pink balls are BCPs of weaker Cu-Sb bonds.

Fig. 7 .

 7 Fig. 7. Hb/ρb vs. |Vb|/Gb for interactions in Cu3Sb. The dash line passes through the point and that of the shortest Cu-Sb bond (black symbol). The insert shows the charge density deformation of Cu1-Sb and Cu2-Sb. Green and blue contour curves stand for positive and negative values, respectively. The red symbol corresponds to the weaker Cu2-Sb bond and the associated number corresponds to its interatomic distance in Å.

Fig. 8 .

 8 Fig. 8. Electrical conductivity as a function of chemical potential in different directions at 300 K for Cu3Sb.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Thermoelectric properties as a function of chemical potential at 300 K for Cu3Sb.

Fig. 11 .

 11 Fig. 11. Thermoelectric power factor as a function of temperature for various carrier concentrations (e/cm 3 ) in Cu3Sb.

Fig. 12 .

 12 Fig. 12. Hb/ρb vs. |Vb|/Gb for interactions in Cu4Sb. The dash line passes through the point (1,0) and that of the shortest Cu-Sb bond. The numbers correspond to interatomic distances in Å.

Fig. 13 .

 13 Fig. 13. Supercell of Cu4Sb. The small pink balls are BCPs of weaker Cu-Sb bonds and the red ones are BCPs of Sb-Sb bonds.

Fig. 14 .

 14 Fig. 14. Electrical conductivity as a function of chemical potential in different directions at 300 K for Cu4Sb.

Fig. 15 .

 15 Fig. 15. Thermoelectric properties as a function of chemical potential at 300 K for Cu4Sb.

Table 1 : Calculated structural parameters (Å) and volume (Å 3 ) of Cu-Sb system compounds. The corresponding experimental data taken from Ref. 27 are given in parenthesis.

 1

	Structure	a	b	c	Volume
	Cu2Sb	4.260 (4.002)	4.260 (4.002)	5.184 (6.104)	94.06 (97.8)
	Cu3Sb	5.934 (6.000)	5.934 (6.000)	5.934 (6.000)	209.0 (216)
	Cu4Sb	2.695 (2.688)	2.695 (2.688)	4.240 (4.328)	28.51 (27.08)

ACKNOWLEDGEMENTS

The PhD thesis of Mr. H. Yang is financially supported by the China Scholarship Council (CSC). This work was granted access to the HPC resources of the Centre Informatique National de l'Enseignement Supérieur (CINES), Montpellier, France under allocation A0070806881 made by the Grand Equipement National de Calcul Intensif (GENCI). It was also granted access to the HPC resources of Aix-Marseille University financed by the project Equip@Meso (ANR-10-EQPX-29-01) of the program "Investissements d'Avenir" supervised by the Agence Nationale de la Recherche.