First-principles study of electronic transport and structural properties of Cu$_{12}$Sb$_4$S$_{13}$ in its high-temperature phase

We present an ab initio study of the structural and electronic transport properties of tetrahedrite, Cu$_{12}$Sb$_4$S$_{13}$, in its high-temperature phase. We show how this complex compound can be seen as the outcome of an ordered arrangement of S-vacancies in a semiconducting fematinite-like structure (Cu3SbS4). Our calculations confirm that the S-vacancies are the natural doping mechanism in this thermoelectric compound and reveal a similar local chemical environment around crystallographically inequivalent Cu atoms, shedding light on the debate on XPS measurements in this compound. To access the electrical transport properties as a function of temperature we use the Kubo-Greenwood formula applied to snapshots of first-principles molecular dynamics simulations. This approach is essential to effectively account for the interaction between electrons and lattice vibrations in such a complex crystal structure where a strong anharmonicity plays a key role in stabilising the high-temperature phase. Our results show that the Seebeck coeffcient is in good agreement with experiments and the phonon-limited electrical resistivity displays a temperature trend that compares well with a wide range of experimental data. The predicted lower bound for the resistivity turns out to be remarkably low for a pristine mineral in the Cu-Sb-S system but not too far from the lowest experimental data reported in literature. The Lorenz number turns out to be substantially lower than what expected from the free-electron value in the Wiedemann-Franz law, thus providing an accurate way to estimate the electronic and lattice contributions to the thermal conductivity in experiments, of great significance in this very low thermal conductivity crystalline material.


I. INTRODUCTION
Within the family of copper-based semiconductors, the ternary Cu-Sb-S system has attracted great interest in recent years due to a variety of appealing structural, electronic and thermal transport properties.Indeed, Cu-Sb-S compounds display a rich structural variety, a large range of band gaps and are characterised by extremely low thermal conductivities.These features combined with the non-toxicity and abundance of the constituent elements make the Cu-Sb-S system an ideal playground to explore and optimise materials for sustainable thermoelectric devices.
Within this system, the tetrahedrite Cu 12 Sb 4 S 13 is a compound of particular interest as in its pristine form it has a remarkably high zT value, approximately 0.6 at 700 K, that is the result of a very low thermal conductivity (below 1 Wm −1 K −1 from 300 to 700 K) and a high power factor 1 .
In spite of a large effort in characterising the material (see, for instance, Refs.1-5), a detailed understanding of the origin of these transport characteristics is still missing-possibly preventing further optimisations of the compound as well as the development of design ideas for novel thermoelectric materials.For instance, the relative electronic and lattice contributions to such a low thermal conductivity are still under debate.In addition, it is not clear how the carrier scattering mechanisms determine the electronic transport and especially the temperature dependence of the electrical resistivity that appears to vary considerably in the available experimental data.In this, of particular interest is the effect on the electronic carriers of a very peculiar lattice dynamics, that shows soft phonons, strong anharmonicity and unusually large anisotropic atomic displacement parameters as shown in diffraction studies 6 .
In order to address these important issues we have initially carried out an extensive study of the structural stability of the high-temperature phase of tetrahedrite.In particular, we have done this by analysing the structural relationship between tetrahedrite and the simpler fematinite structure (Cu 3 SbS 4 ), which under optimal doping shows as well quite good thermoelectric properties 7,8 .As we will show, tetrahedrite can be seen as a fematinite-based structure that accommodates an ordered arrangement of S-vacancies.On one side this sheds new light on the caged nature of tetrahedrite and on its anharmonic structure 6 that are at the origin of the low lattice thermal conductivity.On the other side, the S-vacancies are the natural doping mechanism of the Cu-Sb-S compound; indeed, the starting Cu 3 SbS 4 structure is a semiconducting material with a gap of about 0.8 eV and it shows a zT of 0.6 at 600 K under Ge or Sn doping.
We have then investigated the temperature dependence of the electronic transport properties of the compound.For this we have used the Kubo-Greenwood (K-G) formula 9,10 applied to snapshots of first-principles molecular dynamics simulations at different temperatures.The use of this approach instead of a Boltzmann transport formalism is here motivated not only by the complex crystal structure of the compound but also by the fact that this method inherently incorporates the interaction between electrons and lattice vibrations as well as, very importantly, anharmonicity, that is supposed to stabilise the high-temperature phase of tetrahedrite for which zero-temperature Density Functional Theory (zero-T DFT) predicts unstable phonons 1 .As we will show, this study reveals the intrinsic transport properties of the compound and their temperature dependence: the Seebeck coefficient is in very good agreement with experimental data and the lattice-limited electrical resistivity shows values that are fairly temperature independent, remarkably low but not too far from the lowest experimental value reported in literature 3 .The Lorenz number turns out to be significantly lower than the free electron value, an important result to quantify lattice and electronic contributions to thermal conductivity.
The structure of the paper is as follows.In the Sec.II, we present the theoretical framework and computational details.The structural characterization of tetrahedrite is presented and discussed in Sec.III.In Sec.IV we analyse the electronic transport properties of the compound.

II. THEORETICAL APPROACH AND COMPUTATIONAL DETAILS
The theoretical calculations were performed within DFT in the Kohn-Sham scheme 11,12 (DFT-KS) and using the projector augmented wave method (PAW) 13 as implemented in the Vienna ab-initio simulation package (VASP) 14,15 .We used the Perdew, Burke and Ernzerhof (PBE) exchange-correlation functional 16 for all the calculations and a plane wave energy cutoff of 550 eV; we used a 3 × 3 × 3 grid of k-points to sample the Brillouin zone of the conventional unit cell.Lattice parameters and internal positions were fully relaxed (but constraining the cell to be tetragonal).
Formation energies are calculated via the formula where E 0 is the ground state total energy of the system under consideration that has N a atoms in the unit cell, N t is the number of different types of atoms in the unit cell, and n i and E i are the number of atoms of type i and their ground state total energy respectively.
Core binding energies were computed within the ∆ Self-consistent Field (∆SCF) approximation 17,18 at the DFT-KS level of theory.The ∆SCF approximation includes effects from the relaxation of the valence electrons after the core ionization has taken place, reflecting the screening of the core hole created in the photoemission process.Although a core level binding energy is evaluated as the total energy difference between a first core ionized state and the ground state in a spin-polarized calculation, here we follow a spin-paired atomic approximation.To make the unit cell neutral, a compensating electron charge is introduced into the first conduction band.This seems like a good approximation for metals where core-hole screening is efficient.Furthermore, a relativistic Full Potential Linear Muffin-Tin Orbital (FP-LMTO) treatment of the magnetic quantum number 19 has been implemented in the QUESTAAL package 20,21 and used to evaluate the core binding energies of the Cu 2p 1/2 and 2p 3/2 electronic states.The smoothed LMTO basis used in this work includes atomic orbitals with l ≤ l max = 4 with orbitals 4d of Sb and Cu added in the form of local orbitals.
We compute the dynamic Onsager coefficients L ij (ω) in the static limit for each ionic configuration in the supercell according to the Kubo-Greenwood (K-G) formula [22][23][24] : where f ( ) is the Fermi-Dirac distribution function, µ is the chemical potential, V is the volume of the simulation cell, α = x, y, z and ψ l,k are the KS orbitals (at band l and wave vector k) with corresponding energies l,k .When n,k − m,k goes to zero (the case of intra-band transitions and degeneracies), The electrical conductivity is given by σ = L 11 while the Seebeck coefficient is and the electronic contribution to the thermal conductivity is The temperature dependence of the transport coefficients is obtained by averaging over snapshots computed via molecular dynamics (MD) simulations 25 .We used simulation cells of 232 atoms (2 × 2 × 1 supercells).Initailly we performed NPT runs using the Parrinello-Rahman scheme with implementation of the Nosé-Poincaré approach for isothermal sampling 26 ; we used a timestep of 1 fs and a period for the thermostat of 1.11 ps; the barostat is used with fictitious masses of 10 −3 amu.These MD simulations were run for 5 ps and this time was enough to calculate the averaged cell parameters at different temperatures.As a second step, we performed NVT runs using a Langevin thermostat.We used a drift parameter of 1.0 ps −1 and performed simulations of about 8.5 ps with a timestep of 1 fs; the first 4.5 ps have been used to safely reach thermal equilibration.After equilibration we extracted snapshots of nuclear positions every 500 MD steps.For every snapshot we have performed static electronic DFT calculations and used the output to evaluate Eq. 1.While in the ab inito MD simulations we used a k-points grid of 1 × 1 × 3, for every snapshot we used a denser grid of 3 × 3 × 6; for the smearing we used the Fermi-Dirac function with the electronic temperature corresponding to the temperature of the system.The chemical potential obtained from the ab initio calculation was then consequently used in Eq. 1.The Dirac delta functions in Eq. 1 have been approximated with gaussians of spread 40 meV.The Gaussian broadening is chosen as small as possible in order to remove the small oscillations in the optical conductivity that are due to the discretization of band structure. 27For instance, at 300K the resistivity changes by less than 5% when we change the smearing from 20 meV to 70 meV.It is also worth to stress that we used a cell of 232 atoms that is tetragonal; the conductivity tensor is however diagonal with the proper symmetry for a cubic system, a quite clear indication of a good convergence of the transport coefficients with respect to the cell size.
The electronic transport parameters were also computed using the Boltzmann transport equation within the constant relaxation time approximation (BTE-CRT) using the Wannier interpolation scheme implemented in the BoltzWann code 28 .For these calculations we used the DFT band structure of the symmetric high-temperature phase of tetrahedrite.
atoms are four-fold coordinated as they are at the centre of tetrahedrons of S(1) atoms, whereas the Cu(2) atoms are three-fold coordinated and are located at the vertices of octahedrons at the centre of which there are S(2) atoms.The Sb atoms are on [001] planes at the apex of trigonal pyramids with S(1) atoms.
In order to gain a better understanding of the structure of this complex compound, it is useful to establish a connection between the crystal structure of tetrahedrite and the simpler one of fematinite, a very stable mineral in the Cu-Sb-S system.To do this we consider a 2 × 2 × 1 supercell of fematinite structure, containing 64 atoms (Fig. 2a).This supercell is almost cubic as the primitive cell of fematinite is tetragonal (I 42m) with experimental cell parameters a = 5.39 Å and c = 10.75 Å8 ; in addition, the average lattice parameter of the supercell is quite close to the one of tetrahedrite.At this point, in order to establish a structural link between the tetrahedrite and the fematinite structures, one can i) swap Cu and Sb atoms in the highlighted (001) plane as shown in Fig. 2 (b) (or, equivalently, slip the highlighted plane in the < 110 > direction), and then ii) remove 6 S atoms from the supercell [see Fig. 2 (c)], thus recovering the stechiometry of tetrahedrite.When the system created in this way is relaxed within zero-T DFT, the resulting structure (Fig. 2 (d)) is cubic (with a lattice parameter of 10.39 Å) and quite similar to the one of tetrahedrite.The only difference is a small distortion of the octahedral structures that is not observed in the high-T phase of tetrahedrite (see Fig. 1).This result is not surprising since, in tetrahedrite at low temperatures, there are hints of a structural transition happening around 70 K that modifies the size of the cell and introduces distortions in the local symmetry for T < 70 K; this transition is object of an active experimental investigation 2,31-34 .In addition, the instability at low temperatures emerges also from DFT calculations of the phonon spectrum of the high-T symmetric phase 1 .
To conclude the analysis of the link between tetrahedrite and fematinite, it is interesting to compare the stability of the structures suggested in Fig. 2. Our calculations show that the formation energy of the initial and final structures in Fig. 2 are the same, -0.241 eV, while the formation energy of the high-T symmetric phase (Fig. 1) is 2 meV higher.It is also interesting to observe that the formation energy of the intermediate step where Cu and Sb are exchanged [the structure of Fig. 2 (b)] is only 6 meV higher than the one of fematinite, and that in this structure the clustering of sulphur vacancies in neighbouring positions turns out to be very convenient.These results not only strengthen the connection between tetrahedrite and fematinite but also directly show that the Cu-Sb-S crystalline framework can be quite flexible and open to a variety of energetically competitive structures.

IV. ELECTRONIC STRUCTURE AND XPS
Our calculation for the electronic structure of the symmetric high-temperature phase of tetrahedrite is in agreement with the DFT results reported in literature [see for instance Ref. 1].The density of states (DOS) in Fig. 3 clearly shows that the compound is metallic but it can be seen as an heavily doped p-type semiconductor.Indeed, the Fermi level lies almost at the top of the highest occupied band-the valence band that is mostly formed by sulfur 3p and copper 3d hybridised states; the bottom of the next available empty band-the conduction band-lies at about 1.25 eV from the Fermi level.This band gap is in agreement with other published theoretical results 1 , but it is lower than the experimental gaps reported in literature that vary between 1.7 eV and 1.9 eV. 35,36he electronic structure of tetrahedrite turns out to be very different from the one of fematinite.Indeed, as discussed in Refs.7 and 8, fematinite is a semiconductor with a band gap of around 0.6 eV and the DOS at the top of the valence band is rather smooth and it can be easily fitted with a proper density of states effective mass that can be used in a simple parabolic band model to predict the thermoelectric properties of the compound when p-doped.The symmetric phase of tetrahedrite displays instead a rather complex and spiky DOS around the Fermi level, making very difficult the use of simple models to predict the transport properties.
In Fig. 3 we also show the effect of thermal motion of atoms on the electronic structure of the compound.For this, we present the DOS of tetrahedrite at 300 K and 600 K, obtained by averaging the DOS of snapshots of the MD simulations.Our results clearly show that thermal motion leads to a strong flattening of the sharp peaks of the DOS of the zero-T symmetric structure as well as to a sensible decrease of the energy gap between the Fermi level and the bottom of the conduction band.
For this compound it is also important to point out that the chemical composition is often described by the formula Cu + 10 Cu 2+ 2 Sb 3− 4 S 2− 13 37 , however, the existence of two so different copper atoms (and in different proportion) is still an open question.Indeed, studies of the chemical nature of copper have been carried out by many groups using photoemission spectroscopy (XPS) and X-ray absorption spectroscopy measurements, and some of them report the presence of a Cu 2+ ions 5,38 , while others show that both Cu(1) and Cu(2) atoms behave like monovalent ions 31 .
Our DFT-KS calculations within the ∆SCF scheme show that the two copper atoms in the two inequivalent sites have binding energies of the 2p core states that are very similar and are almost unaffected by the structural distortion discussed in the previous section.In particular, the DFT binding energies of the Cu 2p 1/2 and 2p 3/2 states are 951.7 eV and 930.5 eV for the Cu(1) atom and 953.5 eV and 932.3 eV for the Cu(2) atom.These values are in very good agreement with the two main peaks at 952 eV and 932 eV observed in XPS experiments.These energies are considered the fingerprint of monovalent Cu ions, as they are very similar to what is observed in Cu 2 O, but they are very different from the spectral peaks in the ranges 940-945 eV and 960-965 eV due to divalent Cu ions observed in CuO.Our results for pristine tetrahedrite in its room temperature phase is therefore compatible with the experimental findings of Tanaka et al. 31 and this suggests that the existence of divalent copper atoms reported in literature might be of extrinsic nature.

V. THERMOELECTRIC TRANSPORT COEFFICIENTS
In this section we discuss the transport properties of the compound and their temperature dependence.In Fig. 4 we compare the intrinsic resistivity ρ = 1/σ of tetrahedrite computed with the K-G approach with the experimental results.Most of the experimental data show a slow increase of the resistivity above room temperature and this trend is in good agreement with our theoretical results.The experimental data are distributed on a fairly broad range of values as a result of extrinsic effects (e.g.defects, impurities and polycrystallinity) that in different samples can be of different nature and concentrations.However, our values for the phonon-limited resistivity (about 0.5 × 10 −5 Ω m at 600 K) are not too far from the lowest experimental values reported in literature 3,5 , about 1 × 10 −5 Ω m.These values are quite low for an undoped mineral in the Cu-Sb-S family.Indeed, for instance, we recently showed that nano structured Sn doped fematinite at optimal doping for thermoelectric efficiency (a carrier concentration of about 5 × 10 20 cm −3 ) displays a resistivity between 1 × 10 −5 Ω m and 1.8 × 10 −5 Ω m between 300K and 600K 8 .
It is also interesting to observe that the temperature dependence of our theoretical results and of the experimental results with the lowest resistivity 3,5 can also be reproduced fairly well by tuning the carrier lifetime τ in the BTE-CRT approach that uses the DFT band structure of the symmetric high-temperature phase of tetrahedrite.For instance, as shown in Fig. 4, the experimental data in Ref. 3 can be fitted in a broad temperature range using a carrier lifetime due to both impurities and phonons of 5.6 fs.The same fitting can be done on our K-G theoretical results and in this case the carrier lifetime is 16 fs; this is due only to interactions between electrons and lattice vibrations and it quantifies the upper bound of the carrier lifetime in pristine tetrahedrite.
Quite interestingly, the fact that the temperature dependence of the resistivity is captured by the BTE-CRT approach suggests that the key ingredient to explain this trend is the electronic structure of the compound, while the carrier scattering time due to the electron-phonon coupling seems to depend very weakly on temperature in this material.
In Fig. 5 we show our results for the Seebeck coefficient as a function of temperature.The theoretical data are in very good agreement with experiments, and the K-G formalism gives values that are very similar to the ones obtained with the BTE-CRT approach.The agreement between the exact and the approximate theoretical approaches is not novel (for instance, see the analysis done for simple semiconductors 39,40 ) and too surprising because, as we discussed above, a transport coefficient as the resistivity can be reproduced in a wide range of temperatures using the equilibrium band structure of the symmetric phase and a carrier relaxation time independent of temperature; the possibility to use a constant relaxation time τ justifies to a good extent the BTE-CRT approximation in which the Seebeck coefficient (as well as the Lorenz number discussed below) turns out to be independent of τ .It is important to notice that the values of the Seebeck coefficient are appreciably high in this pristine system.For instance, the predicted value at 600 K, about 100 µV/K, is almost as high as the values found in fematinite at optimal Sn doping, about 130 µV/K 8 .These values of the Seebeck coefficient and a low resistivity lead to an undoped compound with a quite high power factor, S 2 /ρ.The highest reported experimental values of the power factor at 600 K are between 1.2 and 1.4 mW/(K 2 m) 3,5,41 , but our calculations suggest the possibility of further improvement.Indeed, the intrinsic upper bound for the power factor predicted here at 600 K is about 2 mW/(K 2 m).
In Fig. 6 we show the results for the Lorenz number, L = κ e /(σT ).This is a key quantity to estimate the electronic contribution κ e to the thermal conductivity from the experimental values of the electrical conductivity.Once κ e is known, it is possible to extract the lattice contribution to the thermal conductivity of the material, κ L , from the experimental values (κ exp ), i.e. κ L = κ exp − κ e .
Our calculations show that even though, as discussed before, tetrahedrite displays a metallic character above room temperature, the values of L obtained with the K-G formalism are almost temperature-independent but substantially lower than what expected from the free-electron value in the Wiedemann-Franz law (L 0 = 2.44 × 10 −8 WΩ K −2 ).In addition, contrary to what observed for the Seebeck coefficient, the BTE-CRT approach gives values for L that are slightly higher than the K-G ones and with a slightly different temperature trend.Here, the difference between the results from the K-G and BTE-CRT approaches is not surprising because, as shown in other model 42 and ab initio 39 calculations, the Lorenz number usually tends to be quite sensitive to the details of the carrier scattering mechanisms.
These theoretical predictions for L are interesting and very valuable as they show that the electronic contribution to the thermal conductivity approximated using the free-electron value for L (as done, for instance, in Refs. 1 and 3) can be significantly overestimated, up to about 50%; this thus can lead to a marked underestimation of the intrinsic lattice thermal conductivity of tetrahedrite.Indeed, for instance, in Ref. 1 the reported value of thermal conductivity is of about 1.45 W/(m K) at 600 K and the lattice contribution estimated using L 0 is of about 0.4 W/(m K); if instead the estimation is done using L from the K-G approach, we get a value for κ L of about 0.8 W/(m K) (κ L ≈ 0.6 W/(m K) using L from the BTE-CRT approach), suggesting that in experiments the lattice thermal conductivity of the compound is certainly low, but it can be very similar to the electronic contribution.

VI. CONCLUSIONS
We have presented a detailed ab initio study of the structural and electronic transport properties of the hightemperature phase of tetrahedrite.We have shown that the structural variety of the Cu-Sb-S network allows a description of tetrahedrite in terms of a fematinite-like crystal modified by an ordered arrangement of S-vacancies.While this structure presents two crystallographically inequivalent copper atoms, our calculations have shown that these atoms experience similar local chemical environments.Indeed, they have very similar core level binding energies and this helps the analysis of XPS measurements in this compound.The Kubo-Greenwood approach has allowed us to predict two important thermoelectric quantities, the phonon-limited electrical resistivity and the Lorenz number.The predicted resistivity turns out to be quite low for an undoped compound in the Cu-Sb-S system.The lowest experimental data reported in literature are not too far from the predicted intrinsic values: this clearly shows that the quality of the samples studied in experiments is high, but it also suggests the possibility of further improvements of the electronic transport properties.The Lorenz number has turned out to be substantially lower than what expected for the free-electron value, often used to estimate the electronic and lattice contributions to the thermal conductivity in experiments.Thus our result provides a more accurate reference to analyse thermal transport in this compound.Finally, it is important to stress that the K-G approach has been key to predict transport properties in this system.Indeed, the zero-temperature DFT phonon calculations for the high temperature phase result in imaginary frequencies and lattice instabilities, which make hard to apply the exact BTE approach that includes electron-phonon coupling.Our analysis has also allowed us to show that the use of a less computationally demanding BTE-CRT approach is quite effective to reproduce the temperature trends of transport quantities such as the resistivity and the Seebeck coefficient in this complex compound.

FIG. 1 .
FIG. 1. Tetrahedrite cubic conventional cell.S atoms are yellow, Cu atoms are blue and Sb atoms are orange.

FIG. 3 .
FIG. 3. Comparison getween the DOS of the symmetric structure in Fig. 1 (solid blue line) and the snapshot averaged DOS at 300 K (orange dashed line) and at 600 K (red dotted line).EF is the Fermi level.

FIG. 4 .
FIG. 4. Resistivity as a function of temperature.The blue squares (green circles) are the K-G (BTE-CRT) results.The black symbols are the experimental data: the down triangles are from Ref. 2, the up triangles from Ref. 1, the diamonds from Ref. 3, the stars from Ref. 4, and the plus symbols from Ref. 5. The error bars in the K-G results are smaller than the symbols used.

FIG. 5 .
FIG. 5. Seebeck coefficient as a function of temperature.The symbols are as in Fig. 4.