Hailong Yang

Pascal Boulet
email: pascal.boulet@univ-amu.fr

Marie-Christine Record

New Insight Into the Structure-Property Relationships from Chemical Bonding Analysis: Application to Thermoelectric Materials

Keywords: Bond theory, Density functional Theory, Electronic structure, Thermoelectrics

Structure-property relationships are indisputably important to deeply understand the nature of materials.

The main direction pursued to investigate these relations is based on the chemical bonding concept. The integral and local properties of chemical bonding determined from the topological analysis of electron density by combining the quantum theory of atoms in molecules (AIM) and DFT calculations can aid to build a bridge between the atomic structure and the intrinsic properties of materials. In this paper, the electronic and thermal transport properties of thermoelectric materials in the Cu-Sb-Se system have been analysed by coupling energy densities relationships and corresponding electron density Laplacian distribution for specific interactions. The asset of this method is to directly highlight the influence of lone pair electrons and atomic fluctuations on interatomic interactions as well as the coupling effect of weak bonds, which offers a new approach to design and search for materials with low thermal conductivity at low calculation cost.

Introduction

Thermoelectricity has proved to be a promising route for harvesting wasted energy by converting heat into power.

The performance of thermoelectrics is characterized by the dimensionless figure of merit, ZT=(S 2 σ/к)T, where S, σ, к, and T are the Seebeck coefficient, electrical conductivity, thermal conductivity and absolute temperature, respectively. However, these parameters are severely mutually coupled, in such a way that it is hardly possible to improve one of the parameters without degrading others. In order to optimize ZT, many strategies have been developed since the 1960's, such as band-, phonon-and nano-engineering [START_REF] Pei | Convergence of electronic bands for high performance bulk thermoelectrics[END_REF][START_REF] Kim | Strategies for engineering phonon transport in thermoelectrics[END_REF][START_REF] Chen | Manipulation of Phonon Transport in Thermoelectrics[END_REF][START_REF] Boona | Nanomagnets boost thermoelectric output[END_REF], but the understanding of the structure-properties relationships is yet insufficient to rationally design outstanding thermoelectric materials. In order to perceive the distinctions in properties between different crystal structures, various models have been developed on the basis of covalent bonding [START_REF] Phillips | A Posteriori Theory of Covalent Bonding[END_REF], dielectric constants [START_REF] Phillips | Dielectric Definition of Electronegativity[END_REF] partially ionic binding [START_REF] Phillips | Covalent Bond in Crystals. II. Partially Ionic Binding[END_REF][START_REF] Phillips | Dielectric Classification of Crystal Structures, Ionization Potentials, and Band Structures[END_REF] and ionicity [START_REF] Stiles | Trends in the ionicity in the average valence V materials[END_REF]. However these models are insufficient to correctly describe the structure-property relationships. For instance, the lattice thermal conductivity (кL) is not only related to the ionicity of bonds [START_REF] Ioffe | Thermal conductivity of semiconductor solid solutions[END_REF], but also correlated to the ion coordination and vacancies [START_REF] Spitzer | Lattice thermal conductivity of semiconductors: A chemical bond approach[END_REF]. Moreover, in non-normal valence compounds, both coordination and vacancies are difficult to define.

Following the discovery of abnormal interatomic distances, bond angles, large polarizability and unexpected electronic and vibrational properties of IV-VI compounds and group-V elements attributed to the resonance of bonding [START_REF] Lucovsky | Effects of Resonance Bonding on the Properties of Crystalline and Amorphous Semiconductors[END_REF], the latter one was confirmed to have an important role in the low thermal conductivity, which results from phonon-phonon scattering [START_REF] Lee | Resonant bonding leads to low lattice thermal conductivity[END_REF]. The giant phonon anharmonicity in phosphorene resulting from long-range interactions also originates from the resonant bonding [START_REF] Qin | Resonant bonding driven giant phonon anharmonicity and low thermal conductivity of phosphorene[END_REF]. Recently, the analysis of low thermal conductivity of thermoelectrics has been performed by using the approach of chemical-bond hierarchy and it was shown that the phonon-phonon anharmonic scattering originates from weak chemical bonding, while the extend of weak chemical bonds determines the lattice thermal conductivity [START_REF] Qiu | Part-crystalline part-liquid state and rattling-like thermal damping in materials with chemical-bond hierarchy[END_REF][START_REF] Qiu | Diverse lattice dynamics in ternary Cu-Sb-Se compounds[END_REF]. This chemical approach has even been proposed for designing thermoelectric materials [START_REF] Zeier | Thinking Like a Chemist: Intuition in Thermoelectric Materials[END_REF].

Trends in electronic properties and lattice thermal conductivity can be explained from electronegativity, bandwidth, orbital overlap, bond energy, bond length, and bond strength. As a consequence, the chemical approach is undoubtedly a route for developing thermoelectric materials.

Even though the required tools for analysing structure-property relationships of solids are still limited, in a previous work we analysed the chemical bonding in various systems by combining the integral with local properties of bonding [START_REF] Yang | A rapid method for analyzing the chemical bond from energy densities calculations at the bond critical point[END_REF]. This method is a simple one to characterize specific bondings in real space through density functional theory (DFT) calculations and topological analysis of the electron density on the basis of the quantum theory of atoms in molecules (QTAIM).

Since the chemical bonding is responsible for the physical and chemical properties of materials and that the bonding information can be directly obtained from the charge density, which is an observable in physics, the relationships between bonding and thermoelectric properties could be clarified by analysing the charge density of the compound of interest.

The Cu-Sb-Se ternary system is a suitable platform for analysing bonding and properties relationships since it contains three stable thermoelectric materials, namely CuSbSe2, Cu3SbSe3 and Cu3SbSe4, having the same constitutive elements but discriminative thermoelectric, both electronic and thermal, properties. Cu3SbSe4 has a high electrical conductivity in the whole range of temperatures; Cu3SbSe3 has an anomalously low lattice thermal conductivity, which is almost independent of temperature; CuSbSe2 has a high Seebeck coefficient in the whole range of temperatures [START_REF] Skoug | High thermoelectric figure of merit in the Cu3SbSe4-Cu3SbS4 solid solution[END_REF][START_REF] Skoug | Lattice thermal conductivity of the Cu3SbSe4-Cu3SbS4 solid solution[END_REF][START_REF] Skoug | Role of lone-pair electrons in producing minimum thermal conductivity in nitrogen-group chalcogenide compounds[END_REF]. In addition, Cu4SbSe5 and Cu12Sb4Se13, which were predicted by theoretical studies [START_REF] Zhang | Prediction of new stable compounds and promising thermoelectrics in the Cu-Sb-Se system[END_REF], have been found as promising thermoelectric materials, in particular for Cu12Sb4Se13, which presents an extremely low thermal conductivity. It should however be mentioned that, these compounds have never been experimentally evidenced. The Cu12Sb4S13-xSex solid solution only exists for 0<x<3 [START_REF] Lu | Phase Stability, Crystal Structure, and Thermoelectric Properties of Cu12Sb4S13-xSex Solid Solutions[END_REF].

In this paper, the investigation of structure-property relationships focusing on electronic and thermal transport, for the three existing stable compounds CuSbSe2, Cu3SbSe3, Cu3SbSe4 and for the two predicted ones Cu4SbSe5 and Cu12Sb4Se13 [START_REF] Zhang | Prediction of new stable compounds and promising thermoelectrics in the Cu-Sb-Se system[END_REF], hereafter named Se2, Se3, Se4, Se5 and Se13, respectively, has been performed by analysing the chemical bonding derived from the charge density.

Computational Methods

Density functional theory (DFT) calculations were performed using the QUANTUM ESPRESSO (QE) simulation package [START_REF] Giannozzi | QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials[END_REF][START_REF] Giannozzi | Advanced capabilities for materials modelling with Quantum ESPRESSO[END_REF]. The electronic exchange-correlation energy was treated by the generalized gradient approximation (GGA) using the Perdew-Burke-Ernzerhof (PBE) [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] functionals and the interaction between ions and electrons was described by the projector augmented wave (PAW) method [START_REF] Blöchl | Projector augmented-wave method[END_REF]. All the electron density was calculated with a kinetic energy cutoff for wavefunctions of 48 Ry and the Monkhorst-Pack procedure was used to generate a 4x4x4 k-points mesh for the Brillouin-zone sampling. The charge density was subsequently analyzed by the program Critic2 [START_REF] Otero-De-La-Roza | Critic2: A program for real-space analysis of quantum chemical interactions in solids[END_REF] that implements Bader's Quantum Theory of Atoms in Molecules (QTAIM) [START_REF] Bader | Atoms in Molecules: A Quantum Theory[END_REF]. In this work, the model used in Critic2 for the kinetic energy density calculation is based on the Thomas-Fermi equation with the semiclassical gradient correction proposed by Kirzhnits [START_REF] Kirzhnits | Quantum corrections to the Thomas-Fermi equation[END_REF][START_REF] Kirzhnits | Field Theoretical Methods in Many-Body Systems[END_REF]. The local energy densities of chemical bonding are interrelated by the local Virial theorem through the electron density Laplacian Ñ 2 r(r) [START_REF] Bader | Atoms in Molecules: A Quantum Theory[END_REF],

! " ∇ # 𝜌(𝑟) = 2𝐺(𝑟) + 𝑉(𝑟) (1)
𝐺(𝑟) + 𝑉(𝑟) = 𝐻(𝑟) (2)
where G, V and H are the kinetic, potential and total energy densities, respectively. At the bond critical point (BCP) eq.( 2) can be recast as:

(1 -

|( ! | & ! ) (3)
Eq.( 3) displays the variation of Hb/rb along with |Vb|/Gb: the former indicator stands for the total energy per electron related to the total electronic pressure and the latter one exhibits the competition between potential energy and kinetic energy for bonding formation. These two indicators corresponding to specific characteristics of bonding can be used for the bonding classification. At the same time, the kinetic energy per electron Gb/rb is related to the polarization and the rigidity of the interatomic interactions, which can reveal the bonding environment in the crystal structure [START_REF] Yang | A rapid method for analyzing the chemical bond from energy densities calculations at the bond critical point[END_REF].

In order to further investigate the long-range or weak interactions, the source function (SF) [START_REF] Bader | A Green's function for the density[END_REF][START_REF] Gatti | Chemical information from the source function[END_REF][START_REF] Gatti | Source function description of metal-metal bonding in d-block organometallic compounds[END_REF] has been used

to estimate the relative contribution of an atom or group of atoms to the electron density at any point in a system, which can disclose the local and nonlocal character of the electron density distributions.

The SF, which is implemented in a modified version of Critic2 [START_REF] Tantardini | Source function and plane waves: Toward complete bader analysis[END_REF], is an interesting complementary tool to provide additional chemical information from the whole crystal.

In the following, the crystal structures are visualized with VESTA [START_REF] Momma | VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data[END_REF] and the figures depicting the electronic shells come directly from Critic2.

Results and Discussion

Electronic transport

The excellent electrical properties of Se4 are attributed to its crystal structure, consisting of the three-dimensional Cu/Se framework [Cu3Se4] acting as the hole conduction pathway [START_REF] Tantardini | Source function and plane waves: Toward complete bader analysis[END_REF]. The valence band maximum (VBM) of Se4 contains mainly a mixing of Cu-3d and Se-4p and the conduction band minimum (CBM) is dominated by the mixing of Sb-5s and Se-4p [START_REF] Yang | New stannite-like p-type thermoelectric material Cu3SbSe4[END_REF][START_REF] Do | Physics of bandgap formation in Cu-Sb-Se based novel thermoelectrics: the role of Sb valency and Cu d levels[END_REF]. This Cu-3d with Se-4p mixing stabilizes holes and provides high mobility Cu/Se framework paths [START_REF] Yang | New stannite-like p-type thermoelectric material Cu3SbSe4[END_REF]. As the crystal structure of Se4 is tetragonal with high symmetry (space group I-42m [START_REF] Pfitzner | Crystal structure of tricopper tetraselenoantimonate (V), Cu3SbSe4[END_REF], According to this analysis, the location of the Cu-Se bonds for Se5 and Se13 in Fig. 3 suggests a good electrical conductivity for these compounds, which agrees well with their semi-metallic character as reported in Ref [START_REF] Zhang | Prediction of new stable compounds and promising thermoelectrics in the Cu-Sb-Se system[END_REF], that of Se13 being less pronounced than that of Se5.

According to Yang et al. [START_REF] Yang | New stannite-like p-type thermoelectric material Cu3SbSe4[END_REF], the anisotropy, contributing to appropriate transport properties, in electronic and phononic transport of Se4 originates from different Cu-Se bond lengths. But the inconsistent relationship between interatomic distances and local properties of Cu-Se and Sb-Se interactions indicates that bond length only is not sufficient to evaluate the properties of bonds, especially for the bonds in semiconductors with charge transition. As a consequence the topological analysis of VSCC and local, energy density properties of bonding is the most appropriate method to explore the bonding behaviour in complex systems. Both electronic and thermal transports are severely reduced in dispersive and non-local ionic bonds, however in Se4 there is a high similarity in bonding properties for all Cu-Se interactions as well as for all Sb-Se ones, which is the key to both electrical and thermal transport efficiency. The VSCCs analysis can explain the bonding properties but in order to compare the transport properties for different compounds we combine this analysis to the viewpoint of electronic conduction units (ECU) introduced by Bu. et al. [START_REF] Bu | Observation of High Seebeck Coefficient and Low Thermal Conductivity in [SrO]-Intercalated CuSbSe2[END_REF].

For this, two kinds of building units are proposed in this work to describe the electronic transport in the semiconducting compounds of interest, namely one corresponding to the ionic bonds with overlapped VSCCs and the other corresponding to the ionic ones with separated VSCCs. All information on atomic interactions in Se4, Se3

and Se2 has been listed in Tab.1, Tab. In Tab.1, Tab.2 and Tab.3 the values of the kinetic energy per electron Gb/𝜌b are also reported. One can see that, the largest values of Gb/𝜌b correspond to overlapped VSCC, which can be correlated with the bond order as defined by Cremer and Kraka [START_REF] Cremer | Chemical Bonds without Bonding Electron Density? Does the Difference Electron-Density Analysis Suffice for a Description of the Chemical Bond?[END_REF][START_REF] Cremer | A Description of the Chemical Bond in Terms of Local Properties of Electron Density and Energy[END_REF]. Indeed, the bond order linearly increases with Gb/𝜌b. Table .1 Information on the interactions in Se4.

Atom Interaction No. [a]

VSCC [b] Type Connectivity [c] Gb/𝜌b [d] Cu

Se( 1)

Se(2) Se (2) Cu(1) Cu(2)

Thermal transport

The concept of a strong vibrational mode anharmonicity induced by the presence of Sb lone-pair electrons in Se3, leading to a low thermal conductivity, is generally accepted [START_REF] Tyagi | Thermoelectric properties of Cu3SbSe3 with intrinsically ultralow lattice thermal conductivity[END_REF][START_REF] Skoug | Structural effects on the lattice thermal conductivity of ternary antimony-and bismuth-containing chalcogenide semiconductors[END_REF][START_REF] Zhang | First-principles description of anomalously low lattice thermal conductivity in thermoelectric Cu-Sb-Se ternary semiconductors[END_REF]. This idea was originally proposed by Petrov and Shtrum [START_REF] Petrov | Heat Conductivity and the Chemical Bond in ABX2-Type Compounds[END_REF], stating that the lone-pair 5s electrons could interact with the valence electrons of adjacent atoms upon thermal agitation, causing nonlinear repulsive forces; this nonlinearity of forces manifests itself as an anharmonicity [START_REF] Morelli | Intrinsically Minimal Thermal Conductivity in Cubic I-V-VI2 Semiconductors[END_REF]. The existence of these Sb-5s lone-pair electrons was also evidenced in Se2, indeed, it was shown that the Sb 3+ ions were the main contributors to the averaged Grüneisen parameters [START_REF] Qiu | Diverse lattice dynamics in ternary Cu-Sb-Se compounds[END_REF], and that the very low thermal conductivity and distorted coordination of Se2 mainly arose from the occurrence of the lone-pair electrons of Sb as visualized by ELF [START_REF] Zhang | Ternary CuSbSe2 chalcostibite: facile synthesis, electronic-structure and thermoelectric performance enhancement[END_REF]. As for Se4, all Sb electrons participated to the bonding since each Sb was tetrahedrally coordinated to Se [START_REF] Zhang | First-principles description of anomalously low lattice thermal conductivity in thermoelectric Cu-Sb-Se ternary semiconductors[END_REF], but the valence state of Sb in Se4 is still debated [START_REF] Do | Physics of bandgap formation in Cu-Sb-Se based novel thermoelectrics: the role of Sb valency and Cu d levels[END_REF]. Actually, Skoug et al. [START_REF] Skoug | Role of lone-pair electrons in producing minimum thermal conductivity in nitrogen-group chalcogenide compounds[END_REF] showed that in group V chalcogenides the anharmonicity extent in the crystal lattice is not only affected by the morphology of the lone-pair electron orbit but also by the coordination environment of the group V atom. However, some questions such as how the lone-pair electrons affect the atomic interactions, to what extent weak interactions can be involved in softening vibration modes and to what extent the lattice thermal conductivity can be determined by weak bonds remain unanswered. As can be seen in the following, the topological bonding analysis established in this work can help drive things forward.

The influence of lone-pair electrons on atomic interactions can be clearly seen in Fig. 7. Indeed the Sb lone-pair electrons in Se3 induce two destabilized interactions (Hb>0) dominated by kinetic energy, namely Sb-Se(1) and Sb-Cu(2), leading to the connectivity of five with adjacent atoms. In Se2 the Sb lone-pair electrons induce two stabilized Sb-Se( 2) interactions (Hb<0) dominated by potential energy, and one Sb-Se( 2) and two Sb-Sb destabilized interactions, leading to the connectivity of eight with adjacent atoms. All these interactions have relatively low Ñ 2 rb with a large void between VSCC and boundary of atomic basins (see Fig. 6 and Fig. 8 for Se3 and Se2, respectively).

The atoms involved in the weak interactions play the role of "sink" with negative source contributions [START_REF] Gatti | Source function description of metal-metal bonding in d-block organometallic compounds[END_REF] (see Fig. 9), indicating that these interactions are very strongly non-local ones with charge depletion. Moreover for a specific interaction involving the same atoms, when its characteristics at BCP is close to the (|Vb|/Gb=1, Hb/rb=0) point, the non-local degree of this interaction is high. That means that for Hb>0, the non-local degree of the interaction increases with its strength. This behaviour is opposite to that of interactions for Hb<0 (see Fig. 3). In this case a more local character of bonding is also observed when the distance to the point (1,0) increases but the corresponding bond strength increases. In Se2, the large contribution of the Sb atoms to the anharmonic vibrations [START_REF] Qiu | Diverse lattice dynamics in ternary Cu-Sb-Se compounds[END_REF] can be related to the non-local character of the weak bonds induced by the Sb-5s lone-pair electrons. Indeed the source contributions for the Sb-Sb Van der Waals interaction [START_REF] Bu | Observation of High Seebeck Coefficient and Low Thermal Conductivity in [SrO]-Intercalated CuSbSe2[END_REF] from the adjacent four groups, namely two

[SbSe3] pyramids and two [CuSe4] tetrahedrons, reach 80% (inset of Fig. 9). This indicates a coupling between the non-local Sb-Sb weak interaction and atoms belonging to the adjacent groups, which confers to the atoms a certain extent of contribution to the anharmonic vibrations. The low atomic displacement parameter (ADP) of Sb in Se2 in comparison to that of Sb in Se3 [START_REF] Qiu | Diverse lattice dynamics in ternary Cu-Sb-Se compounds[END_REF] can be explained as shown in Fig. 7 by the distribution of lone-pair electrons in these compounds which also agrees with the fact that lone-pair electrons in Se3 is farther removed from the Sb nucleus [START_REF] Skoug | Role of lone-pair electrons in producing minimum thermal conductivity in nitrogen-group chalcogenide compounds[END_REF]. This dispersive distribution of lone-pair electrons in Se3 weakens both its influence on adjacent atoms and the non-local character of interatomic interactions. By contrast, in Se2 the stabilized Sb-Se interactions induced by the lone-pair electrons confine interactions between Sb and other adjacent atoms. The unit of interatomic distance is Å.

In addition to the weak bonds induced by Sb lone-pair electrons in Se3, other kinds of weak bonds also exist that affect the anharmonic vibrations. The out of plane motion of Cu(1) atoms induces fluctuations in the Cu(1)-Se [START_REF] Kim | Strategies for engineering phonon transport in thermoelectrics[END_REF] interatomic distance that goes from 2.718 Å. to 2.873 Å. (Fig. 10). This Cu(1)-Se(2) interaction is a relatively weak bond, since the source contributions from directly bonded atoms involved in this interaction are 44.5% even with the shortest distance. In fact, the most important anharmonic vibrations in Se3 do not originate from Sb atoms but from Cu(1) atoms, the ADP of Cu(1) atoms being much larger than those for other species [START_REF] Qiu | Diverse lattice dynamics in ternary Cu-Sb-Se compounds[END_REF]. This feature is mainly attributed to the very non-local Cu(1)-Cu(1) interaction almost located at point (1,0) in Fig. 9. However, the unusual temperature dependence behavior and relatively small Grüneisen parameters of all the projected acoustic phonon branches [START_REF] Qiu | Diverse lattice dynamics in ternary Cu-Sb-Se compounds[END_REF] indicate that intrinsic nonlinear phonon scattering is not sufficient to reach the minimal lattice thermal conductivity. Beyond the phonon-phonon anharmonic scattering, the low lattice thermal conductivity is related to phonon scattering covering relatively wide frequency range due to hierarchical chemical bonds distribution in a set of atoms coupled with each other [START_REF] Qiu | Part-crystalline part-liquid state and rattling-like thermal damping in materials with chemical-bond hierarchy[END_REF]. This collective low-frequency vibration from a set of atoms is usually observed in amorphous or liquid materials, which originates from a certain group of atoms within relatively loosely packed regions that undergo a liquid-like motion from one configuration to another [START_REF] Wang | Roles of minor additions in formation and properties of bulk metallic glasses[END_REF]. Qiu et al. [START_REF] Qiu | Part-crystalline part-liquid state and rattling-like thermal damping in materials with chemical-bond hierarchy[END_REF][START_REF] Qiu | Diverse lattice dynamics in ternary Cu-Sb-Se compounds[END_REF][START_REF] Qiu | Phonon" scattering beyond perturbation theory[END_REF] described Se3 as a part-crystalline part-liquid material with high atomic-level heterogeneity, which agrees with the above analysis. The presence of weak and non-local chemical bonding coupled with other adjacent atoms or groups of atoms (the greater the non locality of the bonding the larger the number of involved atoms) substantially widens the frequency range of phonons scattering, therefore leading to low thermal conductivity. Since in Se4 there is no weak interaction (Tab.1), the phonon-phonon scattering only occurs in a narrow range of frequencies leading to high lattice thermal conductivity. For Se3, there are two kinds of weak interactions, one originating from the Sb lonepair electrons (Cu(2)-Sb and Sb-Se(1)), another coming from the motion of Cu( 1) atoms (Cu(1)-Se(2) and Cu(1)-Cu( 1)) (Tab.2).

Since these weak interactions involve all atoms in Se3, they lead to a wide range of low frequency vibrations in the crystal and thus to abnormally low lattice thermal conductivity. As for Se2, all weak interactions originate from Sb lone-pair electrons, namely Sb-Sb and Sb-Se(2) (Tab.3), and the collective low-frequency vibrations mainly arise from the groups of atoms surrounding the Sb-Sb interactions. This situation can explain a higher lattice thermal conductivity for Se2 than for Se3.

Weak bondings (Hb>0) are also observed in Se13. The location of these bonds in Fig. 9 suggests a low thermal conductivity, though higher than that of Se2 and Se3 due to a lower non-local character of these bonds. This low thermal conductivity has been predicted in Ref. [START_REF] Zhang | Prediction of new stable compounds and promising thermoelectrics in the Cu-Sb-Se system[END_REF] although it has been found lower than that of Se3. The reason could be that these predictions account only for the lattice contribution to the thermal conductivity, whereas the electronic one may be substantial in semi-metal such as Se13.

Seebeck coefficient

It is reported for both Se2 and Se4 that the ionic bonding between Cu and Se atoms provide substantial DOS variation near the VBM, showing a large slope at the band edges, which leads to large Seebeck coefficient [START_REF] Yang | New stannite-like p-type thermoelectric material Cu3SbSe4[END_REF][START_REF] Zhang | Ternary CuSbSe2 chalcostibite: facile synthesis, electronic-structure and thermoelectric performance enhancement[END_REF].

Furthermore, by intercalating an ionic layer, consisting of alkaline earth metals and oxygen ions, into the VDW gap of Se2, a significant enhancement of the Seebeck coefficient has been fulfilled [START_REF] Skoug | Structural effects on the lattice thermal conductivity of ternary antimony-and bismuth-containing chalcogenide semiconductors[END_REF]. Hence, the Seebeck coefficient has a very close-knit relationship with ionic bonding.

As seen in Tables 1, 2 and 3, all Cu-Se bonds are ionic with overlapped VSCCs in Se4 and Se2, whereas only part of Cu(1)-Se bonds in Se3 have the same features, the other ones being ionic separated or metallic.

Regarding the Sb-Se bonds, they are all ionic with separated VSCCs in Se4, separated or single VSCCs in Se3 and Se2, the number of single VSCCs being larger in Se2. In spite of the same kind of Cu-Se bonds, the Seebeck coefficient of Se2 is higher than that of Se4, hence the difference could be due in this case to the Sb-Se bonds. It seems that the highest Seebeck coefficient is obtained when VSCC difference is well marked between Cu-Se and Sb-Se bonds, namely in Se2 where all the Cu-Se bonds are with overlapped VSCC and part of Sb-Se bonds are with single VSCC. This assumption is in agreement with the evolution of the Seebeck coefficient, which increases from Se3 to Se4 and Se2.

Conclusions

The structure-property relationships of materials can be analysed in the light of the topological analysis of the chemical bonding, which is sensitive to structural perturbation at atomic level. In this paper, this approach has been used to analyse electronic and thermal transport properties in thermoelectric materials. The characteristics of local energy densities for specific interactions, corresponding to the evolution of the Laplacian distribution of the electron density, are shown to be related to electrical and thermal properties of materials. By localizing the electronic transport path the dispersive and non-local ionic bonds severely affect the electrical conductivity. The asset of this approach is to directly highlight the influence of lone pair electrons and atomic fluctuations on interatomic interactions as well as the coupling effect of weak bonds. The non-locality of weak bonds determinates the extend of involved atoms, directly related to the collective frequency range of phonon scattering, which is the key to the lattice thermal conductivity of materials. This general approach based on chemical bonding analysis from the electron density is a relatively low-cost calculation one compared to that of band structure and phonon calculations, which can be easily performed as a preliminary step to design and search for materials with specific properties.

Fig. 1 (

 1 Fig.1(a)), Se4 bears one Wyckoff position for each Sb and Se atom and two ones for Cu atom. However, Cu(1)-Se and Cu(2)-Se bonds own the same distribution of Ñ 2 r(r) reported in Fig.2(a) and have almost the same bonding characteristics, as can be seen from their very close locations in Fig.3. Fig.2(a) shows that the highly delocalized valence shell charge concentration (VSCC) of Cu overlaps with the localized VSCC of Se along the bond path (BP) and this charge accumulation lies in the Se atomic basin indicating that the Cu-Se bond is a kind of ionic one.Although Fig.2(b) also shows a charge transfer in the Sb-Se bond, the VSCCs are not overlapping. Therefore the high mobility Cu/Se framework path should originate from the overlapping VSCCs between Cu and Se whereas the small void between the VSCCs of Sb and Se acts as a barrier to transport. As for Se3 and Se2, which have an orthorhombic crystal structure (space group Pnma[START_REF] Pfitzner | Cu3SbSe3: Synthese und Kristallstruktur[END_REF][START_REF] Imamov | Determination of crystal structure of CuSbSe[END_REF]) as seen in Fig.1(b) and Fig.1(c), respectively, the interatomic interactions in both compounds are more complex and diverse as shown in Figure3. By contrast to the Sb-Se interactions (see Figure4), the Cu-Se interactions characteristics Hb/rb and |Vb|/Gb are not linearly dependent. The deviation from linearity originates from the severe deformation of the charge distribution between Cu and Se corresponding to the overlapped shell-structure. As exhibited in Fig.3, three kinds of Cu-Se bond exist in Se2 and a comparison between the interatomic distances of the Cu-Se bonds show that the shorter the interatomic distance the farther away from the point (|Vb|/Gb=1, Hb/rb=0). Indeed for the bond labelled E in Fig.3, the interatomic distance is 0.03 Å larger than that of the bond labelled F. This observation is in good agreement with the fact that the ratio between the potential energy and the kinetic energy densities increases as the distance on the Hb/rb vs. |Vb|/Gb graph between the bond of interest and the point (|Vb|/Gb=1, Hb/rb=0) increases. This is confirmed by the calculation of the source functions that shows a larger charge sharing magnitude as the distance to the point (1,0) increases, indicating a more local character of bonding, and a gradual VSCCs overlapping from

Figure 1 .

 1 Figure 1. Crystal structures of (a) Se4, (b) Se3, (c) Se2, (d) Se5 and (e) Se13. Blue, brown and green balls stand for Cu, Sb and Se atoms, respectively.

Figure 2 .

 2 Figure 2. Electron density Laplacian of (a) Cu-Se; (b) Sb-Se in Se4. Green and blue contour curves stand for positive and negative values, respectively.

Figure 3 .

 3 Figure 3. |Vb|/Gb vs Hb/rb for Cu-Se bonding of Cu-Sb-Se ternary compounds. Labels A, B, C, D, E and F stand for the different Cu-Se bonds in Se3 and Se2; the related percentage is the source contributions from directly bonded two atoms. R 2 is the regression coefficient of fitted line and a1 is its slope.

  2 and Tab.3, respectively. The former unit, which affords a conduction pathway, is composed by [CuSe4], [Cu(1)Se3] and [CuSe4], for Se4, Se3 and Se2 respectively, whereas the latter unit, which mostly acts as a barrier to transport, is composed by [SbSe4], [Cu(2)Se3] and [SbSe3], and [SbSe3], respectively. In Se4 the conduction pathway extends the whole crystal by means of the [CuSe4] tetrahedral units (Fig.1(a)), whereas in Se2 due to the electron confinement in the [SbSe3] units (see Fig.1(c)), the electronic transport is localized. As for Se3, the nearly planar [Cu(1)Se3] units are surrounded by [Cu(2)Se3] and [SbSe3] ones (see Fig.1(b)), leading to an even more localized transport path. These building units can explain the decrease of the electrical conductivity from Se4 to Se2 and to Se3, due to the increase in the dispersion and the delocalization of interactions, even if Cu(1)-Cu(2) metal-metal bonds are present in Se3 as shown in Fig.6(a).

Figure 4 .

 4 Figure 4. |Vb|/Gb vs Hb/rb for Sb-Se bonding in compounds of the Cu-Sb-Se system.

Figure 5 .

 5 Figure 5. Electron density Laplacian of different Cu-Se bonds in Se3. (a), (b), (c) and (d) correspond to bonds labelled A, B, C, and D in Fig.3, respectively.

Figure 6 .

 6 Figure 6. Electron density Laplacian of various interactions in Se3. (a) Cu-Cu metal bond; (b) Sb-Se weak interaction; (c) Sb-Cu weak interaction; (d) Cu-Cu weak interaction.

Figure. 7

 7 Figure.7 Connectivity of Sb in (a) Se3 and (b) Se2. Pink balls are BCPs.

Figure 8 .

 8 Figure 8. Electron density Laplacian of various interactions induced by lone-pair electrons in Se2. (a) Sb-Se stabilized weak interaction; (b) Sb-Se weak interaction; (c) Sb-Sb weak interaction.

Figure 9 .

 9 Figure 9. |Vb|/Gb vs Hb/rb for all destabilized interactions of Cu-Sb-Se ternary compounds. Data in red and are the source contributions from directly bonded two atoms and interatomic distances (Å) of the corresponding interactions, respectively. The insert shows the total source contributions from adjacent groups for the Sb-Sb interaction in Se2.

Figure 10 .

 10 Figure 10. Vibrational motion of Cu(1) atoms in Se3. SF% is the source contributions from directly bonded atoms.

Table . 2

 . Information on the interactions in Se3

	(1)	Se	4	overlapped	ionic	4	0.70
	Cu(2)	Se	4	overlapped	ionic	4	0.71
	Sb	Se	4	slightly separated	ionic	4	0.54
	[a] Number of interactions with adjacent atoms . [b] Main characteristic of the electron density Laplacian of specific
	interaction. [c] Total number of interactions for one inequivalent atom. [d] Gb/𝜌b corresponds to the kinetic energy
	per electron (Hartree/e).
	Atom Interaction No.	VSCC	Type	Connectivity	Gb/𝜌b
	Cu

Acknowledgements

The authors are grateful to Pr. A. Otero de la Roza for his advices regarding technical aspects of Critic 2 program and its usage and to Pr. Tantardini for offering the C2 code implementing the source functions and the advices for its usage.

This work 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 A0050806881 made by the Grand Equipement National de Calcul Intensif (GENCI). It was also granted access to the HPC resources of Aix-Marseille Université 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.


Hailong Yang

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

Ziya Aliev

Pascal Boulet

Stability investigation of Se-and Te-substituted tetrahedrite

Keywords: Cu-Sb-S system, substituted tetrahedrite, phase stability, DFT, QTAIM

DFT (density functional theory) and QTAIM (quantum theory of atoms in molecules) calculations have been performed to investigate the influence of Se and Te substitutions on the stability of tetrahedrite. The formation energies, the laplacian of the electron density, the total and local volume changes and the site charges ones have been determined. The DFT results show that the single atom substitution is the most preferred one and that selenium is preferred for sulphur replacement whereas antimony can be equally replaced by selenium and tellurium. However as selenium substitutes for either sulphur or antimony, QTAIM calculations evidenced large structural and chemical changes in the strutures. Additional experimental results, in agreement with QTAIM ones, allow us to conclude that the investigation of alloys stability can not be limited to DFT calculations of formation energies and total volume changes. QTAIM appears as a valuable complementary method to get information on phase stabilities.

INTRODUCTION

Copper antimony chalcogenides have high potentialities as semiconductors for electrooptic devices, thermoelectric devices, optical recording media and sustainable solar cell. As examples we can mention Cu3SbSe3 which has enormously drawn attention of the researchers in the thermoelectric field because of its ultra-low lattice thermal conductivity, and the Cu-Sb-S compounds that are modern state-of-the-art thermoelectric materials for midtemperature heat-to-power conversion with the advantage of being environmental friendly. The thermoelectric figure of merit zT for many doped and natural tetrahedrite materials (Cu12Sb4S13) gets close to unity at 450°C ; this value is comparable to that of conventional ptype TE materials. However because of the complexity of the Cu-Sb-S phase diagram around the Cu12Sb4S13 composition, these materials are always constituted of several phases for which the proportions generally differ from a sample to another one. The large scale use of tetrahedrite for thermoelectric applications necessitates the mastering of the properties and thus the mastering of the material constitution. A way could be to extend the composition range of the tetrahedrite stability. To reach this goal, we investigated the stability of various substituted tetrahedrites in which S atoms were replaced by Se and/or Te ones. This work was mainly carried out by coupled DFT and QTAIM calculations and complemented by some experiments.

THE Cu-Sb-S SYSTEM

Information on the constitutive binary systems Sb-S, Cu-S and Cu-Sb can be found in Ref. 1-6, Ref. 7-8 and Ref. 9-11, respectively. The Sb-S phase diagram only shows one condensed stable binary phase Sb2S3 (stibnite) with a congruent melting and two regions of liquid immiscibility. The Cu-S binary system contains three intermediate phases, namely Cu2S (chalcocite), Cu9S5 (digenite) and CuS (covellite) and a solid solution between digenite and chalcocite extending from 83 to 1129°C. Numerous phases have also been identified at low temperatures (<103°C) but considerable uncertainties still remain both on the stability range of the digenite and on the temperatures of phase transitions in the low temperature range. Moreover no studies on the stability of the new minerals as roxbyite (Cu1.78S), geerite (Cu1.6S) and spionkopite (Cu1.4S) have been performed. There are several intermediate phases in the Cu-rich region but only two (d and h) exist at room temperature.

Phase stability and phase relations in the Cu-Sb-S ternary system have been studied by numerous authors. [START_REF] Skinner | Studies of the sulfosalts of copper III.Phases and phase relations in the system Cu-Sb-S[END_REF][START_REF] Ross | The formation of intermediate sulfide phases in the solid state[END_REF][START_REF] Wernick | New semiconducting ternary compounds[END_REF][START_REF] Cambi | Processi Idrotermali-sintesi di solfosali da ossidi di metalli e metalloidi 2 Cuprosolfoantimoniti[END_REF][START_REF] Kuliev | Phase Equilibrium and Intermolecular Interactions in Systems Formed by Copper and Antimony Chalcogenides[END_REF][START_REF] Kuliev | Phase equilibria between copper and antimony chalcogenides[END_REF][START_REF] Godovikov | Bismuth Sulfosalts : their chemical composition, synthesis and classification[END_REF][START_REF] Il'yasheva | Investigation of the system Cu2S-Sb2S3 at 320-400[END_REF][START_REF] Bryndzia | Liquidus phase relations on the quasi-binary join Cu2S-Sb2S3: Implications for the formation of tetrahedrite and skinnerite[END_REF] There are four compounds in the Cu-Sb-S ternary system, namely CuSbS2 (chalcostibite), Cu3SbS3 (skinnerite), Cu3SbS4 (famatinite), and Cu12Sb4S13 (tetrahedrite). Chalcostibite, CuSbS2, forms readily by reaction from the elements and did not show any composition range. [START_REF] Skinner | Studies of the sulfosalts of copper III.Phases and phase relations in the system Cu-Sb-S[END_REF] By contrast, Golovei et al. [START_REF] Golovei | Phase diagram of the system Cu2S-Sb2S3 in the region of existence of CuSbS2[END_REF] found a homogeneity range for CuSbS2 and a polymorphic transformation. Skinnerite was first reported by Godolikov et al. [START_REF] Godovikov | Bismuth Sulfosalts : their chemical composition, synthesis and classification[END_REF] who prepared it by sintering Cu2S and Sb2S3 at 480°C and 530°C. It has the ideal composition of Cu3SbS3 with a small composition field along the line Cu2S-Sb2S3 at 500°C. Karup-Møller et al. [START_REF] Karup-Møller | Skinnerite Cu3SbS3, a new sulfosalt from the IIimaussaq alkaline intrusion, South Greenland[END_REF] have found these ore aggregates crystallized at about 400°C from metal-rich residual magmatic fluids markedly deficient in sulfur. Cu3SbS3 is stable in the Cu-Sb-S system in the temperature range 359-607°C. Below 359°C this phase decomposes into a mixture of native antimony, chalcostibite and tetrahedrite, 3 but if it is cooled slowly, it decomposes to give a mixture of tetrahedrite, antimony and digenite. [START_REF] Whitfield | Polymorphism in skinnerite Cu3SbS3[END_REF] At 607°C, Cu3SbS3 congruently melts.

Pfitzner [START_REF] Pfitzner | Disorder of Cu + in Cu3SbS3:structural investigations of the high-and lowtemperature modification[END_REF] have investigated both the high and low-temperature phases. The high-temperature modification is disordered and stable above 122°C, and is supposed to show high copper ion conductivity or at least a high copper diffusion coefficient. The low-temperature form is ordered and stable below -9°C. Famatinite, Cu3SbS4, derived from sphalerite, melts congruently at 627°C, and is not known to have any stable polymorph forms. [START_REF] Skinner | Studies of the sulfosalts of copper III.Phases and phase relations in the system Cu-Sb-S[END_REF] The only compound with an extensive composition field in the Cu-Sb-S system is tetrahedrite. The basic structure of minerals of the tetrahedrite group was solved by Pauling and Newman 24 and now the widely accepted formula Cu12Sb4S13 represented its composition. As for the Pauling and Neuman formula, Cu12+xSb4+yS13, where 0≤x≤1.92 and -0.02 ≤y≤ 0.27, there is no clear evidence so far proving the substitution schemes implied in this formula are correct. [START_REF] Pauling | The crystal structure of binnite, (Cu, Fe)12As4S13, and the chemical composition and structure of minerals of the tetrahedrite group[END_REF] Cambi and Elli 14 have given an extreme composition range for tetrahedrites going from 4Cu2S•Sb2S3 to 3Cu2S•Sb2S3 (i.e. Cu8Sb2S7 to Cu3SbS3). Since Skinner et al. [START_REF] Skinner | Studies of the sulfosalts of copper III.Phases and phase relations in the system Cu-Sb-S[END_REF] have not managed to prepare a single homogenous phase with the composition Cu12Sb4S13, they believed that this phase has a stability temperature range below 300°C. These authors also found that the unit cell edges of tetrahedrite decrease with the sulfur content. Tatsuka and Morimoto 25 reported the formation of pseudo-tetrahedrite (T'3) which is stable between 250 and 350°C and is a superlattice of Cu12+xSb4+yS13 (T3) (lattice parameter of T'3 is twice that of T3). Along the Cu2S-Sb2S3 section reported by Chang et al., [START_REF] Chang | Phase Diagrams and Thermodynamic Properties of Ternary Copper-Sulfur-Metal Systems, Incra Monograph VII[END_REF] the composition of pseudo-tetrahedrite ranges from Cu3SbS3 to Cu12.39Sb4.54S13. Below 250°C, a prolonged annealing of T3 or T'3 results in a phase mixture consisting of Cu3SbS4, Cu2-xS and Sb. Makowicky and Skinner [START_REF] Makovicky | Studies of the sulfosalts of copper VII. Crystal structures of the exsolution products Cu12.3Sb4S13 and Cu13.8Sb4S13 of unsubstituted synthetic tetrahedrite[END_REF] have confirmed that tetrahedrite has two room-temperature phases Cu-poor phase Cu12.3Sb4S13 and Cu-rich phase Cu13.8Sb4S13 produced by a low-temperature exsolution of tetrahedrite Cu12.59Sb4.03S13. However when melts, whose compositions lie near that of tetrahedrite, are cooled at the rate of 50-100°C/h, instead of separation of the phase Cu3SbS3, metastable solid tetrahedrite solutions crystallize out and then decompose to phases with a complex structure, which are absent in the phase diagram. [START_REF] Il'yasheva | Peculiarities of tetrahedrite crystallization in the system Cu-Sb-S[END_REF] The crystal structure of stoichiometric tetrahedrite Cu12Sb4S13 was determined for the first time by Skinner et al. [START_REF] Skinner | Studies of the sulfosalts of copper III.Phases and phase relations in the system Cu-Sb-S[END_REF] There is no phase transitions for Cu12Sb4S13 up to 300°C, and the two exsolution phases Cu12.3Sb4S13 and Cu13.8Sb4S13 show a high copper-ion mobility even at room temperature. According to [29] tetrahedrite is stable from 4Cu2S•Sb2S3 to Cu3SbS3 at 300°C and from 4Cu2S•Sb2S3 to Cu12Sb4S13 at 400°C. Whereas Skinner et al. [START_REF] Skinner | Studies of the sulfosalts of copper III.Phases and phase relations in the system Cu-Sb-S[END_REF] reported above 500°C a rapid shrinkage of the composition range for tetrahedrite followed by its decomposition into Cu3Sb1S4, Cu3SbS3 and CuxS (x going from 1.6 to 2) at 543°C, Il'yasheva 28 evidenced an incongruent melting for this phase at 574±2°C. All the ternary phases are stable at 500°C. The crystallographic data and the thermodynamic properties of the ternary solid phases are given in Table 1 andTable 2

Computational details

The calculations have been performed in the frame of the density-functional theory [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF][START_REF] Kohn | Self-Consistent Equations Including Exchange and Correlation Effects[END_REF] using the gradient-corrected Perdew-Burke-Ernzerhof (PBE) exchange-correlation functionals. [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] All the structures have been relaxed to their equilibrium state and the energy and force thresholds have been fixed to 10 -7 Ry and 10 -4 Ry.bohr -1 , respectively. The projector augmented wave (PAW) method 38 has been used to treat nucleus and core and valence electrons interactions with cutoffs 48 Ry for the energy and 480 Ry for the density. A 4x4x4 k-point mesh has been utilized. The Quantum Espresso program package has been used for these calculations. [START_REF] Giannozzi | Quantum Espresso : a modular and opensource software project for quantum simulations of materials[END_REF] The analysis of the bonding structure of the alloys has been performed using Bader's quantum theory of atoms in molecule [START_REF] Bader | Atoms in Molecules -A Quantum Theory[END_REF] as implemented in the Critic2 program. 41

Crystalline description of the tetrahedrite

The tetrahedrite crystal structure with stoichiometric formula Cu12Sb4S13 is composed The formation energies have been determined from the energies E calculated by DFT according to the following formula:

Δ𝐸 ! = "#Cu x Sb ! S z X t $%x"(Cu)%y"(Sb) -z"(S)%t"(X) x)y)z)t (1)
Recalling that the number of formula in the unit cell of tetrahedrite is 2, a one-atom sulphur and antimony substituted tetrahedrite on the one hand and a two-atoms sulphur and antimony substituted tetrahedrite on the other hand correspond to Cu12Sb4S12.5X0.5, Cu12Sb3.5X0.5S13, and Cu12Sb4S12X, Cu12Sb3XS13, respectively.

Results

Formation energies and total lattice volume changes

Substitutions of Se and Te for S. The calculated values are given in Table 3. For a single substitution, both selenium and tellurium atoms prefer occupying site 24g, the formation energies being -160.7 and -151.6 meV/atom respectively. The substitution at site 2a is much less favourable by about 8 meV/atom. This preference for the 24g site has already been experimentally evidenced for Se by Lu et al. [START_REF] Lu | Phase Stability, Crystal Structure, and Thermoelectric Properties of Cu12Sb4S13-xSex Solid Solutions[END_REF] For two Te or Se atoms substitutions as mentioned above three possible situations have been envisaged: two atoms at site S 24g, one atom at site S 24g and the other at site S 2a, and 2 atoms at site S 2a. For the two former possibilities we observe that for a given substituting atom the energies are very close and actually lie in the error bar that we estimate to be 3-4 meV/atom. The 2-atoms substitution at site 2a is the least preferential one for both Se and Te. We note that, the substitution with selenium is more preferred than that with tellurium, as the formation reactions pertaining to Se are more exothermic than those pertaining to Te.

These results show that whatever the substituting atom (Se or Te), the most favorable situation corresponds to a single atom substitution at site 24g and that irrespective of the considered site, the replacement of S by Se is the most favorable one. The total volume of the substituted structures changes sizeably with respect to that of pristine tetrahedrite as can be seen in Table 3. We observe that, except for 2-Se-atoms substitution, the more energetically favorable the site is, the smaller the volume change. In agreement with the respective atomic volumes of Se and Te atoms, sites substituted with Te are likely to be more affected, as we observe a larger total volume change, than when Se is the substituting element. It appears that, substitutions on site S 2a induces more drastic changes in the total volume than those on site S 24g, the volume change being about two to three times as large. As a consequence, the volume increase for 2-atoms substitution at site S 2a is as large as +6 and nearly +2 Bohr 3 for Te and Se, respectively. These results on the volumes are in line with the energetic ones as they can lead to the same inferences as above.

Substitutions of Se and Te for Sb. The calculated values are given in Table 4.

Remarkably enough, the formation energies are quite similar, suggesting therefore that both Regarding the total volume change of the lattice upon substitutions, it is remarkable that they are all negative irrespective of the substituting atom and configurations. The smallest volume change corresponds to the one-atom substitution. For two-atoms substitutions the situation is contrasted between Se and Te. For Se, the volume changes remain rather small (less than 2 Bohr 3 ). For Te, the two-atoms substitutions at the edge of the pseudo-cube still lead to the modest volume change of -2.6 Bohr 3 . By contrast, a huge change is observed when two Sb atoms are replaced by Te on the pseudo-cube face and body diagonals.

Notwithstanding this observation, the structures are not energetically destabilized as mentioned above. The lattice volume change is positive when substitutions occur at the sulphur sites, whereas it is negative at the antimony ones, the one atom substitution leading in all cases to the lower volume variation. For two-atoms substitutions, irrespective to the substituted element and site, Te leads to higher volume changes than Se. One can note that, whereas these volume variations are extremely large when Te substitutes for Sb, the energetic variations are quite small.

Local volume and atomic charge changes and electron density laplacian

In order to get better insight into the influence of substitutional atoms in the tetrahedrite, the atomic site volume and charge changes and electron density laplacian have been calculated using Bader's theory of atoms in molecules. Only one-atom substitutions are considered here since these lead to the structures bearing highest formation energies and smallest volume changes.

Pure tetrahedrite. The electron density laplacian of pure tetrahedrite in the (110) plane is depicted in figure 2. This plane contains S 2a atoms (e.g. sulphur located at abscissa ca. 22), which are bonded to two copper 12e (at abscissa ca. 20 and 24); the four remaining copper 12e atoms that complete the octahedra around sulphur ones are located in two planes rotated by 45° from the (110) one and are thus not visible.

In the (110) plane, one can see that in agreement with the structure description already reported in literature, Cu 12e is bonded to two S 24g atoms and despite the long distance between copper 12e and antimony atoms (3.536 Å, see above), these atoms are chemically bonded to each other since a bond path is found to connect them and a bond critical point (BCP) is found at the atomic basins boundary. Indeed, in the Bader's description the bonds are characterized by a path connecting two atoms, with b (see Fig. 2) depicting the position of the BCP where the electron density gradient vanishes. The effect of this replacement on the atomic charges and atomic volumes can be seen in Table 5. Interestingly, sulphur 2a is the most affected atom with a positive change of both properties; a volume increase as large as +5 bohr 3 is observed. The Cu 12e atom to which Se is bonded is also affected, though in a lower extent and in the opposite way with a decrease of both the volume and the charge. The site of the antimony atom to which Se is bonded undergoes a volume decrease although its charge remains unchanged. The volume decrease is about the same in absolute value as the volume increase of the S 24g site where the substitution occurs.

As Te substitutes for sulphur 24g the same overall trends are observed as for Se. The graph of the electron density laplacian (Fig. 3b) shows that one of the Cu12e-Sb is broken (3.592 Å, +1.2%) while the other one is strengthened (3.181 Å, -10%) and the octahedron around S 2a is distorted (the Cu-S-Cu angle is 163°). However, we note that these structural changes are less drastic than for Se, although the Te substitution is energetically less favorable than the Se one. The volume change of the site 24g of sulphur, where the Te substitution occurs, increases sizably (Table 6) by +3.3 Bohr 3 (2.9 times that for Se substitution) and the volume of the site 12d (Cu) is almost twice as large as that for the Se substitution. Only the site 12e of Cu undergoes a volume decrease of about the same amount as for Se substitution. These local volume changes explain why the overall volume of the structure increases (Table 3) as Te substitutes for sulphur.

We note finally that the atomic site charges vary in the same way and same extent as for the Se substitution.

Substitutions of Se and Te for Sb. The influence of the Se and Te substitutions for

Sb on the electron density laplacian and bonding network of the tetrahedrite is depicted in figure 4a and 4b, respectively. Interestingly, as Se substitutes Sb a new chemical bond appears between the selenium atom and the sulphur one of the 2a site in spite of the large interatomic distance (4.027 Å); The bond is characterized by a bond path between both atoms and a bond critical point located ca. x=24 (Figure 4a). The Se-S2a interatomic distance is 1.156 Å shorter than the Sb-S2a one of the pure tetrahedrite. The appearance of this new bond implies a distortion of the octahedron around the S2a atom; the Cu-S-Cu angle amounts to 153.07° and the sulphur atom moves out of the Cu4 plane by 17.9° in the direction of Se. As to the case of the Te for Sb substitution, the bonding network of the structure is not affected compared to that of the pure tetrahedrite (Figure 4b). Some differences in the atomic charge and volume changes between these two types of substitutions can be observed (Tables 7 and8). First, the atomic charges are barely varying upon Te substitution. By contrast, S2a and Sb undergo significant increase and decrease, respectively, of their charge upon Se substitution, which can be correlated to the appearance of the new bond between these atoms.

The volumes are affected in various ways. For both Se and Te substitutions, the volume of S2a increases drastically by up to 5 and 7 Bohr 3 , respectively, and that of Cu12e

(bonded to S2a) decreases by 2.5 Bohr 3 , on average. The Cu12d atoms being far from the substitution site, their volumes are not affected. As to the substituted site Sb, the volume decreases (increases) drastically (slightly) when Se (Te) is the substituting atom, these variations being in line with the differences in atomic sizes between Sb, Se and Te. Since the substitution at the Sb site is energetically equally plausible by Se and Te, we can infer that the total contributions to the formation energy of the various structural, electronic and bonding changes occuring in the structures are comparable. Partial conclusions. Irrespective of both the substituting and the substituted atoms, the atomic charge changes are quite weak on all sites, the maximum being 0.06 electron charge unit. When substitution occurs on the S24g site, irrespective of the substituting atom, charge change is only noticeable on the S2a site. When substitution occurs on the Sb8c site, the largest charge change occurs on this site while being weaker as Te substitutes for Sb.

Regarding the local volume changes, the S2a site is the most impacted by both S24g and Sb replacements. In all cases, the volume of the site S2a increases. One can see that the volume of the site Cu12e, which is adjacent to the site S2a, decreases, though to a lesser extent. As one could expect, the impact of the substitution on the local volume changes is weaker when the site of interest is far from the substituted one (e.g. volume of the Cu12d site when the replacement is on the Sb site).

The electron density laplacian shows few or even no chemical and structural modifications when Te substitutes for S24g or Sb, respectively. In the case of Se, bigger structural distortions have been evidenced with the breaking of a Cu12e-Sb bond when the substitution occurs on S24g site and the formation of a Se-S2a one when the substitution occurs on Sb site.

SUBSTITUTED TETRAHEDRITE: EXPERIMENTAL PART

DFT calculations of formation energies and total volume changes showed that : i) the most favourable situation corresponds to a single atom substitution, ii) selenium is preferred for the sulphur substitution whereas tellurium and selenium are equally likely to replace antimony, whereas QTAIM calculations evidenced larger structural and chemical distortions when Se is the substituting element irrespective of the substituted site. Since these results seem inconsistent, we experimentally investigated the following compositions:

Cu12Sb3.5Se0.5S13, Cu12Sb4S12.5Te0.5 and Cu12Sb3.5Se0.5S12.5Te0.5. We have not investigated Cu12Sb3.5Te0.5S13 and Cu12Sb4S12.5Se0.5 as these works were already done by Bouyrie et al. [START_REF] Bouyrie | Crystal structure, electronic band structure and high-temperature thermoelectric properties of Te-substituted tetrahedrites Cu12Sb4-xTexS13 (0.5≤x≤2.0)[END_REF] and Lu et al. [START_REF] Lu | Phase Stability, Crystal Structure, and Thermoelectric Properties of Cu12Sb4S13-xSex Solid Solutions[END_REF] , respectively, and the authors found single-phased materials with a tetrahedrite structure.

Synthesis and characterization

The three aforementioned ternary and quaternary alloys have been prepared by mixing their constitutive elements with high purity grade in required proportion in evacuated (~10 - 2 Pa) sealed silica ampoules. The mixture was heated up to 680 °C and maintained at this temperature for 8 h. In order to get the equilibrium state, the samples have been quenched from melt in cold water and annealed at 400°C for about 150 h. The constitutive phases of the samples were evidenced by X-ray powder diffraction using a Bruker D8 Advance diffractometer with Cu-K radiation in the 2θ range from 5° to 80° and their proportion were evaluated by Rietveld refinement.

Results

The X-ray diffraction patterns for the three investigated samples are given in Figure 5 and the constitution of the alloys with the percentage of each phase determined by Rietveld refinement are given in Table 9.

When one Te atom is substituted for a S one, the alloy is constituted of one phase only which crystallizes in a tetrahedrite structure. As expected the lattice parameter of this phase (a = 1036.1 pm) is larger than that of the pure sulphur tetrahedrite (a=1032.4 pm, see above).

When one Se atom is substituted for a Sb one, the alloy is constituted of two phases in an equivalent ratio. One of these phases crystallizes in a tetrahedrite structure with a lattice parameter a= 1034.3 pm, the other one crystallizes in a famatinite structure with the lattice parameters a=539.0 pm and c=1076.1 pm which are larger than those of the pure sulphur famatinite (a=538.5 pm and c=1075.4 pm, see above). When one Se atom is substituted for a Sb one and one Te is substituted for a S one, the alloy is mainly constituted (95%) of a tetrahedrite phase, the secondary phase (5%) crystallizing in a famatinite structure.

The whole results show that in spite of low formation energies and high volume changes, alloys containing a tellurium atom at a sulphur site stabilize the tetrahedrite structure. This statement is in agreement with the QTAIM calculations, which find lower structural and chemical distortions when Te is the substituting element.

CONCLUSION

In this work, we investigated Se and Te substitutions in tetrahedrite by coupling DFT and QTAIM calculations. DFT calculations of formation energies and total volume changes showed that the most favourable situation corresponds to a single atom substitution and that selenium is preferred for the sulphur substitution whereas tellurium and selenium are equally likely to replace antimony. By contrast, QTAIM calculations evidenced irrespective of the substituted site, larger structural and chemical distortions when Se is the substituting element.

In order to clarify these inconsistencies, we experimentally investigated three alloys and found that by contrast to a Se substitution for Sb, a Te susbtitution for S stabilizes the tetrahedrite structure. This investigation shows that a DFT calculation of formation energies and total volume changes is not sufficient to put forward conclusions on alloys stability. By giving chemical and structural information on a compound, in addition of being a good method to explain the structure-properties relationships, QTAIM appears as a valuable method to get information on phase stabilities.

Fig. 1 :

 1 Fig. 1: Structural features of the tetrahedrite. a) Environment of copper in 12e Wyckoff

  Se and Te may equally substitute for Sb. The most preferred substitution is by far obtained with one Se or Te atom only in place of antimony. The energy difference between the oneatom and two-atoms substitutions is above 20 meV/atom for both Se and Te, the largest energy difference being 23 meV/atom for Se. For two-atoms substitutions, the three configurations, namely two atoms along the pseudo-cube edge, two atoms along the face diagonal of the pseudo-cube and two atoms along the principal diagonal of the pseudo-cube, are energetically undistinguishable since the differences in formation energies lie in the error bar (3-4 meV/atom) ; they spread in the ranges [-122.2; -123.8] meV/atom and [-120.4; -122.5] meV/atom for Te and Se, respectively.

  From the formation energies viewpoint the most favourable situation corresponds to a single atom substitution, both Te and Se preferring replacing sulphur at the 24g Wyckoff position to sulphur 2a and antimony 8c ones. Selenium is preferred for the sulphur substitution whereas tellurium and selenium are equally likely to replace antimony.

Fig. 2 : 2 (Fig. 3 :

 223 Fig. 2: Electron density laplacian of pure tetrahedrite. Legend. b: bond critical point; r: ring

Fig. 4 :

 4 Fig. 4: Electron density laplacian of substituted tetrahedrite. a) Se substitution for Sb 8c ; b)

Fig. 5 :

 5 Fig. 5: X-ray diffraction patterns of Cu12Sb4S12.5Te0.5, Cu12Sb3.5Se0.5S13 and

, respectively. Table I. Crystallographic data of the Cu-Sb-S ternary solid phases

	Phase / Name	Crystal	Space	Lattice parameters (pm)	Reference
		system	group
	Cu3SbS4	Tetragonal	I-42m	a=538.5,c=1075.4	[30]
	/Famatinite
	CuSbS2	Orthorhombic Pnma	a=601.8,b=379.58,c=1449.5 [32]
	/Chalcostibite
	Cu12Sb4S13	Cubic	I-43m	a=1032.4	[32]
	/Tetrahedrite
	Cu3SbS3	Monoclinic	P21/c	a=781.5,b=1025.2,c=1327.0 [33]
	/Skinnerite

Table II. Thermodynamic properties of the Cu-Sb-S ternary solid phases from Ref. 34

	SUBSTITUTED TETRAHEDRITE: COMPUTATIONAL PART
	Phase Temperature	ΔfG 0	ΔfH 0	S 0
			(J/mol of atoms)	(J/mol of atoms)	(J/K.mol of atoms)
	Cu3SbS4	25°C	-207500 ± 3900	-200000 ± 6300	298 ± 18
	CuSbS2	25°C	-121400 ± 3400	-119000 ± 3500	148 ± 6.4
	Cu3SbS3	25°C	-197800 ± 3800	-189300 ± 6500	269.5 ± 13.7

Table III : Formation energies (in meV/atom) and total volume change of the cell (in bohr 3 ) for 1 and 2 atoms Te or Se substitutions for sulphur. Columns 1 and 8 correspond to the ranking w.r.t. formation energies.

 III

	Te substituted for S	ΔEf	ΔV	ΔEf	ΔV	Se substituted for S
	1 1 Te for S 24g	-151.6 +1.297 -160.7 +0.093 1 Se for S 24g	1
	2 1 Te for S 2a	-144.0 +2.452 -153.1 +0.312 1 Se for S 2a	4
	3 2 Te for 2 S 24g	-136.7 +3.066 -155.5 +0.943 2 Se for 2 S 24g	2
	4 1 Te for S 24g and	-136.2 +3.976 -153.5 +0.685 1 Se for S 24g and	3
	1 Te for S 2a					1 Se for S 2a
	5 2 Te for 2 S 2a	-127.3 +6.200 -144.7 +1.828 2 Se for 2 S 2a	5

Table IV : Formation energies (in meV/atom) and total volume change of the cell (in bohr 3 ) for 1 and 2 atoms Te or Se substitutions for antimony. Columns 1 and 8 correspond to the ranking w.r.t. formation energies.

 IV

		Te substituted for Sb	ΔEf	ΔV	ΔEf	ΔV	Se substituted for Sb
	1	1 Te for Sb 8c		-144.8 -1.2	-143.4	-1.24	1 Se for Sb 8c	1
	2	2 Te for 2 Sb 8c		-123.8 -2.6	-122.5	-1.91	2 Se for 2 Sb 8c	2
		(pseudo-cube edge)					(pseudo-cube edge)
	3	2 Te for 2 Sb 8c		-123.6 -26.5	-120.7	-1.71	2 Se for 2 Sb 8c	3
		(pseudo-cube						(pseudo-cube
		diagonal)						diagonal)
	4	2 Te for 2 Sb 8c		-122.2 -36.1	-120.4	-1.93	2 Se for 2 Sb 8c
		(pseudo-cube	face
		diagonal)

Table V : Variation of the atomic charge (in electron charge unit) and atomic sites volume (in bohr 3 ) upon Se substitution for S 24g.

 V

		24g (S)	2a (S)	12e (Cu)	12d (Cu)	8c (Sb)
	Dq	+0.01	+0.05	-0.01	-0.01	0
	DV	+1.15	+5	-2.26	+0.45	-1.31

Table VI : Variation of the atomic charge (in electron charge unit) and atomic sites volume (in bohr 3 ) upon Te substitution for S 24g

 VI

		24g (S)	2a (S)	12e (Cu)	12d (Cu)	8c (Sb)
	Dq	+0.01	+0.04	-0.01	-0.01	-0.01
	DV	+3.3	+5.6	-2.1	+0.7	+0.2

Table VII : Variation of the atomic charge (in electron charge unit) and atomic sites volume (in bohr 3 ) upon Se substitution for Sb 8c.

 VII

		24g (S)	2a (S)	12e (Cu)	12d (Cu)	8c (Sb)
	Dq	+0.03	+0.05	-0.01	0	-0.06
	DV	-0.90	+5.1	-2.32	0	-4.07

Table VIII : Variation of the atomic charge (in electron charge unit) and atomic sites volume (in bohr 3 ) upon Te substitution for Sb 8c.

 VIII

		24g (S)	2a (S)	12e (Cu)	12d (Cu)	8c (Sb)
	Dq	+0.01	+0.01	0	+0.01	-0.02
	DV	+0.60	+7.2	-2.7	-0.2	+0.6

ACKNOWLEDGEMENTS

This work is financially supported by China Scholarship Council (CSC N°201608530165). This work was granted access to the HPC resources of the Centre Informatique National de l'Enseignement Supérieur (CINES), Montpellier, France under allocation A0050806881 made by the Grand Equipement National de Calcul Intensif (GENCI). It was also granted access to the HPC resources of Aix-Marseille Université 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.


Wei Liu

Xiao-Gang Lu

Pascal Boulet

Marie-Christine Record

Hao Wang

Qing-Miao Hu

Influence of atomic mixing and atomic order on molar volume of the binary sigma phase

The present work systematically clarifies the influence of the atomic mixing and the atomic order (i.e. atomic constituent distribution or site occupancy preference on inequivalent sites of a crystal) on the molar volume of the binary sigma phase. We calculated the molar volumes of 21 binary sigma compounds with different atomic order by using EMTO-CPA (Exact Muffin-Tin Orbitals -Coherent Potential Approximation) method. The cal-culation results follow a similar trend as the experimental data. Moreover, the calculated molar volumes with respect to the composition of the sigma phase in disordered states show that the sigma compounds do not comply with the linear volume-composition relationship. This is due to the influence of the atomic mixing and the deviation can be either negative or positive depending on the tendency of electron loss or gain of the two constituent elements. Additionally, we compared the molar volume difference (ΔVm = Vdisorder-Vorder) between the hypothetically disordered and ordered states of the binary sigma compounds. The results indicate that the effect of the atomic order on the molar volume of the sigma phase depends on the size difference and the electron configurations of the two constituent elements.

Introduction

The sigma phase is a nonstoichiometric intermetallic compound with tetragonal structure (space group P42/mnm) containing 30 atoms distributed on five inequivalent sites [1][2][3] as shown in Table 1. This phase exists in 49 different binary systems, among which Cr-Ni, Ta-V, Ti-Mn, Zr-Ir and Zr-Re are not confirmed. It is a hard brittle phase mostly formed between transition elements and when presenting in stainless and other high-alloy steels as well as Ni-based superalloys, it deteriorates the properties of products [1,2].

As an important physical property, the molar volume (or lattice parameters) of the sigma phase has attracted great interests of investigations [1][2][3][4][5], by using XRD (X-ray diffraction), ND (neutron diffraction), CALPHAD [6][7][8] or first-principles calculation methods.

From the 1950s-2010s, there are numerous measurements (see e.g. Ref. [5] and references therein) on the lattice parameters of the sigma compounds by using XRD or ND methods. The measured data from the 2000s-2010s are more accurate. Most measured volumes with respect to composition present linear relationship within the homogeneity range [5]. Besides,Crivello et al. [4] calculated the cell volumes of the stoichiometric end-member compounds (i.e. the complete set of ordered configurations) of Re-X (X = Co, Cr, Fe, Hf, Mn, Mo, Nb, Ni, Os, Ru, Ta, Tc, Ti, V, W and Zr) by using first-principles calculations. They indicated that for most systems, the cell volume with respect to composition present an approximately linear relationship, while for Hf-Re, Nb-Re, Ta-Re, Ti-Re and Zr-Re, a negative deviation from linearity was observed (see Fig.5 in Ref. [4]), which indicates a denser packing with atomic mixing [4]. However,Crivello et al. [4] did not further discuss the origin of the deviation.

In a recent work, Liu et al. [5] modeled and evaluated the molar volume of the binary sigma phase at room temperature and atmospheric pressure by using the CALPHAD method and establishing a CALPHAD-like database. The adopted model expresses the molar volume of a nonstoichiometric sigma phase as a linear average of volumes of the constituent elements in their hypothetic sigma structure [5]. The proposed linear model can well reproduce most experimental data. It has interesting features in developing databases, such as reducing the number of assessed parameters and separating the volume data from the corresponding thermodynamic databases. The proposed linear model works well as a phenomenological method, which however does not present any physical and chemical insights.

In fact, two influencing factors were neglected in the above linear model, namely atomic mixing and atomic order (i.e. atomic constituent distribution or site occupancy preference on inequivalent sites of a crystal). Both the atomic mixing and atomic order may lead to a deviation from the linearity mentioned above. Unlike solution phases, the molar volume of the sigma phase can vary with the atomic order even at a fixed composition. The purpose of the present work is to study how the atomic mixing and atomic order affect the molar volume of the sigma phase.

Hall and Algie [1], Joubert [2] and then Liu et al. [9] systematically investigated the atomic order of the binary sigma phase. Liu et al. [9] indicates that between the two constituent elements, the one with larger size, smaller number of valence electrons or smaller total number of electron shells prefers occupying large coordination number (CN) sites and the other one prefers occupying small CN sites.

In the present work, we calculated the molar volume of the A-B binary sigma phase with different atomic order by using EMTO-CPA approach. The atomic order involves the completely ordered and dis-ordered states as defined in Table 1, as well as the experimental site occupancies. We used the molar volume difference between the A2B sigma compounds in the ordered and disordered states to describe the influence of atomic order on molar volume. On the other hand, we describe the influence of atomic mixing on molar volume by studying deviation of the volume of the hypothetically disordered sigma phase from linear volume-composition relationship, which describes volume variations of the hypothetically disordered sigma phase spanning from pure sigma A to pure sigma B. Notably, the disordered state is the closest state to solid solution. Thus, it excludes the influence of atomic order on molar volume.

To facilitate explanation, all the binary sigma compounds investigated are designated as A-B where atom A always has a larger size than atom B (i.e. V(A) > V(B). The size of atom A or B is determined by the molar volume of pure A or B in their hypothetic sigma phase structure from the CALPHAD assessments [5,9]. Besides, we introduced different categories of the sigma phase divided by Liu et al. [9] as shown in Table 2. It is classified according to the electron configurations of the two constituent elements. There exist six different categories, i.e. SMe, SMs-SMe, LAs, LAs-SMe, LAs-LAe, A-Al (A = Nb, Ta), where SMe or LAe indicates that atom A has smaller or larger number of valence electrons than atom B, respectively; SMs or LAs indicates that atom A has smaller or larger total number of electron shells than atom B, respectively.

Methodology and calculation details

First-principles calculations were performed using the exact muffin-tin orbitals (EMTO) method [10,11] based on the density-functional theory (DFT) [12]. With this method, Green's function technique is used to solve the one electron Kohn-Sham equation. The effective potential in the one-electron equation is treated with the optimized overlapping muffin-tin approximation. Besides, the total energy is corrected with the full charge density (FCD) method [10]. The basis sets of the exact muffin-tin orbitals [10,11] are used to expand the wave function. The coherent potential approximation (CPA) [13][14][15] is incorporated within the EMTO code, which facilitates the calculations dealing with chemical disordered alloys. With these techniques, the EMTO-CPA method holds a great higher accuracy than the linear muffin-tin orbital (LMTO) method [16] and ensures the reliability of calculations dealing with atomic and magnetic disordered systems [17].

In the present work, Green's function was calculated for 16 complex energy points distributed exponentially on a semicircular contour. The scalar-relativistic and soft-core approximations were employed. We used the electronic exchange-correlation functional within the generalized-gradient approximation (GGA) parameterized by Perdew et al. [18]. The Brillouin zone was sampled by a uniform k-point mesh without any smearing technique. The k-point mesh was set to 3×3×6. For Cr-Fe and V-Fe, the magnetic parameters for ferro-magnetism were set.

Results and discussion

The binary sigma compounds chosen to perform EMTO-CPA calcu-lations are Cr-Co, Cr-Fe, Cr-Mn, Mo-Co, Mo-Mn, Mo-Os, Mo-Re, Nb-Al, Nb-Ir, Os-Cr, Re-Cr, Re-Mn, Re-V, Ru-Cr, Ta-Al, Ta-Pd, Ta-Rh, V-Co, V-Fe, V-Mn and V-Ni. The molar volumes of the pure elements in a hypothetic sigma phase structure obtained from first-principle calculations are compared in Fig. 1 with those obtained from the CALPHAD assessments. This figure shows that the volumes calculated by EMTO-CPA and VASP methods follow a similar trend. Besides, it indicates that for most elements, the differences between the first-principles calculations and the CALPHAD assessed ones are within 10%. The small deviations validate our calculations. Fig. 2 shows the calculated molar volumes of the A-B (V(A) > V(B)) binary sigma compounds compared with the experimental data. The calculated molar volume data with respect to the composition follow similar changing trends as compared with the measured ones, espe-cially for Cr-Co, Mo-Re, Re-Cr and Ta-Al. In addition, Fig. 2 indicates that for some systems, the changing trends of the calculation results are less obvious than the experimental ones so they have to be magnified as in Fig. 2 (a,q), respectively, for Cr-Co and Re-Cr.

In the present work, the calculations accounting for ferromagnetism were conducted on several sigma compounds of A-Fe (A = Cr, V) to determine the molar volume as shown in Fig. 2 (b,e). The measured magnetic moment for V39.9Fe60.1 compound is 0.63 µB and the Curie temperature is 168.4 K [26]. The calculated magnetic moment by using the corresponding measured site occupancy taken from Ref.

[3] is 0.56 µB. The good agreement between these results further validates our calculation results. For the other systems in Fig. 2, ferromagnetism was not considered. It should be noted that since no sigma compounds known to exist in binary alloy systems exhibit magnetic order above room temperature [26], thus during the whole discussion in the present work, no magnetism is considered.

Definition of order degree

For an A-B binary ordered phase, the order degree is quantitatively expressed by longrange order parameter S [27]:

where PA is the probability of occurrence of atom A on the right crystal site for atom A (i.e. the crystal site is occupied by atom A when the crystal is in totally ordered state); XA is the mole fraction of atom A. When in ordered state, PA = 1 and thus S = 1; When in disordered state, PA =XA and thus S=0. However, Eq. ( 1) is unable to describe the situation when all the right crystal sites for atom A are occupied by atom A (i.e. PA = 1) but there are extra atom A occupying the wrong site (i.e. the crystal site should be occupied by atom B when the crystal is in totally ordered state). Hence we introduce the following description for long-range order parameter S:

where PA and PB are the probabilities of occurrence of atoms A and B (i.e. the site occupancies), respectively on their right crystal sites; XA and XB are the mole fractions of atom A and B, respectively; n, nA and nB are the total number of atoms, the numbers of atom A and B, in a unit cell respectively when the crystal is in totally ordered state. When S = 1, the compound is totally ordered; when the probabilities PA and PB equal the mole fractions XA and XB respectively, the compound is totally disordered with S = 0.

In the present work, we treat the completely ordered compound A2B (B2aA4fA8i1B8i2A8j) defined in Table 1 as a reference ordered state of the sigma phase and extend Eq. ( 2) to the sigma phase as follows:

where atom A has a larger size than atom B; PA and PB are the prob-abilities of occurrence (i.e. the site occupancies) on 2a, 4f, 8i1, 8i2 and 8j sites of atoms A and B, respectively; XA and XB are the mole fractions of atom A and B, respectively. The five terms on the right hand side represent the order degrees of the five inequivalent sites, respectively, designated as S2a, S4f, S8i1, S8i2 and S8j. The long-range order parameter S is hereafter considered as the order degree.

When S=1, the compound is totally ordered; when the prob-abilities PA (i.e. PA2a, PA4f, PA8i1, PA8i2 and PA8j) or PB (i.e. PB2a, PB4f, PB8i1, PB8i2 and PB8j) on the five inequivalent sites equal the mole fractions XA or XB, respectively, the compound is totally disordered with S = 0. When the order degree on the five inequivalent sites (i.e. S2a, S4f, S8i1, S8i2 or S8j) is negative with PB2a < XB, PA4f < XA, PA8i1 < XA, PB8i2 < XB, or PA8j < XA, one deals with an anti-order de-gree. It should be noted that for a specific value of S, there are many possible values of site occupancy. However, for a specific sigma compound, the site occupancy on every site changes within a limited range, and thus the values of S still convey an important piece of information, i.e. compounds with larger values of S are in more ordered state. The calculated order degrees of several sigma compounds based on experimental site occupancies are listed in Table A.1.

Influence of atomic mixing on molar volume of the sigma phase

The influence of atomic mixing on molar volume is described by the volume deviation of the hypothetically disordered sigma phase from the linear composition-volume relationship. It can be seen from Fig. 2 that the sigma compounds in their disordered states do not comply with the above mentioned linearity and present positive or negative deviations, which is summarized in Table 3.

For A-Al, SMe, SMs-SMe categories and two systems (i.e. Ta-Pd and Ta-Rh) in LAs-SMe category, the disordered calculation results present a negative deviation from the linear composition-volume relationship. For LAs, LAs-LAe categories and two systems (i.e. Mo-Co and Mo-Mn) in LAs-SMe category, a positive deviation is observed. It is expected that the deviation from the linearity is related to the tendency of electron loss or gain between the two constituent elements. The tentative dis-cussion is as follows.

For SMe and SMs-SMe categories, large atom A tends to lose elec-trons and small atom B tends to gain the corresponding amounts of electrons [9]. The smaller the atomic size, the stronger the bond be-tween the nucleus and the valence electrons [28]. Thus the transferred electrons to atom B are expected to be more bound. This is consistent with the observed negative deviation for the disordered sigma phase in Fig. 2 and Table 3.

Vice versa, for LAs and LAs-LAe categories, large atom A tends to gain electrons and small atom B tends to lose the corresponding elec-trons [9]. Thus the transferred electrons to atom A are expected to be less bound. This is consistent with the observed positive deviation for the disordered sigma phase in Fig. 2 and Table 3.

With regard to LAs-SMe category, a competition of LAs (i.e. the total number of electron shells of atom A is larger than atom B) and SMe (i.e. the number of valence electrons of atom A is smaller than atom B) factors determines the electron transfer and thus the deviation from the linear composition-volume relationship. For Mo-Co and Mo-Mn, LAs factor is more crucial; for Ta-Pd and Ta-Rh, SMe factor is more crucial. It should be mentioning that the electron transfer behavior for A-Al category is reported to be similar as LAs-SMe category [9].

From the above discussion, it is understandable that among all the systems studied by Crivello et al. [4] as mentioned in Section 1, for Hf-Re and Ta-Re (both belonging to the SMe category), and Nb-Re, Ti-Re and Zr-Re (all belonging to the SMs-SMe category), the calculated cell volumes of the 32 ordered configurations present a negative deviation with respect to composition. The calculations involve the full set of ordered configurations, including the ones in anti-ordered states, for which the cell volumes can be very large due to the size mismatch. For the calculated cell volume of each compound, the negative deviation is presented no matter what the atomic order of the compound is. Thus the crucial factor leading to the negative deviation is the atomic mixing rather than the atomic order. Further, we can conclude the negative deviation is due to a significant effect of SMs (i.e. the total number of electron shells of atom A is smaller than atom B) or SMe factor, i.e. the difference of the total number of electron shells or the number of valence electrons of the constituent elements A and B is large. As concerned to W-Re, Re-Os (both belonging to SMe category), Mo-Re (be-longing to SMs-SMe category), no obvious negative deviation is observed since for these systems the effects of SMs and SMe factors are less remarkable, i.e. the differences of the total number of electron shells and the number of electrons of electrons of the constituent ele-ments A and B are small.

Influence of atomic order on molar volume of the sigma phase

Experimental evidence

In addition to the atomic mixing, the atomic order also affects the molar volume of the sigma phase. Fig. 3 (a) and (b) present the cal-culated order degree of the sigma phase based on experimental site occupancy and the corresponding measured molar volume of the sigma phase from the literature with respect to the composition. As the measured site occupancy is sensitive to alloy processing, instrument accuracy and systematic error etc., for a specific system, only the data measured by the same authors are selected to conduct the comparison. Cr-Fe and Mo-Re are not presented here as the experimental results [3,20,23] indicate that the influence of atomic order on molar volume is small, which will be discussed in Section 3.3.2.

When the rate of change for order degree increases, the rate of change for molar volume decreases as presented in Fig. 4. The re-processing of data from Figs. 3 to 4 excludes the influence of atomic mixing on molar volume of the sigma phase. It thus further indicates that for compounds with a specific composition, when the order degree (S) increases, the molar volume decreases. This is consistent with the calculation results in Fig. 2, i.e. for most binary systems the compounds with smaller order degree have a larger molar volume. One exception is Mo-Re, which will be discussed in Section 3.3.2. Notably, Fig. 2 (r) indicates that for Re-V, the calculated molar volume based on experi-mental site occupancy is even larger than the disordered one. It is because the experimental site occupancy leads to an anti-order degree (i.e. S = -0.236) and thus causes a large mismatch between atomic size and lattice site.

First-principles calculations

The influence of atomic order on molar volume varies with different sigma systems. Fig. 5 indicates that among all the sigma systems, the volume difference between disordered and ordered compounds is the largest for Ta-Al. It indicates that for Ta-Al the atomic order affects the molar volume more dramatically.

Table 4 and Fig. 5 present the volume difference (ΔVm = Vdisorder-Vorder) between disordered and completely ordered states of the A-B (V (A) > V(B)) binary sigma phase at 66.7 at.% A calculated by EMTO-CPA method. It is found that when increasing the total number of electron shells of atom A, s(A) as in Fig. 5a (or decreasing s(B) as in Fig. 5b), or decreasing the number of valence electrons of atom A, e(A) as in Fig. 5c (or increasing e(B) as in Fig. 5d), the volume difference increases. Moreover, increasing s(A) (or decreasing s(B)) indicates en-hancing LAs factor or weakening SMs factor; decreasing e(A) (or in-creasing e(B)) indicates enhancing SMe factor or weakening LAe factor; Thus we can deduce that LAs and SMe factors increase the volume difference while SMs and LAe factors decrease the volume difference. It should be noted that for the four systems (Mo-Co, Mo-Mn, Ta-Pd and Ta-Rh) belonging to LAs-SMe categories, LAs factor is the same, i.e. s(A)-s(B) = 1, and it is expected that the systems (i.e. Ta-Pd and Ta-Rh) with larger SMe factor should have larger volume differences (ΔVm = Vdisorder-Vorder) based on the deduction in the last paragraph. However, Fig. 5 indicates that the systems (i.e. Mo-Co and Mo-Mn) with smaller SMe factor have larger volume differences. This is due to the fact that the volume difference is also affected by the size factor, which in the present work is described by the size difference of the constituent elements, i.e. (VA-VB)VA. Fig. 6 (d) indicates that the larger size factors for Mo-Co and Mo-Mn make their volume difference larger than Ta-Pd and Ta-Rh. For the other systems in Fig. 6, either the effect of size factor on the volume difference is consistent with the electronic factors (e.g. Nb-Al vs Ta-Al, Mo-Mn vs Mo-Re) or the size factor is too small to counteract the effect of the electronic factors (i.e. Cr-Co vs Cr-Fe, V-Co vs V-Fe vs V-Ni, Ta-Pd vs Ta-Rh).

It can also be seen from Table 4 and Fig. 5 that for Mo-Re (belonging to SMs-SMe category), the volume difference is negative, i.e. the or-dered volume is larger than the disordered one. However, Mo-Os also belongs to SMs-SMe category, has the same SMs factor (s(A) -s (B)=-1) as Mo-Re, and the volume difference is positive. The dif-ference is the SMe factor (e(A) -e(B) = -1 for Mo-Re and -2 for Mo-Os) and thus increases the volume difference of Mo-Os to be positive. Among all the stable sigma systems, Mo-Re is the only case, for which the volume difference (ΔV=Vdisorder-Vorder) is negative.

In Section 3.3.1, according to the experimental evidence [3,20,23], we conclude that the influence of atomic order on molar volume is relatively small for Cr-Fe (3d54s1-3d64s2) and Mo-Re (4d55s1-5d56s2) belonging to SMe and SMs-SMe categories respectively. This is inter-pretable by considering the small SMe (i.e. e(A) -e(B) = -2 for Cr-Fe, e (A) -e(B) = -1 for Mo-Re) and SMs (i.e. s(A) -s(B) = -1 for Mo-Re) factors for the these systems, which increase and decrease the volume difference respectively.

To support our arguments, we conducted the calculations for hy-pothetic Nb-W (4d45s1 -5d46s2) and Au-W (5d106s1 -5d46s2) sigma compounds. The two systems belong to SMs-SMe and LAe categories respectively. The volume difference for Nb-W is -0.496 × 10-7 m3/ mol (ΔVm = Vdisorder-Vorder = 10.959e-6 m3/mol -11.008e-6 m3/mol) and for Au-W is -0.779 × 10-7 m3/mol (ΔVm = 10.726e-6 m3/mol -10.804e-6 m3/mol). Both volume differences present negative value due to the effect of SMs or LAe factors.

Conclusions

In the present work, the influence of the atomic mixing on molar volume is described by studying the volume deviation of the hypothetic disordered sigma phase from the linear volume-composition relationship; the influence of atomic order on molar volume is described by the molar volume difference (ΔVm=Vdisorder-Vorder) be-tween the A2B sigma compounds in the ordered and disordered states. To facilitate explanation, the long-range order parameter (i.e. order degree) is redefined, which can describe all the situations of site occupation. The main conclusions can be drawn as follows:

1) The A-B (V(A) > V(B)) disordered sigma phase does not follow a linear volumecomposition relationship. It presents a negative or positive deviation depending on the tendency of electron loss or gain of the two constituent elements. For SMe and SMs-SMe cate-gories, larger atom A tends to lose electrons leading to a denser packing with atomic mixing, and thus the compounds present a negative deviation; for LAs and LAs-LAe categories, larger atom A tends to gain electrons leading to a looser packing with atomic mixing, and thus the compounds present a negative deviation. For LAs-SMe category, the atomic mixing on the molar volume depends on the competition of the electronic factors (i.e. LAs and SMe fac-tors).

2) The influence of the atomic order on the molar volume depends on the size difference and the electron configurations of the two con-stituent elements, i.e. large size difference increases the volume difference; SMe and LAs factors increase volume difference (ΔV=Vdisorder-Vorder); LAe and SMs factors decrease volume differ-ence.

3) For most binary sigma compounds, when increasing the order de-gree, the molar volume decreases. The only exception is Mo-Re belonging to the SMs-SMe category, for which the ordered volume is larger than the disordered one due to the SMs factor causing the volume difference (ΔVm = Vdisorder-Vorder) to be negative.

Acknowledgements

This work was supported by National Key R&D Program of China (Grant number: 2017YFB0701502). W. Liu gratefully acknowledges the financial support from the China Scholarship Council. X.-G. Lu ac-knowledges the High Performance Computing Center of Shanghai University for providing the calculation facility.