Investigation of PbSnTeSe High-Entropy Thermoelectric Alloy: A DFT Approach

Thermoelectric materials have attracted extensive attention because they can directly convert waste heat into electric energy. As a brand-new method of alloying, high-entropy alloys (HEAs) have attracted much attention in the fields of materials science and engineering. Recent researches have found that HEAs could be potentially good thermoelectric (TE) materials. In this study, special quasi-random structures (SQS) of PbSnTeSe high-entropy alloys consisting of 64 atoms have been generated. The thermoelectric transport properties of the highest-entropy PbSnTeSe-optimized structure were investigated by combining calculations from first-principles density-functional theory and on-the-fly machine learning with the semiclassical Boltzmann transport theory and Green–Kubo theory. The results demonstrate that PbSnTeSe HEA has a very low lattice thermal conductivity. The electrical conductivity, thermal electronic conductivity and Seebeck coefficient have been evaluated for both n-type and p-type doping. N-type PbSnTeSe exhibits better power factor (PF = S2σ) than p-type PbSnTeSe because of larger electrical conductivity for n-type doping. Despite high electrical thermal conductivities, the calculated ZT are satisfactory. The maximum ZT (about 1.1) is found at 500 K for n-type doping. These results confirm that PbSnTeSe HEA is a promising thermoelectric material.


Introduction
In recent years, the energy problem has become increasingly serious. Thermoelectric materials have attracted extensive attention as they can directly convert waste heat into electric energy [1][2][3][4]. The figure of merit ZT = S 2 σT/κ [5] is used to evaluate the thermoelectric conversion efficiency in which S, σ, T, and κ are the Seebeck coefficient, electrical conductivity, temperature, and thermal conductivity (electronic and lattice), respectively. However, due to the strong coupling between these parameters, it is not easy to improve the thermoelectric efficiency [6]. Therefore, in order to improve the thermoelectric efficiency, various routes have been explored, e.g., using band engineering to improve power factor (S 2 σ) [7,8], and reducing the dimensionality of the material to reduce the lattice thermal conductivity [9][10][11]. In addition to improving the thermoelectric efficiency of existing materials with these strategies, finding new thermoelectric materials is also an important approach [12][13][14].

Computational Details
Our computations were performed within the DFT framework utilizing the projector augmented wave (PAW) [67] technique, as implemented in the Vienna Ab Initio Simulation Package (VASP) [68][69][70]. The Perdew-Burke-Ernzerhof (PBE) functional under the generalized gradient approximation (GGA) [71] was utilized as the exchange-correlation functional. The nonempirical PBE functional is known to yield accurate crystal parameters and properties, and is known to fulfil many sum rules on the exchange-correlation hole [72]. Generally, where GGA functionals fail, local density functionals fail too. Concerning the thermal conductivity, compared to local density functionals, PBE does not over bind structures (it slightly under binds), hence the interatomic force constants are not too soft. Overall, the lattice thermal conductivity is quite the same whether using local density or GGA functionals [73].
The special quasi-random structures (SQS) of high-entropy alloys have been generated using the Monte Carlo SQS (MCSQS) tool as implemented in the alloy theoretic automated toolkit (ATAT) [74]. A 2 × 2 × 2 supercell consisting of 64 atoms was built for calculations. A 3 × 3 × 3 Monkhorst-Pack k-point mesh was used, and the kinetic energy cutoff was set to 400 eV. The geometric structures were totally relaxed until the Hellmann-Feynman forces were less than 0.01 eV/Å. The electronic transport properties were computed with the Boltztrap2 code, which implements the semiclassical Boltzmann theory under the relaxation time approximation (RTA) [75]. As the calculation of transport coefficients is very demanding in term of the band structure accuracy, a much denser, 9 × 9 × 9 k-point mesh was employed to obtain the electrons energy eigenvalues for the subsequent electronic transport properties calculation. As within the RTA the lifetime of electrons has to be determined separately, which constitutes an inherent limitation to the linearized Boltzmann transport equation approach, the relaxation time (τ) was determined using the deformation potential (DP) theory [76]. Due to heavy elements present in the structure, the spin-orbit coupling (SOC) effect was considered in our calculations. The lattice thermal conductivity was calculated using the Green-Kubo theory for which the heat flux was obtained from the on-the-fly machine-learned force fields [65,66] module of VASP. The evident advantage of using this approach is that classical molecular dynamics allows for catching the phonon dynamics at any order and is applicable to potentially extremely large structures. The difficulty is that interatomic parameters of good quality have to be available. The first step of the MLFF elaboration consists in training the force field by machine learning through molecular dynamics (MD) simulations in the NVT ensemble, where N, V, and T are the number of particles, volume and temperature, respectively, that are kept constant during the simulations (canonical ensemble). A supercell of 512 atoms was used for the training, and the MD simulations were run with a time step of 1 fs. Then, after training, the force field was used to equilibrate the system in the NVT ensemble at the desired temperature with a time step of 1 fs for 100 ps. Finally, the heat flux was calculated through molecular dynamics simulations in the NVE ensemble, where N, V and the energy E are kept constant (microcanonical ensemble), with a time step of 1 fs for 100 ps. For the ensemble average, 10 independent molecular dynamics simulations were performed for the calculation of the heat flux.
The design of the force-field (FF) is described in detail in Ref. [59]. We give a brief summary here. The fitting of the force-field parameters relies on the availability of a database (DB) of structures, the quality of the fitting being assessed from energies, atomic forces, and stress tensors. This DB is built on-the-fly through the ab initio molecular dynamics (AIMD) simulations. At each step of the AIMD a decision is made as to whether a new AIMD step should be run to add a new structure to the DB or if the MD step is run with the FF. The decision is made after the estimated errors between the ab initio atomic forces and the FF ones based on a Bayesian inference. Hence, the algorithm relies on the Bayesian linear regression to assess the quality of the FF parameters. In our case, the force-field parameters were built from a DB containing more than 2100 structures. The Bayesian error on the atomic forces, energies and stress are below 0.006 eV/Å, 0.5 meV and 0.04 kB, respectively. The FF parameters quality is assessed based on the capability of the FF to reproduce properly the two-body and three-body distributions (Equations (2) and (3) in Ref. [59]) that represent the likelihood to find, around a given atom, an atom at a certain distance or a pair of atoms at a certain distance and angle, respectively. Thereupon, the calculated cell parameters with the FF yields are the same as those obtained ab initio.

Structural and Electronic Properties of PbSnTeSe High-Entropy Alloy
The PbSnTeSe high-entropy alloy is based on the NaCl-type face-centered cubic (FCC) crystal structure and a 2 × 2 × 2 supercell containing 64 atoms was constructed. Figure 1a shows the schematic illustration of the crystal structure of the PbSnTeSe HEA. The HEA structure was built from the PbTe crystal structure for which the same atomic ratio of 50% was employed for both Pb and Sn on the Pb Wyckoff positions, and Te and Se on the Te Wyckoff positions. Figure 1b  Bayesian linear regression to assess the quality of the FF parameters. In our case, the forcefield parameters were built from a DB containing more than 2100 structures. The Bayesian error on the atomic forces, energies and stress are below 0.006 eV/Å , 0.5 meV and 0.04 kB, respectively. The FF parameters quality is assessed based on the capability of the FF to reproduce properly the two-body and three-body distributions (Equations (2) and (3) in Ref. [59]) that represent the likelihood to find, around a given atom, an atom at a certain distance or a pair of atoms at a certain distance and angle, respectively. Thereupon, the calculated cell parameters with the FF yields are the same as those obtained ab initio.

Structural and Electronic Properties of PbSnTeSe High-Entropy Alloy
The PbSnTeSe high-entropy alloy is based on the NaCl-type face-centered cubic (FCC) crystal structure and a 2 × 2 × 2 supercell containing 64 atoms was constructed. Figure 1a shows the schematic illustration of the crystal structure of the PbSnTeSe HEA. The HEA structure was built from the PbTe crystal structure for which the same atomic ratio of 50% was employed for both Pb and Sn on the Pb Wyckoff positions, and Te and Se on the Te Wyckoff positions. Figure 1b Figure 2. PbSnTeSe is a semiconductor with a direct band gap. Without spin-orbit coupling, the valence-band maximum (VBM) and the conduction-band minimum (CBM) are both located at the  point. Because of the heavy elements Pb, Te and Sn, the spin-orbit coupling effect was accounted for. The band gap decreases from 0.21 eV to 0.08 eV when considering the SOC effects, and the direct bandgap shifts slightly away from . In addition, the SOC lifts the degeneracy of the crystal states leading to two states (j = l ± ½ ) between the high symmetry k-points, which can be explained by the main contribution of the p atomic orbitals of the atoms.  Figure 2. PbSnTeSe is a semiconductor with a direct band gap. Without spin-orbit coupling, the valence-band maximum (VBM) and the conduction-band minimum (CBM) are both located at the Γ point. Because of the heavy elements Pb, Te and Sn, the spin-orbit coupling effect was accounted for. The band gap decreases from 0.21 eV to 0.08 eV when considering the SOC effects, and the direct bandgap shifts slightly away from Γ. In addition, the SOC lifts the degeneracy of the crystal states leading to two states (j = l ± 1 2 ) between the high symmetry k-points, which can be explained by the main contribution of the p atomic orbitals of the atoms.

Seebeck Coefficient, Electrical Conductivity, and Power Factor of PbSnTeSe
Based on the calculated electronic structure, the Seebeck coefficient (S) of PbSnTeSe was determined. Figure 3 shows the Seebeck coefficient at 300 K, 500 K, and 700 K as a

Seebeck Coefficient, Electrical Conductivity, and Power Factor of PbSnTeSe
Based on the calculated electronic structure, the Seebeck coefficient (S) of PbSnTeSe was determined. Figure 3 shows the Seebeck coefficient at 300 K, 500 K, and 700 K as a function of carrier concentration. For both n-type and p-type doping, the Seebeck coefficients first rise then decrease as the carrier concentration increases, which can be interpreted from the Mott formula [77] where h, k B , m* d , T and n are the Planck constant, Boltzmann constant, density of states effective mass and carrier concentration, respectively. According to this expression, apart from the effect of temperature, the Seebeck coefficient is governed by the ratio m* d n −2/3 . Assuming the simple evolution of the density of states (DOS) for 3D materials as the square-root of the state energies (see, e.g., Figure 39.1 in Ref. [78]), at low doping level the curvature radius increases drastically and hence the DOS mass, and overall, the Seebeck coefficient increases sharply. As the doping level increases further, the DOS mass becomes roughly constant, and the n −2/3 term starts dominating in the Mott formula, leading to a decrease of S. For both low and highly doped compound (n~10 17 cm −3 and n ≥ 10 21 cm −3 ), the Seebeck coefficient is improved with the increase in the temperature, which can also be understood by this formula. In the meantime, the maximum S values are reduced with the temperature increase. The peak values of S for PbSnTeSe are all in the range of 160-190 µV/K with little difference between n-type and p-type.

Seebeck Coefficient, Electrical Conductivity, and Power Factor of PbSnTeSe
Based on the calculated electronic structure, the Seebeck coefficient (S) of PbSnTeSe was determined. Figure 3 shows the Seebeck coefficient at 300 K, 500 K, and 700 K as a function of carrier concentration. For both n-type and p-type doping, the Seebeck coefficients first rise then decrease as the carrier concentration increases, which can be interpreted from the Mott formula [77] where h, kB, m*d, T and n are the Planck constant, Boltzmann constant, density of states effective mass and carrier concentration, respectively. According to this expression, apart from the effect of temperature, the Seebeck coefficient is governed by the ratio m*d n −2/3 . Assuming the simple evolution of the density of states (DOS) for 3D materials as the square-root of the state energies (see, e.g., Figure 39.1 in Ref. [78]), at low doping level the curvature radius increases drastically and hence the DOS mass, and overall, the Seebeck coefficient increases sharply. As the doping level increases further, the DOS mass becomes roughly constant, and the n −2/3 term starts dominating in the Mott formula, leading to a decrease of S. For both low and highly doped compound (n  10 17 cm −3 and n  10 21 cm −3 ), the Seebeck coefficient is improved with the increase in the temperature, which can also be understood by this formula. In the meantime, the maximum S values are reduced with the temperature increase. The peak values of S for PbSnTeSe are all in the range of 160-190 μV/K with little difference between n-type and p-type.  As in RTA the electrical conductivity is scaled by the carrier relaxation time τ, this parameter has to be determined to get the values of σ. For this, the deformation potential (DP) theory [76] was utilized, from which the expression of τ reads: where and T are the reduced Planck constant and temperature, respectively. The effective mass of the carrier is calculated by m * = 2 / ∂ 2 E/∂k 2 , and the elastic constant is defined as C = ∂ 2 E/∂(∆a/a 0 ) 2 /V 0 where E, ∆a and V 0 are the total energy of the system, the change of the lattice parameter and the equilibrium volume, respectively. The DP constant E 1 corresponds to the shift of the band edge energy and is given by where E edge is the band edge energy. The calculated τ values are listed in Table 1. At the temperatures of 300 K, 500 K, and 700 K, the relaxation times for n-type and p-type PbSnTeSe vary from 158 to 44.3 fs and from 17.3 to 4.78 fs, respectively. The relaxation time for the n-type PbSnTeSe is larger than that of the p-type because the effective mass and deformation potential of the n-type compound are smaller than those of the p-type. The decrease in the relaxation time with increasing temperature means that the scattering of carriers is gradually enhanced. The scattering inhibits the transport of carriers, which should lead to a decrease of the conductivity with increasing temperature. Based on the calculated relaxation time, the electrical conductivity (σ) of PbSnTeSe was obtained. Figure 4 shows the calculated σ at 300 K, 500 K and 700 K as a function of the carrier concentration. For both n-type and p-type doping, σ increases with the increase in carrier concentration. For low carrier concentrations σ increases with temperature, and for high concentrations (>10 19 cm −3 ) σ decreases with the increase in temperature, which can be understood from the Drude-Sommerfeld formula [79][80][81]: where n is the carrier concentration and µ is the mobility of the carriers. As mentioned above, with the temperature increase the relaxation time decreases, which leads to a decrease in mobility, hence the conductivity decreases as the temperature increases. In each case, the n-type PbSnTeSe has a better electrical conductivity than the p-type, which is caused by a larger relaxation time for the n-type than for the p-type PbSnTeSe.  Based on the calculated and , the power factor PF ( 2 ) was determined. Fig  5 shows the calculated PF at 300 K, 500 K, and 700 K as a function of carrier concentrat For both n-type and p-type doping, the PF first increases then decreases with the incr in carrier concentration. The peak values of PF for p-type and n-type PbSnTeSe are in range of 0.9-1.2 mW/(m K 2 ) and 7-8 mW/(m K 2 ), respectively. In each case, the optima values for n-type doping are larger than those for p-type doping because of the hig electrical conductivity, indicating that n-type doping is more efficient than the p-typ improving the TE performance of PbSnTeSe. Based on the calculated S and σ, the power factor PF (S 2 σ) was determined. Figure 5 shows the calculated PF at 300 K, 500 K, and 700 K as a function of carrier concentration. For both n-type and p-type doping, the PF first increases then decreases with the increase in carrier concentration. The peak values of PF for p-type and n-type PbSnTeSe are in the range of 0.9-1.2 mW/(m K 2 ) and 7-8 mW/(m K 2 ), respectively. In each case, the optimal PF values for n-type doping are larger than those for p-type doping because of the higher electrical conductivity, indicating that n-type doping is more efficient than the p-type at improving the TE performance of PbSnTeSe. 5 shows the calculated PF at 300 K, 500 K, and 700 K as a function of carrier concentra For both n-type and p-type doping, the PF first increases then decreases with the incr in carrier concentration. The peak values of PF for p-type and n-type PbSnTeSe are in range of 0.9-1.2 mW/(m K 2 ) and 7-8 mW/(m K 2 ), respectively. In each case, the optima values for n-type doping are larger than those for p-type doping because of the hi electrical conductivity, indicating that n-type doping is more efficient than the p-typ improving the TE performance of PbSnTeSe.

Electronic Thermal Conductivity of PbSnTeSe
Thermal conductivity is composed of two parts, the electronic thermal conducti ( ) and the lattice thermal conductivity ( ). The electronic thermal conductivity was culated using the Wiedemann-Franz law [82,83] where L was approximated by the Lorenz number L0 that takes the value 2.44 × 10 8 W Whereas the deviation of the L/L0 ratio from one is still an open question for nanos materials [84], it seems that the deviation from the Wiedemann-Franz law occurs ma at low temperature (well below 300 K) where lattice vibrations increase the L/L0 r above 1. However, this tendency can also be counteracted by electronic corrections, l ing finally to a small change of the L/L0 ratio (between 0.8 and 1.2, at most). For bulk c pounds, the same effects can occur. Therefore, we are confident that the conclus Figure 5. Power factor (PF) for n-type (a) and p-type (b) PbSnTeSe at 300 K, 500 K and 700 K as a function of carrier concentration.

Electronic Thermal Conductivity of PbSnTeSe
Thermal conductivity is composed of two parts, the electronic thermal conductivity (k e ) and the lattice thermal conductivity (k L ). The electronic thermal conductivity was calculated using the Wiedemann-Franz law [82,83], where L was approximated by the Lorenz number L 0 that takes the value 2.44 × 10 8 WΩK −1 .
Whereas the deviation of the L/L 0 ratio from one is still an open question for nanoscale materials [84], it seems that the deviation from the Wiedemann-Franz law occurs mainly at low temperature (well below 300 K) where lattice vibrations increase the L/L 0 ratio above 1. However, this tendency can also be counteracted by electronic corrections, leading finally to a small change of the L/L 0 ratio (between 0.8 and 1.2, at most). For bulk compounds, the same effects can occur. Therefore, we are confident that the conclusions presented hereafter should not change drastically with L. Figure 6 shows the calculated k e at 300 K, 500 K, and 700 K as a function of the carrier concentration. Due to the linear correlation between k e and σ, the impact of carrier concentration and temperature on the electronic thermal conductivity is the same to that on the electrical conductivity. For both n-type and p-type doping, k e increases with the increase in carrier concentration, and increases (decreases, resp.) with temperature for low (high, resp.) carrier concentrations. In each case, the n-type PbSnTeSe has a larger k e than the p-type. presented hereafter should not change drastically with L. Figure 6 shows the calcul at 300 K, 500 K, and 700 K as a function of the carrier concentration. Due to the li correlation between and , the impact of carrier concentration and temperature on electronic thermal conductivity is the same to that on the electrical conductivity. For n-type and p-type doping, increases with the increase in carrier concentration, increases (decreases, resp.) with temperature for low (high, resp.) carrier concentrati In each case, the n-type PbSnTeSe has a larger than the p-type. Figure 6. Electronic part of the thermal conductivity (κe) for n-type (a) and p-type (b) PbSnTe 300 K, 500 K, and 700 K as a function of carrier concentration.

Machine-Learned Force-Field Potential
To build an interatomic potential force field (FF) for PbSnTeSe the on-the-fly chine-learned FF algorithm integrated in the VASP code was used. The training stra of MLFF consists in constructing the force field on the fly during MD simulation and predicted Bayesian error is used at every MD step to judge whether additional first-p ciple calculations need to be performed and a new structure be included in the datas not. When the force field is trained by the on-the-fly machine-learning algorithm, m Figure 6. Electronic part of the thermal conductivity (κ e ) for n-type (a) and p-type (b) PbSnTeSe at 300 K, 500 K, and 700 K as a function of carrier concentration.

Machine-Learned Force-Field Potential
To build an interatomic potential force field (FF) for PbSnTeSe the on-the-fly machinelearned FF algorithm integrated in the VASP code was used. The training strategy of MLFF consists in constructing the force field on the fly during MD simulation and the predicted Bayesian error is used at every MD step to judge whether additional first-principle calculations need to be performed and a new structure be included in the dataset or not. When the force field is trained by the on-the-fly machine-learning algorithm, many of the MD steps are carried out with the force field, and first-principles calculations are executed only when the predicted Bayesian error is large. Figure 7 depicts the estimated Bayesian error of the MLFF, which shows that it is consistently lowering. Table 2 shows the root-means-square errors (RMSE) in the energies, forces, and stress tensors predicted by MLFFs for the training dataset. The predicted errors are low, which agrees with the general observation that the Bayesian linear regression (BLR), which was adopted to obtain the regression coefficients of kernel-based methods, leads to lower errors than other methods [59].

Lattice Thermal Conductivity of PbSnTeSe
According to the Green-Kubo theory, the thermal conductivity and the heat flux are related by [64]: where k B , T and V denote the Boltzmann constant, the temperature, and the volume of the system, respectively. j(t) is the heat flux and the symbols < · > represent the ensemble average over every MD simulation. Based on the heat flux calculated by the MLFF, the heat-flux autocorrelation function (HFACF) of PbSnTeSe was obtained. Figure 8a shows the calculated normalized averaged HFACF at the temperature of 300 K as a function of correlation time. At the beginning, the HFACF starts at one and rapidly drops to oscillate around zero. As correlation time increases, the oscillation gradually decreases, and finally approaches zero, indicating that the HFACF has converged. The HFACF is then used to determine the thermal conductivity (κ L ) of PbSnTeSe. Figure 8b shows the calculated κ L at the temperature of 300 K as a function of correlation time. The trend of κ L over correlation time is the same as that of HFACF. Initially, the oscillation is relatively large and finally tends to stabilize converging towards a constant value, hence asserting the proper convergence of our simulations.
Materials 2023, 16, x FOR PEER REVIEW 10 ensemble average over every MD simulation. Based on the heat flux calculated by MLFF, the heat-flux autocorrelation function (HFACF) of PbSnTeSe was obtained. Fi 8a shows the calculated normalized averaged HFACF at the temperature of 300 K function of correlation time. At the beginning, the HFACF starts at one and rapidly d to oscillate around zero. As correlation time increases, the oscillation gradually decrea and finally approaches zero, indicating that the HFACF has converged. The HFAC then used to determine the thermal conductivity ( ) of PbSnTeSe. Figure 8b shows calculated at the temperature of 300 K as a function of correlation time. The tren over correlation time is the same as that of HFACF. Initially, the oscillation is relati large and finally tends to stabilize converging towards a constant value, hence asser the proper convergence of our simulations. Using the same approach, the value of was computed for various temperat in the range 300-700 K (Figure 9). At 300 K, the value is 0.4 W K −1 m −1 , which is a low value of the lattice thermal conductivity, favorable for thermoelectric materials. is due to the lattice distortion effect of high-entropy alloys that can reduce phonon velo and enhance phonon scattering, resulting in low thermal conductivity. Additionally, can observe that when the temperature rises, the thermal conductivity of the lattice creases, which is also due to an increase in phonon scattering at higher temperatures Using the same approach, the value of κ L was computed for various temperatures in the range 300-700 K (Figure 9). At 300 K, the κ L value is 0.4 W K −1 m −1 , which is a very low value of the lattice thermal conductivity, favorable for thermoelectric materials. This is due to the lattice distortion effect of high-entropy alloys that can reduce phonon velocity and enhance phonon scattering, resulting in low thermal conductivity. Additionally, one can observe that when the temperature rises, the thermal conductivity of the lattice decreases, which is also due to an increase in phonon scattering at higher temperatures.
in the range 300-700 K (Figure 9). At 300 K, the value is 0.4 W K m , which is a very low value of the lattice thermal conductivity, favorable for thermoelectric materials. This is due to the lattice distortion effect of high-entropy alloys that can reduce phonon velocity and enhance phonon scattering, resulting in low thermal conductivity. Additionally, one can observe that when the temperature rises, the thermal conductivity of the lattice decreases, which is also due to an increase in phonon scattering at higher temperatures.

Figure of Merit of PbSnTeSe
Based on the electronic and thermal transport coefficients, the ZT value can be determined. Figure 10 shows the calculated ZT at 300 K, 500 K, and 700 K as a function of the carrier concentration. For both p-type and n-type PbSnTeSe, the ZT value first increases to an optimal value then decreases with increasing concentration. The peak value of ZT for n-type does not vary significantly as the temperature rises, only the ideal carrier concentration shifts a little towards higher values. However, for p-type, the peak of ZT increases significantly with increasing temperature, and the ideal carrier concentration shifts to higher values. The best ZT value for n-type PbSnTeSe is 1.1 at 500 K, while for p-type doping it amounts to 0.75 at 700 K. At each temperature, the ZT value for n-type compound is greater than that of p-type, mainly resulting from the large electrical conductivity for n-type doping, indicating that n-type PbSnTeSe exhibits better thermoelectric properties than the p-type.

Figure of Merit of PbSnTeSe
Based on the electronic and thermal transport coefficients, the ZT value can be determined. Figure 10 shows the calculated ZT at 300 K, 500 K, and 700 K as a function of the carrier concentration. For both p-type and n-type PbSnTeSe, the ZT value first increases to an optimal value then decreases with increasing concentration. The peak value of ZT for n-type does not vary significantly as the temperature rises, only the ideal carrier concentration shifts a little towards higher values. However, for p-type, the peak of ZT increases significantly with increasing temperature, and the ideal carrier concentration shifts to higher values. The best ZT value for n-type PbSnTeSe is 1.1 at 500 K, while for ptype doping it amounts to 0.75 at 700 K. At each temperature, the ZT value for n-type compound is greater than that of p-type, mainly resulting from the large electrical conductivity for n-type doping, indicating that n-type PbSnTeSe exhibits better thermoelectric properties than the p-type.

Comparison with Available Data on PbSnTeSe and Other HEA
As mentioned in the introduction section, Fan et al. [61] and Raphel et al. [62,63] have reported experimental investigations on the thermoelectric properties of PbSnTeSe. The transport coefficients reported by these groups are presented in Table 3. We first note that the experimental figures of merit ZT are lower than ours by a factor of around two. To trace back the origin of this difference, we compare the predicted transport coefficients at 700 K and for the carrier concentration in electrons of 6 × 10 19 e/cm −3 with the experimental results from Fan et al. The calculated total thermal conductivity (2.3 W/(mK)) is higher than that obtained experimentally, which should degrade the theoretical ZT value, but this is not what we observe. By contrast the calculated power factor S 2  is much higher (almost eight times as high). As our Seebeck coefficient is the same as the experimental one, the reason for the difference is to be found in the electrical conductivity. Indeed, theo amounts to 15 × 10 4 S/m. This large value can be explained by two factors. First, the calculated gap that includes the spin-orbit interaction is smaller by a factor of three than the

Comparison with Available Data on PbSnTeSe and Other HEA
As mentioned in the introduction section, Fan et al. [61] and Raphel et al. [62,63] have reported experimental investigations on the thermoelectric properties of PbSnTeSe. The transport coefficients reported by these groups are presented in Table 3. We first note that the experimental figures of merit ZT are lower than ours by a factor of around two. To trace back the origin of this difference, we compare the predicted transport coefficients at 700 K and for the carrier concentration in electrons of 6 × 10 19 e/cm −3 with the experimental results from Fan et al. The calculated total thermal conductivity (2.3 W/(mK)) is higher than that obtained experimentally, which should degrade the theoretical ZT value, but this is not what we observe. By contrast the calculated power factor S 2 σ is much higher (almost eight times as high). As our Seebeck coefficient is the same as the experimental one, the reason for the difference is to be found in the electrical conductivity. Indeed, σ theo amounts to 15 × 10 4 S/m. This large value can be explained by two factors. First, the calculated gap that includes the spin-orbit interaction is smaller by a factor of three than the experimental result, and second, we are modeling a pure, defect-free compound. In real compounds, electrons are scattered by impurities, defects, and grains boundaries, which are not accounted for in our model. Table 3. Experimental thermoelectric properties of PbSnTeSe (from Refs. [61,62]). The data from Fan et al. [61] are at 700 K and for n = 6 × 10 19 e/cm −3 , and those from Raphel et al. [62] are at 625 K.

Property
Fan et al. [ Recently, Bafekry et al. [85] have reported thermoelectric properties for the GeSnPb-SSeTe HEA. Their investigation is limited to the electronic transport properties (Seebeck coefficient and electrical and electronic thermal conductivities) but the electrons relaxation time was not determined. Interestingly, the Seebeck coefficient for the n-doped compound is very similar to that of PbSnTeSe (around −150 µV/K). Assuming about the same value of the electron relaxation time the electrical conductivity of GeSnPbSSeTe is of the same order of magnitude (~20 × 10 4 S/m) as that of PbSnTeSe. Assuming further that the thermal conductivity is similar for both compounds, one can infer that these compounds perform equally. From the experimental side [86], GeSnPbSSeTe shows a similar Seebeck coefficient, but the electrical conductivity is much lower, being of the order of 600 S/m.
As a different HEA, Sn 0.25 Pb 0.25 Mn 0.25 Ge 0.25 Te was investigated both experimentally and theoretically by Wang et al. [87]. A ZT value of 1.0 was found at 700 K, probably due to a drastic decrease in the thermal conductivity down to 0.76 W/(mK) by entropy engineering, compared to SnTe (4 W/(mK)). Compared to PbSnTeSe, the Seebeck coefficient of Sn 0.25 Pb 0.25 Mn 0.25 Ge 0.25 Te is twice as small (~100 µV/K), but the power factor is notably higher (14 × 10 −4 W/(mK 2 )). It was indeed observed that the electrical conductivity increases with alloying with more elements. This result shows that, PbSnTeSe could be an interesting candidate for thermoelectric application as a high-entropy alloy materials, but there is probably room for improvement, in particular on the electrical conductivity by further alloying with other elements.

Conclusions
In summary, by combining first-principles calculations and on-the-fly machine learning technique with the semiclassical Boltzmann transport theory and Green-Kubo theory, the thermoelectric transport properties of PbSnTeSe high-entropy alloy have been thoroughly investigated. The electronic and thermal transport coefficients of PbSnTeSe high entropy-alloy have been discussed in detail. The results indicate that PbSnTeSe has very low lattice thermal conductivities, below 0.4 W K −1 m −1 . It has been found that the PF values for n-type doping are always larger than those for p-type doping because of the higher electrical conductivity. The n-type PbSnTeSe exhibits better thermoelectric properties than the p-type. The maximum ZT (≈1.1) is found at 500 K for n-type doping. These results confirm that the PbSnTeSe HEA is a promising thermoelectric (TE) material.