Enargite (Cu₃AsS₄): a ductile mid-temperature thermoelectric material

Abstract

Despite their intrinsic large band gap, low carrier concentration, and electrical conductivity, sulfide-based thermoelectric materials have been explored extensively because of their abundance and feasibility. This report uncovers the thermoelectric performance of Cu₃AsS₄ by combining density functional theory, the modified Debye–Callaway model, and phonon Boltzmann transport equations. The overall assessment of thermal, mechanical, and dynamical stability is confirmed through the ab initio molecular dynamics simulations, elastic constants, and phonon dispersion computations. The ultra-low lattice thermal conductivity of 0.327 and 1.020 W m⁻¹ K⁻¹ at 900 K obtained through various approaches can be attributed to the scattering of phonons induced by bonding heterogeneity and large lattice anharmonicity. Also, to improve the reliability of electronic transport properties, the carrier relaxation time is calculated by including acoustic, optical, and impurity phonon scattering mechanisms. The favourable band features and electron and phonon characteristics collectively facilitate a larger optimum power factor accompanied by the figure-of-merit of 1.07 to 2.31 at 900 K for p-type Cu₃AsS₄. These results highlight the potential applicability of Cu₃AsS₄ for mid-temperature thermoelectric applications. Also, this work elaborates the relationship between physical and mechanical characteristics of the crystal structure, which intensifies the understanding from materials to devices.

---

1. Introduction

The substantial enlargement and advancement of renewable energy sectors are essential to overcome global warming and the energy crisis. The solar and wind power sectors are evolving day-by-day and contributing their maximum to society.^1,2 Along these lines, the unique thermoelectric energy conversion technique has the potential to address both issues. Briefly, thermoelectric (TE) devices can convert heat into electricity (generator) and vice versa (refrigerator). The performance of TE devices is determined by the material's dimensionless figure of merit, zT = (σS²/κ)T, where S, σ, κ and T are the Seebeck coefficient, electrical conductivity, thermal conductivity and temperature, respectively.³ That is, a good TE material should have high S and σ with low κ. Conceptually, it is difficult to realize high S and σ simultaneously as they are intrinsically coupled through carrier concentration (n). The complex relationship between σ and κ also resists the tailoring process. Besides these complications, several approaches are developed to enhance the zT of the TE materials.

At least zT > 1 or an efficiency above 15% is required to commercialize the TE devices effectively.⁴ The continuous search in this aspect led to the discovery of several well-performing materials like Bi_0.5Sb_1.5Te₃ (1.47 at 350 K), SnSe (1.4 at 300 K), PbTe (1.5 at 773 K), and GeTe-based compounds (2.14 at 670 K).^5–7 But, the high cost of the precursor elements and the need for sophisticated synthesis approaches drastically increase the cost of electricity production using TE devices.⁸ It is also well-known that ductile materials are crucial to prevent a device's mechanical fractures and failures. However, recent reports disclose the shortage of ductile TE materials.^9,10 Hence, over the past two decades, numerous techniques have been developed and employed to enhance the TE properties of earth-abundant sustainable materials, particularly, (i) to tailor power factor (PF = S²σ) – resonant level doping and band, strain and pressure engineering and (ii) to reduce thermal conductivity (κ_Total = κ_e + κ_L, where κ_e – electronic thermal conductivity and κ_L – lattice thermal conductivity) – nano-structuring, introducing complex crystals, and defect and grain boundary engineering.^11–18 Research has been extended to various materials, including semiconductors, polymers, and nanostructured materials, with considerable improvements in nanostructuring that helps in reducing thermal conductivity via phonon scattering. Additional interest in low-dimensional materials, like 2D materials and nanowires, which exhibit enhanced thermoelectric properties, has been rising. Further advancements emphasize the significance of doping, interface engineering, and the design of heterostructures to further boost the thermoelectric performance.^19–21 In addition to these strategies, it is essential to explore new functional materials for effective implementation of TE technology. The role of high-throughput computations and first principles calculations are unavoidable in this direction.^22,23

Recently, the pnictogen containing I₃–V–VI₄ (I = Cu; V = As, Sb; VI = S, Se) materials have gained widespread attention in various areas, including photovoltaic, photocatalysis and energy storage applications.^24–26 Some of these candidates exhibit noticeable thermoelectric performance ascribed to their attractive carrier concentration and low thermal conductivity. Particularly, Sn, Bi doped Cu₃SbS₃ synthesized by mechanical alloying exhibits a large zT of 0.76 at 623 K. This behaviour is associated with the low κ_L (∼0.75 W m⁻¹ K⁻¹ at 623 K) ascribed to inherent mass fluctuation. Reports also show the possible enhancement of PF by 150% through hole activation via doping.²⁷ Similarly, the Cu₃SbSe₄–CuAlSe₂ composite yields a peak zT of 1.4 at 723 K with a low κ of ∼0.4 W m⁻¹ K⁻¹. This large zT is ascribed to the increase in carrier concentration boosted by creating cation vacancies.²⁸

Also, one of their family members, Cu₃AsS₄, is investigated for photovoltaic, optical, and other semiconducting applications.^29,30 Nevertheless, the reports on the electron, phonon, and thermoelectric properties of Cu₃AsS₄ are uncommon. The enargite (chemically, Cu₃AsS₄) mineral is rarely found in the Longfellow mine, Red Mountain district, San Juan County, Colorado. The enargite crystal exists in a dimorphic state (orthorhombic and tetragonal). The second phase is very rare, while the orthorhombic phase is commonly available in several mines along with the copper ores.³¹ The orthorhombic phase can also be experimentally synthesized by various physical and chemical methods.^32,33

Favourably, the enargite mineral exhibits semiconducting behaviour with holes as the dominating carriers. Allen Pratt analyzed the oxidation states of Cu (1+), As (5+) and S (2−) through X-ray photoemission spectroscopy (XPS).³⁴ Pauporté et al. reported the electrical and optical properties of the natural mineral (NM) and synthetic sample (SS) of Cu₃AsS₄. The highlights are: (i) the carrier concentrations of NM and SS are 1 × 10¹⁷ and 7.8 × 10¹⁹ cm⁻³ at 295 K; (ii) electrical resistances are 6–8 Ω cm⁻¹ (for NM) and 10.5 Ω cm⁻¹ (for SS); and (iii) display semiconducting nature with a band gap (E_g) of 1.43 eV.³⁵ The inherently large n and σ of Cu₃AsS₄ mineral are expected to result in interesting thermoelectric properties. However, to the best of our knowledge, the thermoelectric properties of Cu₃AsS₄ remain unexplored. Hence, the present work is motivated and intended to systematically investigate the transport and thermoelectric properties of Cu₃AsS₄.

To ensure a comprehensive understanding of the thermoelectric properties, multiple state-of-the-art computational approaches are employed to evaluate phonon transport properties, specifically, the modified Debye–Callaway (mDC) model, Cahill–Pohl (CP) model, Slack model and phonon Boltzmann transport equation under single-mode relaxation time approximation (PBTE-RTA) as well as iterative solution (PBTE-Iterative) approaches. This enables the cross-validation of results obtained from different underlying methodologies. Consequently, four different values of zT are presented, each corresponding to a combination of electron and specific phonon transport. Notably, for all combinations, the zT exceeds 1 at 900 K, indicating the potential of Cu₃AsS₄ for mid-temperature thermoelectric applications. The use of multiple approaches not only facilitates cross-comparison but also enhances the reliability of the predicted thermoelectric performance.

2. Computational methods

2.1. Electronic structure calculation

The electronic band structure and related characteristics are calculated using the Vienna ab initio simulation package (VASP) under the plane augmented wave (PAW) method.^36–38 To treat the core and valence electrons, the Perdew–Burke–Ernzerhof (PBE) exchange correlational functional based on generalised gradient approximation (GGA) is used.³⁹ The kinetic energy cutoff of 550 eV along with the converged 7 × 7 × 7 k-mesh is utilized for the entire calculations. The crystal structure is relaxed under the energy and force convergence condition of 1 × 10⁻⁸ eV and 1 × 10⁻⁷ eV Å⁻¹, respectively. Even though the lattice parameters match well with the experimental values, the band gap is strongly underestimated by the PBE functional. In order to resolve this issue, the Heyd–Scuseria–Ernzerhof (HSE06) hybrid functional is adopted for electronic structure calculation.⁴⁰ Moreover, the Effmass code is used to compute the effective mass of the carriers and the Elastic code to calculate the elastic constant of the present system.^41,42

2.2. Relaxation time and electron transport calculation

The electron transport coefficients are obtained by solving the Boltzmann transport equation (BTE), which is employed in the highly efficient AMSET code.⁴³ It solves the BTE based on momentum relaxation time approximation to calculate scattering rates. This code computes the relaxation time (τ) by considering various band and k-point-dependent scattering approaches. In this work, the τ_tot is calculated by accounting: (i) acoustic deformation potential (τ_aco); (ii) polar optical phonon (τ_pop); and (iii) impurity (τ_imp) scattering mechanisms. To obtain the aforesaid scattering rates, the band structure, deformation potential, polar phonon frequency dielectric, and elastic constants are computed using VASP. Furthermore, they are averaged by Matthiessen's rule i.e., 1/τ_tot = 1/τ_aco + 1/τ_pop + 1/τ_imp to obtain τ_tot. Based on the performed convergence test, the interpolation factor of 21 is used to control the interpolation density (refer to Fig. S1, ESI†).

2.3. Phonon properties and lattice thermal conductivity

The phonon dispersion of Cu₃AsS₄ is obtained using a supercell approach under the finite displacement method (FDM) using the phonopy program.⁴⁴ To calculate lattice dynamical properties, the 2 × 2 × 2 supercell with 128 atoms is used. Using the finite displacement method, the supercell is displaced (according to its symmetry). The force of each displaced structure is obtained by performing single point calculation using 3 × 3 × 3 k-point mesh. The dynamical matrix is constructed by collecting all the forces. This dynamical matrix is used to calculate the phonon properties such as phonon dispersion spectra, group velocity, Grüneisen parameters and specific heat.

The κ_L is calculated by various approaches like the modified Debye–Callaway (mDC) model, Cahill–Pohl (CP) model, Slack model and phonon Boltzmann transport equation under single-mode relaxation time approximation (PBTE-RTA) as well as iterative solution (PBTE-Iterative). The mDC, PBTE-RTA and PBTE-Iterative approaches are implemented in AICON2, Phono3py and ShengBTE programs, respectively.^44–46 In the mDC model, the second order IFCs required for calculation are obtained for three different volumes (including equilibrium, slightly smaller and larger volumes). In addition, the group velocity and Grüneisen parameters are also essential. After computing all the necessary parameters, the AICON2 is triggered to calculate κ_L and phonon relaxation time. Usually, the AICON2 calculates κ_L based on the following relation:

κ = κ_TA + κ_TA′ + κ_LA + κ_O (1)

Here, κ_TA, κ_TA′, κ_LA and κ_O represent the thermal conductivity contribution of the transverse, longitudinal acoustic and optical phonon modes. The CP and Slack models are applied manually using elastic constants to estimate minimum and temperature-dependent κ_L, respectively. The detailed description of the relations employed can be found in the ESI,† S1.1. For the PBTE-RTA approach, Phono3py is used, incorporating third-order IFCs and employing the tetrahedron method for the Brillouin zone integration.⁴⁴ The convergence test and direction dependent lattice thermal conductivities are shown in Fig. S2 (ESI†). Here, the converged q-mesh of 15 × 15 × 15 is used. Subsequently, the iterative solution of full linearized phonon Boltzmann transport equation (PBTE-Iterative) is obtained using ShengBTE. Here, the third order IFCs are computed using the thirdorder.py script. The convergence test of κ_L performed as a function of the number of neighbours (nn) is shown in Fig. S3(a) (ESI†). Based on this test, the third order IFCs up to the 22nd nearest neighbour (0.56 nm) are considered, which requires single-point calculations on 3056 displaced structures. A converged q-mesh grid of 13 × 13 × 13 and a scaleboard value of 0.5 are used. All convergence tests are performed with a convergence threshold for a Δκ_L of 0.05 W m⁻¹ K⁻¹. The results of convergence tests and direction dependent lattice thermal conductivities are given in Fig. S3(a–d) (ESI†). The non-analytical term corrections are included both in PBTE-RTA and PBTE-iterative approaches.

3. Results and discussion

3.1. Crystal structure and bonding nature

The enargite mineral in orthorhombic symmetry with the Pmn2₁ space group (no. 31) is studied.⁴⁷ The unit cell fully relaxed using the PBE exchange–correlational functional is displayed in Fig. 1(a), generated with the help of the VESTA program.⁴⁸ The relaxed unit cell parameters a = 6.485 Å, b = 7.470 Å, and c = 6.191 Å match well with earlier results (refer to Table 1). The unit cell consists of six copper, two arsenic, and eight sulfur atoms. Both the cations (Cu and As) are individually surrounded by the four anions (S) in tetrahedral coordination; they are represented in Fig. 1(b) and (c).

Fig. 1 (a) Crystal structure along the a-axis, tetrahedral coordination of (b) Cu–S and (c) As–S and (d) 2D electron localization function for the 2 × 2 × 2 supercell projected along the (0 1 0) plane.

Table 1 The calculated and experimental lattice parameters of enargite mineral

		Lattice parameter (Å)
		a	b	c
Present work	Calc.	6.485	7.470	6.191
Robert R. Seal, II et al.^51	Expt.	6.432	7.402	6.149
Pauling, L. & Weinbaum, S.^52	Expt.	6.460	7.430	6.180
D. J. Vaughan et al.^53	Expt.	6.470	7.440	6.190

---

In addition, the distance between the cation and the anion is analysed. The bond lengths of Cu–S₁, S₂, S₃, and S₄ are 2.312 Å, 2.888 Å, 2.296 Å, and 2.306 Å, whereas the bond lengths of As–S₁, S₂, S₃, and S₄ are 2.245 Å, 2.293 Å, 2.272 Å and 2.293 Å. The unequal bond lengths between the cation and the anion are observed in both the tetrahedra. In order to identify the bonding characteristics, the 2D electron localization function (ELF) is computed for the supercell and illustrated in Fig. 1(d). The ELF values 0 (blue), 0.5 (green) and 1.0 (red) indicate no, equal and complete localization of electrons. Moreover, the ELF value between 0.7 and 0.8 is observed in As–S, which manifests the stronger bonding (mostly covalent) characteristics. On the other hand, lower electron population is noted between Cu and S, indicating lower bonding strength (probably ionic bonding). These unequal bond lengths and bonding heterogeneity can trigger the lattice anharmonicity, which resists phonon propagation and thus lowers the lattice thermal conductivity.

To further understand the thermodynamic stability of the studied compound, the formation and cohesive energies are calculated using eqn (2) and (3). Formation energy gives the idea about the possibility of forming the required phase experimentally and the cohesive energy describes the stability of the compound after synthesis.

Here, E(Cu₃AsS₄) represents the total energy of the system and E(Cu_solid), E(As_solid), and E(S_solid) represent the total energy of the respective elements in their bulk form (in terms of eV per atom), while E(Cu_atom), E(As_atom), and E(S_atom) correspond to the total energy of constituent elements as isolated atoms.^49,50 The calculated formation and cohesive energies of Cu₃AsS₄ are −0.20 eV per atom and −4.08 eV per atom, respectively. The negative sign in both these parameters indicates that Cu₃AsS₄ can be experimentally synthesized and is thermodynamically stable.

3.2. Thermal stability

To examine the thermal stability of enargite mineral, AIMD simulations are performed using the Nóse–Hoover thermostat for 5 ps. The structures are initially equilibrated for 1 ps at the corresponding temperatures. After neglecting the initial 1 ps, the evaluation of free energy as a function of time (up to 4 ps) is shown in Fig. 2(a). The free energy remains almost the same with negligible fluctuations at 300, 600 and 900 K. The evolved structures after the simulation time of 5 ps are shown for both 300 and 900 K in Fig. 2(b) and (c), which indicates that the structure is stable up to 900 K. At 300 K, the structure retains its original atomic arrangements with small distortion, whereas this distortion is intense at 900 K.

Fig. 2 (a) Evolution of free energy at different simulation temperatures. The final structures at (b) 300 K and (c) 900 K. (d) Averaged total and partial pair distribution functions. The atomic trajectory of atoms at (e) 300 K, (f) 600 K and (g) 900 K, respectively.

Fig. 2(d) depicts the calculated radial pair distribution function (PDF) to examine the crystallinity of Cu₃AsS₄. The peaks are well consistent and they broaden at higher temperature up to 600 K. Even though the observed peaks are present up to 900 K, the peaks above 5 Å radius start disappearing. This behaviour can be attributed to the solid to liquid phase transformation. However, the experimental melting point of Cu₃AsS₄ is 944 K.⁵¹ Hence, it is believed that the melting starts around 900 K. To get further insight into this argument, the partial PDF for each possible pair is calculated and depicted along with the total PDF. The first peak at 2.25 Å arises from the superposition of Cu–S and As–S pairs, whereas the second peak (highlighted in yellow colour) at 3.75 Å can be due to Cu–As, Cu–Cu and S–S pairs. It is noted that the fourth peak (highlighted in purple colour) at 5.3 Å is definite both at 300 K and 600 K, while at higher temperature (900 K) it completely disappears. This peak is contributed by Cu–Cu, As–As and S–S pairs. However, the peak is detectable for all pairs even at 900 K, except for the Cu–Cu pair. This behaviour is also observed near 10 Å (highlighted in green). Hence, it can be expected that the solid to liquid phase transition begins with the diffusion of Cu atoms. These characteristics are quite common in chalcogenides such as AgCrSe₂ and Cu₃SbSe₃.^54,55

In order to verify this, the simulated trajectory of atoms in the a–b plane is presented in Fig. 2(e–g). It is noted from Fig. 2(e) that the movement of all atoms in the lattice is similar at 300 K. Upon further heating to 600 K, As and S atoms vibrate at a similar amplitude to the previous one. In contrast, the Cu atoms cover a relatively larger area. This movement of Cu atoms becomes more vibrant at a higher temperature of 900 K. It can also be associated with the weakly bonded (i.e., ionic bonding) Cu–S tetrahedra. In other words, the Cu atoms are weakly bound to the lattice, resulting in large vibrations upon heating. Hence, the Cu atoms move back and forth more vibrantly. This can be the possible reason for the absence of definite peaks above 5 Å in the total PDF. However, there is no bond breakage in Cu–S and As–S tetrahedra, ensuring the crystallinity. This vibration of Cu atoms at higher temperatures results in enhanced tetrahedral distortion and large anharmonicity, which drives the phonon scattering process frequently.

3.3. Elastic constants and mechanical stability

The study of elastic constants plays a vital role in understanding the mechanical stability of the solids. For Cu₃AsS₄, these constants are calculated by the conventional strain approach. According to orthorhombic symmetry, there exist 9 independent elastic constants whose values are listed in Table 2. The Born's criterion is a well-known condition to ensure the mechanical stability of the crystalline solids. It depends on the crystal symmetry. Eqn (4) represents the Born's stability criteria for orthorhombic structure.⁵⁶

C₁₁ > 0; C₁₁C₂₂ > C₁₂²; C₁₁C₂₂C₃₃ + 2C₁₂C₁₃C₂₃ − C₁₁C₁₂² − C₂₂C₁₃² − C₃₃C₁₂² > 0; C₄₄ > 0; C₅₅ > 0; C₆₆ > 0 (4)

Table 2 The calculated elastic constants and mechanical properties of Cu₃AsS₄. Here, B – bulk modulus, G – shear modulus, Y – Young's modulus, υ – Poisson's ratio and CP – Cauchy's pressure

Elastic constants (GPa)		Mechanical properties
C11	146.02	B (GPa)	81.52
C12	44.35	G (GPa)	27.43
C13	27.55	Y (GPa)	73.99
C22	34.52	B/G	2.97
C23	25.88	υ	0.35
C33	173.81	CP (GPa)	111.50
C44	34.52
C55	36.11
C66	49.52

---

The enargite mineral satisfies all these six conditions, which ensures its mechanical stability. Using Voigt–Reuss–Hill approximation, bulk modulus (B), Young's modulus (Y) and shear modulus (G) are calculated (refer to Table 2).^57–59 The calculated bulk and shear moduli are 81.52 GPa and 27.43 GPa, respectively. The larger B, G and Y values indicate the resistive nature of Cu₃AsS₄ towards deformation, fracture, and stiffness. The value of Pugh's ratio (B/G) > 1.75, υ > 0.26 and positive Cauchy's pressure (CP) suggest the ductile nature of the system.⁶⁰ The Pugh's ratio of Cu₃AsS₄ (2.97) is higher than that of Ag₂S (∼2.8) and AgCuS (2.3).^61,62 The compared sulphur based thermoelectric materials are well-known for their excellent mechanical stability (ductility). But, Ag₂S has the constraint of adjusting the carrier concentration via doping without affecting its ductility, while zT of AgCuS requires significant improvement. In this aspect, Cu₃AsS₄ is regarded as an advantageous candidate, which intrinsically possesses excellent metallic like ductility and can outperform well in the device aspect.

3.4. Phonon dispersion and dynamical stability

The phonon dispersion spectra along with phonon DOS are displayed in Fig. 3(a) and (b). The phonon dispersion is calculated along high symmetry points i.e., Γ (0, 0, 0) − X (½, 0, 0) − S (½, ½, 0) − Y (0, ½, 0) − Γ (0, 0, 0) − Z (0, 0, ½) − U (½, 0, ½) − R (½, ½, ½) − T (0, ½, ½) − Z (0, 0, ½) − X (½, 0, 0) − U (½, 0, ½) − Y (0, ½, 0) − T (0, ½, ½) − S (½, ½, 0) − R (½, ½, ½). The presence of positive (real) frequency in the phonon spectrum reveals the dynamical stability of Cu₃AsS₄. Furthermore, 48 branches are observed in the spectra with 3 acoustic and 45 optical modes. Generally, the acoustic mode of vibrations is the primary contributor to the heat transport in many solids.⁶³ Hence, transverse acoustic (TA, TA′) and longitudinal acoustic (LA) modes are highlighted with red, green, and violet colours, respectively. The observed highest acoustic mode frequency is 1.97 THz. And the branches spread up to a maximum of 10.12 THz. Generally, the A–O phonon mode coupling promotes the phonon scattering process and decreases the phonon relaxation time. The acoustic and low-frequency optical modes are primarily contributed by the Cu and As atoms, whereas high-frequency optical modes by the S atoms, as inferred from the phonon density of states (PHDOS) (refer to Fig. 3b). The low frequency A–O coupling in the Cu₃AsS₄ system is due to the vibrations of Cu and As elements, with nearly equivalent atomic masses. The optical–optical (O–O) frequency gap is noted near 5 THz. The temperature dependent PHDOS at various temperatures are extracted from the AIMD simulations and shown in Fig. 3c–e. Also, the population of phonon modes in this O–O phonon band gap region increases with temperature (refer to Fig. 3c–e, yellow highlighted region). The slope of the optical modes decreases, indicating the flattening of these branches. Furthermore, the phonon branches shift towards the low frequency regime, manifesting the reduction of group velocity and thus decrease of lattice thermal conductivity with temperature is expected.

Fig. 3 (a) Phonon dispersion spectrum and phonon density of states at (b) 0 K, (c) 300 K, (d) 600 K and (e) 900 K for enargite.

3.5. Phonon transport properties

3.5.1. Group velocity, Grüneisen parameter, specific heat and Debye temperature. To understand the phonon conduction mechanism in Cu₃AsS₄, it is necessary to evaluate the phonon related properties such as phonon group velocity (v_g), Grüneisen parameter (γ), Debye temperature (θ_D) and specific heat (C_v) to compute κ_L of the system. Generally, the propagation speed of phonons in the lattice can be understood from the group velocity of phonons. The group velocity is calculated using .⁶⁴ The frequency-dependent v_g is depicted in Fig. 4(a). It shows that the predominant contribution to the phonon propagation arises from the acoustic mode, particularly LA modes. Thus, vibration of Cu atoms determines the thermal transportation in Cu₃AsS₄. For better clarity, mode dependent group velocity is illustrated in Fig. 4(b). From this, the average v_g of individual acoustic modes can be obtained from the mathematical average of maximum v_g at the respective Γ-point. The calculated average v_g for TA, TA′ and LA is 2635, 1437 and 1865 m s⁻¹, respectively. Comprehensively, the average v_g for acoustic modes is 2186 m s⁻¹, derived from .⁶⁵ This value for Cu₃AsS₄ is comparable with those of the benchmark TE materials, as evident from Table 3.

Fig. 4 Group velocity as a function of (a) frequency and (b) q-points. The Grüneisen parameter as a function of (c) frequency, (d) q-points and (e) temperature. (f) Specific heat against temperature.

Table 3 Comparison of calculated average group velocity (v̄_g), Debye temperature ([small theta, Greek, macron]_D), Grüneisen parameter ([small gamma, Greek, macron]) and lattice thermal conductivity (κ_L) (at 300 K) of enargite with promising TE materials

	v̄g (m s^−1)	[small theta, Greek, macron]D (K)	[small gamma, Greek, macron]	κL (W m^−1 K^−1)	Ref.
Cu3AsS4	2186	245	1.50	1.173 (mDC)	Present work
				1.024 (Slack)
				2.986 (PBTE-RTA)
				3.164 (PBTE-Iterative)
Cu3SbS4	2725	—	—	4.100	67
Cu3SbSe4	2188	—	1.22	2.910	27
PbTe	1920	130	1.56	2.300	68
SnTe	2160	152	1.50	2.130	68
BiCuSeO	2107	243	1.50	0.550	69

---

The Grüneisen parameter, which determines the degree of anharmonicity in the lattice, is calculated and shown in Fig. 4(c). Generally, the positive and negative γ implies the degree of anharmonicity and negative thermal expansion. Since the acoustic modes are important from the TE perspective, the mode dependent γ is calculated with respect to q-points using the relation: and displayed in Fig. 4(d).⁶⁴ Furthermore, the average Grüneisen parameter, calculated for TA, TA′ and LA is 1.91, 1.86 and 0.38, respectively. It is further averaged mathematically to obtain material's γ. The value of [small gamma, Greek, macron] for the Cu₃AsS₄ system is 1.39. To understand the effect of temperature on lattice anharmonicity, the temperature dependent γ is calculated by implementing quasi harmonic approximation (QHA) and shown in Fig. 4(e). Interestingly, the γ value rises from 1.50 at 300 K to 1.76 at 900 K (∼18% increase). This rapid increase in γ can be attributed to tetrahedral distortion, as evident from the analysis of AIMD simulation. Also, the obtained values of Cu₃AsS₄ are compared with other minerals and chalcogenides in Table 3. The observed anharmonicity facilitates the phonon scattering process at high temperatures.

Due to the direct relationship between κ_L and specific heat capacity, understanding the C_v of the present system is also important. Fig. 4(f) displays the computed C_v as a function of T. Initially, C_v abruptly increases with temperature, obeying the T³ law. Then it gets saturated at higher temperatures, reaching the Dulong–Petit limit. The classical limit according to the Dulong–Petit law (= 3NR, N – the number of atoms in the unit cell, R – ideal gas constant) for the enargite mineral is 399.11 J mol⁻¹ K⁻¹. The calculated C_v falls within this limit.

It is equally important to calculate the Debye temperature of the materials from the aspect of thermoelectrics. The θ_D values for individual acoustic modes are calculated. The values for TA, TA′ and LA are 126, 69 and 90 K, respectively. The global θ_D of enargite calculated from the averaged v_g (, where N and V are the number of atoms and volume of the unit cell respectively) is 245 K.⁶⁶ This value is relatively higher than those of most of the reported thermoelectric materials (see Table 3).

3.5.2. Lattice thermal conductivity and phonon relaxation times. The κ_L is the independent variable in determining the thermoelectric behaviour of the materials. The κ_L values computed using various computational approaches (discussed in Section 2) are presented in Fig. 5(a). The PBTE-RTA and PBTE-Iterative approaches calculate κ_L along different crystallographic directions as displayed in Fig. S2(b) and Fig. S3(d) (ESI†). The geometric mean of these values is used to determine the total κ_L reported in Fig. 5a. The κ_L value obtained both by mDC and Slack models are comparable, whereas those derived from the PBTE-RTA and PBTE-Iterative methods are significantly higher. This discrepancy can be attributed to the strong treatment of phonon dispersion and scattering rates in PBTE, in contrast to the simplified approximations employed in the mDC and Slack models. For instance, at 900 K, κ_L values obtained for the different models are 0.327, 0.341, 0.959 and 1.020 W m⁻¹ K⁻¹, respectively. And they are larger than the minimum thermal conductivity of the system i.e., 0.226 W m⁻¹ K⁻¹. Regardless of the model employed, the observed κ_L decreases with increasing temperature. For example, in the mDC model, at 300 K, the observed κ_L is 1.173 W m⁻¹ K⁻¹, which decreases by ∼72% (0.327 W m⁻¹ K⁻¹) at 900 K. Similarly, κ_L from Slack, PBTE-RTA and PBTE-Iterative approaches also decreases by approximately 66%. This can be ascribed to the larger phonon scattering at higher temperatures. Moreover, κ_L decrease with temperature follows 1/T relation in all the models. This characteristic indicates that the Umklapp three phonon scattering process is dominant in the enargite mineral. Further insights into the scattering mechanisms can be understood from the calculated temperature dependent phonon relaxation time shown in Fig. 5(b). Notably, Umklapp scattering is relatively stronger than isotope scattering for acoustic branches. On the other hand, isotope scattering is predominant in optical modes.

Fig. 5 (a) Lattice thermal conductivity and (b) phonon relaxation time as a function of temperature for enargite obtained in mDC.

3.6. Electronic structure

Fig. 6 demonstrates the calculated electronic properties of Cu₃AsS₄. The projected electronic band structure for the Cu₃AsS₄ system is depicted in Fig. 6(a). The enargite mineral is a direct band gap semiconductor with the VBM and the CBM located at the Γ-point (0, 0, 0). The calculated band gap (E_g) is 1.37 eV, which is in reasonable agreement with experimental reports (refer to Table 4). In addition, the total and projected DOS (Fig. 6b and c) reveal that the VBM is mostly composed of the combination of Cu (d-orbital) and S (p-orbital) atoms, whereas the CBM predominantly arises from the S (p-orbital) atoms. The nearly degenerate bands are observed near the CBM rather than the VBM, named CBM₁ and CBM₂. The energy difference between the CBM and adjacent degenerate bands (such as CBM₁ and CBM₂) is 0.030 and 0.086 eV. The observed band convergence effect has a positive impact from the aspect of improving the electrical conductivity. Even with a slight increase in temperature, the carriers in CBM₁ and CBM₂ will also take part in conduction and thus lead to a significant increase in σ for n-type Cu₃AsS₄. To look in detail, the Fermi surfaces are drawn using the IFermi package and given in Fig. 6(d) and (e).⁷⁰ Here, the Fermi surface is plotted 0.1 eV below the VBM and 0.1 eV above the CBM. It is noted that the Fermi surface near the VBM of Cu₃AsS₄ has only one closed sphere-shaped hole pocket at the centre of the BZ, whereas the CBM has three electron pockets due to the presence of degenerate bands at Γ–Z and Z–X. It is also inferred that the group velocity of the electrons is higher than the holes. This degenerate band and higher carrier group velocity may result in larger electrical conductivity for n-type Cu₃AsS₄.

Fig. 6 (a) Orbital-resolved electronic band structure, (b) total and (c) projected density of states. Fermi surface at (d) 0.1 eV below the VBM and (e) 0.1 eV above the CBM. Band decomposed charge density near the (f) VBM and (g) CBM.

Table 4 The calculated and experimental band gaps for enargite

		Band gap (eV)
Present work	Calc. (PBE)	0.004
Present work	Calc. (HSE06)	1.37
Shi et al.^30	Calc. (HSE06)	1.34
Suzanne K. Wallace^71	Calc. (HSE06 + SOC)	1.24
Durant S. & Parkinsaon B. A.^72	Expt.	1.36

---

Furthermore, the quantification of the effective mass of the carriers (m*) plays a significant role in understanding the transport properties of the material. Here m* is calculated by fitting the band edges. Compared to the CBM, a slightly larger band dispersion is observed at the VBM, resulting in 8% heavier holes than electrons , i.e., 0.31 m₀. Here m₀ is the rest mass of the electron. Generally, large mass leads to less mobility. According to this observation, electrons possess higher mobility in the Cu₃AsS₄ lattice. Furthermore, the charge density near the VBM and the CBM is demonstrated to understand the conduction pathways of the carriers (refer to Fig. 6f and g). Even though both band extrema show similar carrier density, the VBM has slightly higher density, leading to larger carrier pathways for holes. The heavier holes may yield a large Seebeck coefficient for p-type Cu₃AsS₄.

3.7. Electron transport properties

3.7.1. Relaxation time and electrical conductivity. Usually, the thermoelectric properties are calculated both for hole and electron-doped systems. It is inferred from the previous work that Cu₃AsS₄ is a p-type semiconductor.³⁵ Hence, in this work, the transport and thermoelectric properties are limited to hole doping. However, the properties of the electron-doped system are given in the ESI,† S1.2, concerning the future scope of the mineral. Usually, CRTA underestimates the electron transport of the system. To overcome the CRTA, the relaxation time (τ) of the carriers is calculated by accounting for various scattering mechanisms using the AMSET program associated with VASP. Fig. 7(a) shows the calculated carrier concentration dependent τ for holes. It is noted that the τ_aco and τ_pop are almost constant, while τ_imp decreases with n. At the same time, these parameters are highly reliant on temperature. For example, at a specific n of 1 × 10¹⁹ cm⁻³, the value of τ_aco, τ_pop and τ_imp for holes is 13.9 fs (33.0 fs), 311.0 fs (177.0 fs) and 118.0 fs (37.6 fs) at 300 K (900 K). These relaxation times are averaged by Matthiessen's rule. The τ_tot of holes in the Cu₃AsS₄ lattice is 11.90 fs and 2.98 fs at 300 K and 900 K, respectively. The τ_tot is strongly dependent on τ_aco rather than polar optical and impurity phonon scattering mechanisms. Furthermore, the electron relaxation time for the Cu₃AsS₄ lattice is displayed in Fig. S4(a) (ESI†).

Fig. 7 (a) Carrier relaxation time and (b) electrical conductivity of hole-doped Cu₃AsS₄ against carrier concentration at different temperatures.

The calculated σ at various temperatures is shown in Fig. 7(b). The σ increases with n following the relation, σ = neμ. For example, at 300 K, the observed σ is 5019 S m⁻¹ for a specific carrier concentration of 1 × 10¹⁹ cm⁻³. And it is increased to 260 015 S m⁻¹ at n = 1 × 10²¹ cm⁻³ in Cu₃AsS₄. At the carrier concentration of 1 × 10¹⁹ cm⁻³ and elevated temperature of 900 K, σ is 1397 S m⁻¹ for the hole-doped system. The drastic decrease of τ with temperature leads to relatively higher carrier scattering, which results in the decrease in σ with temperature. The electron-doped system yields larger σ when compared to its counterpart (refer to Fig. S4b, ESI†). This can be ascribed to the larger m* of holes than electrons and to the presence of degenerate bands near the CBM. The lighter electrons exhibit more oscillations and can undergo frequent collisions and scattering at higher temperatures. This characteristic yield decreased σ at 900 K for n-type Cu₃AsS₄. The details about the magnitude of the calculated parameters are recorded in Table S1 (ESI†).

3.7.2. Seebeck coefficient and power factor. Fig. 8(a) displays the obtained S as a function of n at different temperatures for hole-doped Cu₃AsS₄. Usually, the Seebeck coefficient is highly reliant on the band structure. In simple words, materials with a larger band gap exhibit higher S. The S value decreases with n by satisfying Mott's relation for both hole and electron-doped Cu₃AsS₄ systems.⁷³ For instance, at 300 K, the S value of 341 and 58 μV K⁻¹ is observed at 300 K for the hole concentration of 1 × 10¹⁹ and 1 × 10²¹ cm⁻³, respectively. Under similar conditions, the S value of hole doped Cu₃AsS₄ is 14% larger than that of the electron-doped counterpart (see Fig. S5a, ESI†). This can be due to the larger m* of holes than that of electrons. Conceptually, larger m* of the carrier can manifest larger S due to its direct relationship. Furthermore, S gets increased with the rise in temperature, following an inverse trend with σ. At 900 K and 1 × 10¹⁹ cm⁻³, S is 508 μV K⁻¹ for p-type Cu₃AsS₄.

Fig. 8 (a) Seebeck coefficient and (b) power factor of hole-doped Cu₃AsS₄ against carrier concentration at different temperatures.

Since S and σ are indirectly related and they are in the numerator of zT relation, it is important to understand the compensation between S and σ. This can be done by calculating PF (S²σ). The PF is one of the deciding parameters in TE materials (i.e.): higher TE performance requires larger PF. Hence, the PF is calculated with respect to n at different temperatures and displayed in Fig. 8(b). The PF increases with carrier concentration up to some extent for both carrier doping and then decreases. Furthermore, the peak PF noted for p-type is 1.49 mW m⁻¹ K⁻² at 300 K. This value is higher than that of the n-type counterpart (0.82 mW m⁻¹ K⁻²), which clearly shows that the major contribution to PF is acquired from S. In other words, even though σ is higher for the electron-doped system, the PF is higher for the hole doped system, indicating the higher reliance of PF on S. The PF further increases with temperature. At 900 K, the peak PF reaches 1.80 mW m⁻¹ K⁻² for p-type Cu₃AsS₄, whereas 0.36 mW m⁻¹ K⁻² for n-type Cu₃AsS₄ (refer to Fig. S5b, ESI†). For better comparison, the peak PF and corresponding carrier concentration are recorded in Table S2 (ESI†).

3.7.3. Electronic and total thermal conductivity. The κ_e is calculated as a function of n at different T for hole and electron-doped Cu₃AsS₄ and displayed in Fig. 9(a) and Fig. S6 (ESI†). The value of κ_e increases with increment of carriers, and at the same time decreases with temperature, following a similar pattern to σ. This behaviour can be understood from the Wiedemann–Franz law (κ_e = LσT, where L is the Lorenz number), which shows their direct relationship. The observed κ_e of p-type Cu₃AsS₄ is 0.034 W m⁻¹ K⁻¹ at n = 1 × 10¹⁹ cm⁻³ and T = 300 K. This observed value is lower than that of its counterpart (0.051 W m⁻¹ K⁻¹) under the same conditions. The reason behind this observation can be ascribed to the lighter mass of the electrons that can move faster than heavier holes, which are conductors of both electricity and heat. At 900 K, the value of κ_e decreases to 0.017 W m⁻¹ K⁻¹ and 0.008 W m⁻¹ K⁻¹ for hole- and electron-doped Cu₃AsS₄.

Fig. 9 (a) Electronic and (b) total thermal conductivity of hole-doped Cu₃AsS₄ at different temperatures. The light-colored lines indicate the corresponding lattice thermal conductivity at specific temperatures obtained via the PBTE-Iterative approach.

The computed total thermal conductivity (κ_Total) for hole-doped Cu₃AsS₄ is depicted in Fig. 9(b). Since κ_L is obtained using different approaches, κ_Total is also calculated accordingly. For simplicity, only the κ_Total obtained using the PBTE-Iterative approach is presented here, while the results based on remaining approaches are provided in the ESI† (see Fig. S7). Interestingly, κ_Total is predominately contributed by the lattice part at low carrier concentration, whereas the electronic contribution is larger at higher concentrations. For instance, at 300 K and 1 × 10¹⁸ cm⁻³, the value of κ_Total is almost equal to κ_L i.e., 3.16 W m⁻¹ K⁻¹. It gets increased by nearly ∼61% at higher concentrations (1 × 10²¹ cm⁻³, i.e., 1.89 W m⁻¹ K⁻¹). This change is more vibrant at higher temperatures due to the larger κ_e contributed by higher carrier oscillations. Furthermore, the κ_Total for the electron-doped system is given in Fig. S8 (ESI†). The values are also tabulated in Table S3 (ESI†) for better clarity.

3.8. Figure of merit

As observed from the earlier discussion, the lattice thermal conductivity computed by adopting various models has significant difference. Hence, it is important to calculate the zT of Cu₃AsS₄ by incorporating κ_L values from different models, in addition to the electron transport properties calculated from AMSET. Fig. 10(a)–(d) displays the calculated zT values as a function of temperature and carrier concentration for hole-doped Cu₃AsS₄. Regardless of the models, zT increases consistently with temperature. For example, the peak zT of p-Cu₃AsS₄ (using mDC) at 300 K is 0.26, which increases by 9 times (i.e., 2.31) at 900 K. This behavior can be attributed to the larger carrier scattering and contribution from degenerate bands at higher temperatures. Also, the maximum zT attained at 900 K with κ_L obtained via Slack, PBTE-RTA and PBTE-Iterative is 2.25, 1.12 and 1.07, respectively. The zT for electron doped Cu₃AsS₄ is displayed in Fig. S9(a–d) (ESI†). The peak zT values and corresponding carrier concentrations are summarized in Table S4 (ESI†) for better clarity. The larger zT in Cu₃AsS₄ can be ascribed to: (i) ultralow lattice thermal conductivity: originating from mixed bonding nature (phonon scattering), A–O coupling (reduced phonon relaxation time), low group velocity (limited phonon propagation), larger Grüneisen parameter (enhanced lattice anharmonicity) and (ii) significant power factor: carrier effective masses and nearly degenerate bands (moderate S and σ).

Fig. 10 The calculated zT using κ_L computed using (a) mDC, (b) Slack, (c) PBTE-RTA and (d) PBTE-Iterative for hole-doped Cu₃AsS₄ and (e) zT and κ_L are compared with the state-of-the-art materials and As-based materials. Alphabets indication: (a) SnSe,⁷⁴ (b) PbTe,⁷⁵ (c) Bi₂Te_3−xSe_x,⁷⁶ (d) Sb doped GeTe,⁷⁷ (e) MgAgSb,⁷⁸ (f) CeFe₄As₁₂,⁷⁹ (g) β-As₂Te₃,⁸⁰ (h) Re₃GeAs₆,⁸¹ (i) CoAsSb,⁸² (j) Eu₅Sn₂As₆⁸³ and (k)–(n) Cu₃AsS₄. The unit of κ_L is W m⁻¹ K⁻¹.

The obtained zT and κ_L are compared with other arsenic-based and benchmark candidates in Fig. 10(e). It is noted from the state-of-the-art that the enargite mineral has the tendency to compete with the Te/Se based benchmark materials. Also, the thermoelectric performance of Cu₃AsS₄ is higher than that of CuFeS₂, Cu₁₂Sb₄S₁₃, and Cu₂SnS₃ – sulfur-based compounds.^84–87 In the case of its applicability towards real-time applications, the use of earth-abundant sulphur and arsenic (despite their toxicity) rather than critical elements (like Se/Te) helps in maintaining the low cost of the device. Although As is toxic, the toxicity of the compound should be prioritized at two key stages: (i) during handling at the synthesis stage, (ii) integration into the devices, and (iii) during disposal into the environment. The first and second stages can be managed by carefully reviewing the safety data sheet, adhering to safety protocols, and using appropriate personnel protective equipment. The third stage can be addressed by following proper disposal procedures. Also, Cu₃AsS₄ exhibits zT greater than 1 at 900 K, which is higher than most of the As-based compounds (blue regime in Fig. 10e). This suggests that the enargite mineral can be an appropriate candidate for mid-temperature TE applications.

Besides the significant thermoelectric performance of Cu₃AsS₄, adopting a suitable strategy to reduce its lattice thermal conductivity further improves its zT. Generally, vacancy and impurity engineering hold significant promise for tuning the thermoelectric properties of Cu-related materials. Recent studies have increasingly focused on the role of intrinsic and extrinsic defects in enhancing TE performance. For example, in Cu₃SbS₃, an excess of sulfur content has been shown to increase both hole carrier concentration and mobility.⁸⁸ Additionally, the elevated sulfur levels enhance phonon scattering, leading to a substantial reduction in the κ_L. As a result, the thermoelectric zT improves by nearly fourfold (i.e., increased from 0.2 to 0.71) compared to the pristine material. Similarly, in Cu₂SnSe₄, the introduction of cationic vacancies has been demonstrated to significantly enhance thermoelectric performance, doubling the zT by promoting strong phonon scattering.⁸⁹ These cationic vacancies effectively reduce κ_L to reach the amorphous limit at 900 K. For instance, κ_L is reduced from 1.1 W m⁻¹ K⁻¹ to 0.6 W m⁻¹ K⁻¹. Given these promising results in related materials, similar defect-engineering strategies could be highly effective in optimizing the TE properties of Cu₃AsS₄. These results highlight an important research avenue and invite further experimental and theoretical research on the use of defects to enhance the thermoelectric efficiency of Cu₃AsS₄.

4. Conclusions

This study presents a comprehensive computational investigation of the structural, mechanical, thermal, and thermoelectric properties of the enargite mineral. The experimentally existing and unexplored (from the TE perspective) orthorhombic Cu₃AsS₄ system is chosen. The unit cell comprises two tetrahedra (Cu–S₄) and (As–S₄), with inequivalent bonding. It led to mixed bonding characters between cations and anions. Mechanical and dynamical stabilities of the present system are confirmed through elastic constants and phonon spectra, respectively. The electronic structure reveals the direct band gap semiconducting nature of enargite. The A–O, O–O coupling at high temperatures, intensive Cu–S tetrahedra distortion, and significant Grüneisen parameter ensure the occurrence of enhanced phonon scattering process and lattice anharmonicity. This favourable structural and phonon properties ended up with an ultralow κ_L from 0.327 to 1.020 W m⁻¹ K⁻¹ at 900 K. The moderate S and σ yield a noticeable PF. The combined effect of ultralow κ_L and moderate PF yields the large zT (above 1) at 900 K for hole-doped Cu₃AsS₄. These findings suggest that enargite can be a suitable system for mid-temperature thermoelectric applications, when compared to other As-based thermoelectric materials. Also, the present work illustrates a clear physical and mechanical picture of Cu₃AsS₄, which helps to achieve better understanding from the device aspect.

Author contributions

Prakash Govindaraj: conceptualization, software, formal analysis, investigation, writing – original draft, visualization, data curation; Hern Kim: supervision, resources, project administration, data curation, funding acquisition, conceptualization; and Kathirvel Venugopal: writing – reviewing and editing, supervision, resources, validation, methodology.

Data availability

The computational results that predominantly support our claims were calculated using the Vienna ab initio simulation package (VASP) that can be found at https://www.vasp.at/ with 10.1103/PhysRevB.54.11169 and 10.1103/PhysRevB.59.1758. The version of the code employed for this study is version 5.4.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the National Research Foundation (NRF) grants funded by the Ministry of Education (RS-2020-NR049576), Republic of Korea. Also, the authors thank High Performance Computing Center, SRM Institute of Science and Technology for providing the computational facility.

References

1. W. Akram, X. Li, S. Ahmed, Z. Ouyang and G. Li, RSC Sustainability, 2025, 3, 21–36.
2. Y. Cao, Y. Meng, Z. Zhang, Q. Yang, Y. Li, C. Liu and S. Ba, Renewable Sustainable Energy Rev., 2024, 196, 114351.
3. X.-L. Shi, J. Zou and Z.-G. Chen, Chem. Rev., 2020, 120, 7399–7515.
4. N. V. Morozova, I. V. Korobeinikov and S. V. Ovsyannikov, J. Appl. Phys., 2019, 125, 220901.
5. P. Xu, K. Jin, Y. Yu, M. Huang, Z. Yan, X. Li, X. Gao, L. Fu and B. Xu, Adv. Funct. Mater., 2024, 2401240.
6. B. Qin, D. Wang, T. Hong, Y. Wang, D. Liu, Z. Wang, X. Gao, Z.-H. Ge and L.-D. Zhao, Nat. Commun., 2023, 14, 1366.
7. S. Wang, C. Chang, S. Bai, B. Qin, Y. Zhu, S. Zhan, J. Zheng, S. Tang and L.-D. Zhao, Chem. Mater., 2023, 35, 755–763.
8. Z. Zhang, M. Gurtaran and H. Dong, ACS Appl. Energy Mater., 2024, 7, 5629–5646.
9. Q. Yang, S. Yang, P. Qiu, L. Peng, T.-R. Wei, Z. Zhang, X. Shi and L. Chen, Science, 2022, 377, 854–858.
10. Z. Gao, Q. Yang, P. Qiu, T. R. Wei, S. Yang, J. Xiao, L. Chen and X. Shi, Adv. Energy Mater., 2021, 11, 2100883.
11. X.-Y. Xin, J. Ma, H.-Q. Liu, Y.-J. Gu, Y.-F. Wang and H.-Z. Cui, J. Mater. Sci. Technol., 2023, 151, 66–72.
12. L. Yu, S.-t Wei, L.-j Wang, Z.-p Zhang, Z. Ji, S.-t Luo, J.-x Liang, W.-y Song and S.-Q. Zheng, ACS Appl. Mater. Interfaces, 2023, 15, 53594–53603.
13. R. Gupta, S. Kakkar, B. Dongre, J. Carrete and C. Bera, ACS Appl. Energy Mater., 2023, 6, 3944–3952.
14. J. Yue, S. Guo, J. Li, J. Zhao, C. Shen, H. Zhang, Y. Liu and T. Cui, Mater. Today Phys., 2023, 39, 101283.
15. N. Tabachkova, M. Shtern, A. Sherchenkov, Y. Shtern, M. Rogachev, V. Panchenko, A. Babich, M. Voronov, M. Tapero and E. Korchagin, Solid State Sci., 2024, 107609.
16. T. Ishibe, J. Chikada, T. Terada, Y. Komatsubara, R. Kitaura, S. Yachi, Y. Yamashita, T. Sato, T. Suemasu and Y. Nakamura, Appl. Phys. Lett., 2021, 119, 141603.
17. Y. Zheng, T. J. Slade, L. Hu, X. Y. Tan, Y. Luo, Z.-Z. Luo, J. Xu, Q. Yan and M. G. Kanatzidis, Chem. Soc. Rev., 2021, 50, 9022–9054.
18. S. Li, Z. Huang, R. Wang, C. Wang, W. Zhao, N. Yang, F. Liu, J. Luo, Y. Xiao and F. Pan, J. Mater. Chem. A, 2021, 9, 11442–11449.
19. C. Artini, G. Pennelli, P. Graziosi, Z. Li, N. Neophytou, C. Melis, L. Colombo, E. Isotta, K. Lohani and P. Scardi, Nanotechnology, 2023, 34, 292001.
20. G. J. Snyder and E. S. Toberer, Nat. Mater., 2008, 7, 105–114.
21. J. Mao, Z. Liu, J. Zhou, H. Zhu, Q. Zhang, G. Chen and Z. Ren, Adv. Phys., 2018, 67, 69–147.
22. B. Andriyevsky, I. Barchiy, I. Studenyak, A. Kashuba and M. Piasecki, Sci. Rep., 2021, 11, 19065.
23. X. Jia, Y. Deng, X. Bao, H. Yao, S. Li, Z. Li, C. Chen, X. Wang, J. Mao and F. Cao, npj Comput. Mater., 2022, 8, 34.
24. U. Chalapathi, B. Poornaprakash and S.-H. Park, Ceram. Int., 2017, 43, 5229–5235.
25. G. García, J. J. Fernández, P. Palacios and P. Wahnón, J. Electron. Mater., 2019, 48, 6777–6785.
26. C. Behera, R. Samal, C. S. Rout, R. S. Dhaka, G. Sahoo and S. L. Samal, Inorg. Chem., 2019, 58, 15291–15302.
27. M. Shen, S. Lu, Z. Zhang, H. Liu, W. Shen, C. Fang, Q. Wang, L. Chen, Y. Zhang and X. Jia, ACS Appl. Mater. Interfaces, 2020, 12, 8271–8279.
28. Y. Huang, B. Zhang, J. Li, Z. Zhou, S. Zheng, N. Li, G. Wang, D. Zhang, D. Zhang and G. Han, Adv. Mater., 2022, 34, 2109952.
29. J. Andler, N. Mathur, F. Zhao and C. Handwerker, Assessing the Potential Environmental Impact of Cu3AsS4 PV Systems, 2019 IEEE 46th Photovoltaic Specialists Conference (PVSC), Chicago, IL, USA, 2019, pp. 1669–1674 DOI:10.1109/PVSC40753.2019.8981146.
30. T. Shi, W.-J. Yin, M. Al-Jassim and Y. Yan, Appl. Phys. Lett., 2013, 103, 152105.
31. P. W. Pohwat, Rocks Miner., 2012, 87, 342–351.
32. A. Müller and R. Blachnik, Thermochim. Acta, 2002, 387, 153–171.
33. M. Fantauzzi, D. Atzei, B. Elsener, P. Lattanzi and A. Rossi, Surf. Interface Anal., 2006, 38, 922–930.
34. A. Pratt, Surf. Interface Anal., 2004, 36, 654–657.
35. T. Pauporté and D. Lincot, Adv. Mater. Opt. Electron., 1995, 5, 289–298.
36. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
37. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169.
38. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558.
39. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
40. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2003, 118, 8207–8215.
41. L. D. Whalley, J. Open Source Softw., 2018, 3, 797.
42. R. Golesorkhtabar, P. Pavone, J. Spitaler, P. Puschnig and C. Draxl, Comput. Phys. Commun., 2013, 184, 1861–1873.
43. A. M. Ganose, J. Park, A. Faghaninia, R. Woods-Robinson, K. A. Persson and A. Jain, Nat. Commun., 2021, 12, 2222.
44. A. Togo, J. Phys. Soc. Jpn., 2023, 92, 012001.
45. T. Fan and A. R. Oganov, Comput. Phys. Commun., 2021, 266, 108027.
46. W. Li, J. Carrete, N. A. Katcho and N. Mingo, Comput. Phys. Commun., 2014, 185, 1747–1758.
47. J. Henao, G. D. de Delgado, J. Delgado, F. Castrillo and O. Odreman, Mater. Res. Bull., 1994, 29, 1121–1127.
48. K. Momma and F. Izumi, J. Appl. Crystallogr., 2011, 44, 1272–1276.
49. J. Zhang, Y. Wang, C. Hua, S. Yang, Y. Liu, J. Luo, T. Liu, J. Nai and X. Tao, Phys. Rev. B, 2022, 106, 235401.
50. L. Cheng, M. Li, L. Zhang and H. Wang, Phys. Rev. B, 2024, 110, 155202.
51. R. R. Seal, R. A. Robie, B. S. Hemingway and H. T. Evans Jr, J. Chem. Thermodyn., 1996, 28, 405–412.
52. L. Pauling and S. Weinbaum, Z. Kristallogr. - Cryst. Mater., 1934, 88, 48–53.
53. P. Velásquez, D. Leinen, J. Pascual, J. Ramos-Barrado, R. Cordova, H. Gómez and R. Schrebler, J. Electroanal. Chem., 2000, 494, 87–95.
54. C. Wang and Y. Chen, npj Comput. Mater., 2020, 6, 26.
55. C. Wang, Y. Wu, Y. Pei and Y. Chen, npj Comput. Mater., 2020, 6, 155.
56. F. Mouhat and F.-X. Coudert, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 90, 224104.
57. R. Hill, Proc. Phys. Soc., Sect. A, 1952, 65, 349.
58. W. Voigt, BG Teubner, Lehrbuch der Kristallphysik (Textbook of crystal physics), Leipzig und Berlin, 1928.
59. A. Reuss, J. Appl. Math. Mech., 1929, 9, 49–58.
60. M. Woodcox, J. Young and M. Smeu, Phys. Rev. B: Condens. Matter Mater. Phys., 2019, 100, 104105.
61. Y. Deng, M. Wei, Y. Lei, J. Lu, P. Peng, Y. Zhang and Z. Zheng, CrystEngComm, 2025, 27, 1055–1077.
62. M. Zhu, X.-L. Shi, H. Wu, Q. Liu and Z.-G. Chen, Chem. Eng. J., 2023, 473, 145236.
63. W. Kim, J. Mater. Chem. C, 2015, 3, 10336–10348.
64. X.-L. Zhu, C.-H. Hou, P. Zhang, P.-F. Liu, G. Xie and B.-T. Wang, J. Phys. Chem. C, 2019, 124, 1812–1819.
65. S. Mukhopadhyay, L. Lindsay and D. J. Singh, Sci. Rep., 2016, 6, 37076.
66. R. A. Robie and J. L. Edwards, J. Appl. Phys., 1966, 37, 2659–2663.
67. A. Suzumura, M. Watanabe, N. Nagasako and R. Asahi, J. Electron. Mater., 2014, 43, 2356–2361.
68. R. Knura, T. Parashchuk, A. Yoshiasa and K. T. Wojciechowski, Dalton Trans., 2021, 50, 4323–4334.
69. L.-D. Zhao, J. He, D. Berardan, Y. Lin, J.-F. Li, C.-W. Nan and N. Dragoe, Energy Environ. Sci., 2014, 7, 2900–2924.
70. A. M. Ganose, A. Searle, A. Jain and S. M. Griffin, J. Open Source Softw., 2021, 6, 3089.
71. S. K. Wallace, K. L. Svane, W. P. Huhn, T. Zhu, D. B. Mitzi, V. Blum and A. Walsh, Sustainable Energy Fuels, 2017, 1, 1339–1350.
72. B. Durant and B. A. Parkinson, Photovoltaic response of naturally occurring semiconducting sulfide minerals, 2016 IEEE 43rd Photovoltaic Specialists Conference (PVSC), Portland, OR, USA, 2016, pp. 2774–2779 DOI:10.1109/PVSC.2016.7750157.
73. N. Wang, M. Li, H. Xiao, X. Zu and L. Qiao, Phys. Rev. Appl., 2020, 13, 024038.
74. L.-D. Zhao, S.-H. Lo, Y. Zhang, H. Sun, G. Tan, C. Uher, C. Wolverton, V. P. Dravid and M. G. Kanatzidis, Nature, 2014, 508, 373–377.
75. T. Fu, X. Yue, H. Wu, C. Fu, T. Zhu, X. Liu, L. Hu, P. Ying, J. He and X. Zhao, J. Materiomics, 2016, 2, 141–149.
76. G. Zhang, B. Kirk, L. A. Jauregui, H. Yang, X. Xu, Y. P. Chen and Y. Wu, Nano Lett., 2012, 12, 56–60.
77. R. K. Vankayala, T. W. Lan, P. Parajuli, F. Liu, R. Rao, S. H. Yu, T. L. Hung, C. H. Lee, S. I. Yano and C. R. Hsing, Adv. Sci., 2020, 7, 2002494.
78. H. Zhao, J. Sui, Z. Tang, Y. Lan, Q. Jie, D. Kraemer, K. McEnaney, A. Guloy, G. Chen and Z. Ren, Nano Energy, 2014, 7, 97–103.
79. A. Watcharapasorn, R. Feigelson, T. Caillat, A. Borshchevsky, G. Snyder and J.-P. Fleurial, J. Appl. Phys., 2002, 91, 1344–1348.
80. J.-B. Vaney, J. Carreaud, G. Delaizir, C. Morin, J. Monnier, E. Alleno, A. Piarristeguy, A. Pradel, A. P. Goncalves and E. B. Lopes, J. Electron. Mater., 2016, 45, 1786–1791.
81. N. Soheilnia, H. Xu, H. Zhang, T. M. Tritt, I. Swainson and H. Kleinke, Chem. Mater., 2007, 19, 4063–4068.
82. X. Tan, K. P. Devlin, X. Deng, C.-J. Kang, M. Croft, C. E. Frank, C. Pak, S. Lapidus, S. M. Kauzlarich and V. Taufour, Chem. Mater., 2018, 30, 4207–4215.
83. J. Wang, S.-Q. Xia and X.-T. Tao, Inorg. Chem., 2012, 51, 5771–5778.
84. A. V. Powell, J. Appl. Phys., 2019, 126, 100901.
85. T. Bernard, M. A. Malagutti, K. Lohani, M. D’Incau, N. Ataollahi and P. Scardi, J. Mater. Sci., 2024, 59, 15491–15503.
86. B. S. Shah, J. B. Raval, D. Kumar, S. H. Chaki and M. Deshpande, J. Alloys Compd., 2023, 938, 168566.
87. M. A. Malagutti, K. Lohani, I. Caño Prades, A. Navarro-Güell, T. Bernard, A. Chiappini, E. Saucedo, N. Ataollahi and P. Scardi, ACS Appl. Electron. Mater., 2024, 7, 21–32.
88. L. Li, J. Wan, C. Lin, Y. Wu, Z. Cheng, Y. Lin, S. Yang and M. Luo, Small, 2025, 2503137.
89. W. Li, S. Lin, X. Zhang, Z. Chen, X. Xu and Y. Pei, Chem. Mater., 2016, 28, 6227–6232.

---

Footnote

† Electronic supplementary information (ESI) available: The convergence tests to calculate electron transport properties; the detailed description of Cahill–Pohl and Slack models; the convergence test and lattice thermal conductivity calculated through phono3py and ShengBTE; the relaxation time, electrical conductivity, Seebeck coefficient, thermoelectric power factor, electronic and total thermal conductivity, and figure-of-merit for the electron doped system; and the transport and thermoelectric properties at specific carrier concentration and temperature. See DOI: https://doi.org/10.1039/d5tc00351b

---