Band degeneracy enhanced thermoelectric performance in layered oxyselenides by first-principles calculations

Band degeneracy is effective in optimizing the power factors of thermoelectric (TE) materials by enhancing the Seebeck coefficients. In this study, we demonstrate this effect in model systems of layered oxyselenide family by the density functional theory (DFT) combined with semi-classical Boltzmann transport theory. TE transport performance of layered LaCuOSe and BiCuOSe are fully compared. The results show that due to the larger electrical conductivities caused by longer electron relaxation times, the n-type systems show better TE performance than p-type systems for both LaCuOSe and BiCuOSe. Besides, the conduction band degeneracy of LaCuOSe leads to a larger Seebeck coefficient and a higher optimal carrier concentration than n-type BiCuOSe, and thus a higher power factor. The optimal figure of merit (ZT) value of 1.46 for n-type LaCuOSe is 22% larger than that of 1.2 for n-type BiCuOSe. This study highlights the potential of wide band gap material LaCuOSe for highly efficient TE applications, and demonstrates that inducing band degeneracy by cations substitution is an effective way to enhance the TE performance of layered oxyselenides.


INTRODUCTION
The thermoelectric (TE) effect can realize direct and reversible conversion between heat and electric energy, and thus provides a promising technology for green power generation and refrigeration systems 1 . The conversion efficiency of TE materials can be evaluated by the figure of merit ZT = S 2 σT/κ, in which, S is Seebeck coefficient, σ is electrical conductivity, T is absolute temperature 2 . The thermal conductivity κ consists of electronic thermal conductivity κ e and lattice thermal conductivity κ 1 . Generally, narrow band gap semiconductors are considered suitable for TE applications due to relatively high doping efficiency and excellent electronic transport properties 3,4 .
Recently, layered compounds show great promise in TE applications due to the intrinsic low lattice thermal conductivities with phonon anharmonicity induced by the interactions between adjacent layers 5,6 . Especially, the band structures of layered materials with multiple chemical bonding typically show multiband and/or multi-band valleys behaviors around the Fermi level, which provides the possibility of optimizing their TE transport properties by modifying electronic band structures, i.e., enhancing the power factors resulted from enlarged Seebeck coefficients 5,7 . As illustrated in Fig. 1, band valleys that next to the conduction band minimum (CBM) or valence band maximum (VBM), can be modify to move toward their respective band edges in higher energy direction through band engineering approaches, leading to a highly degenerated band structure with reduced energy difference. This band degeneracy helps to confine extra effective density of state (DOS) near the band edge region beyond the original effective DOS that only determined by VBM or CBM, resulting in enhanced Seebeck coefficient [7][8][9] . Therefore, layered semiconductors with narrow bandgaps and multiple chemical bonding have attracted wide attention in the quest for the TE materials 5,6,10 . Among them, p-type quaternary oxyselenides BiCuOSe (as shown in Fig. 2a, b) is an emerging TE material due to the moderate Seebeck coefficient, narrow band gap (~0.8 eV) 11 and low lattice thermal conductivity (~0.88 W m −1 K −1 ) 12 .
BiCuOSe has been reported as a potential TE material for the first time by Bleijs et al., and the ZT value of single phase BiCuOSe has been measured to be 0.53 at 873 K 12 . It is also found that the σ of BiCuOSe is found to be only 112 S m −1 at room temperature because of the weak interactions between the adjacent layers 13 . Such low σ drives the follow-up research works toward improving its electronic transport property. The most common way to increase the electrical conductivity of BiCuOSe is to dope divalent metal ions into its Bi atoms sites [14][15][16] . Recently, both experimental and theoretical investigations suggest n-type doping of BiCuOSe may achieve a larger σ and thus enhanced ZT value 5,17 . The optimal ZT value can be improved by 32% for n-type BiCuOSe and reach to 1.2 at 920 K, as compared with the p-type doping approach 17 . Yet, n-type BiCuOSe is shown to exhibit an intrinsically lower Seebeck coefficient than p-type BiCuOSe due to lower density of states in the conduction band minimum (CBM) than in the valence band maximum (VBM) 17 . From this perspective, single method of n-type doping is still restricted by low Seebeck coefficient, thus limiting further ZT optimization for this material.
Control of electronic structure, in this case the band degeneracy, might offer another approach to promote the thermoelectric properties of these layered oxyselenides. In BiCuOSe, CBM is occupied by Bi elements with 6 s character 18 , which intrinsically shows larger orbital bandwidth than that of d orbitals due to symmetry restriction [19][20][21] . Thus if more localized metal d-bands can be successfully introduced into the CBM of layered oxyselenides, it might be helpful to enhance the band-edge convergence thus the effective density of state at the CBM, leading to further improvement of Seebeck coefficient and PF value (Fig. 1b, c). Considering this, La element is better than Bi due to its significant contribution of 5d orbitals in the unoccupied states. Layered LaCuOSe, with similar isostructure of BiCuOSe (as shown in Fig. 2a, b), exhibits excellent photoelectric, magnetic and electronic transport properties [22][23][24][25] . It is proposed to be a candidate for TE materials since the work of Yasukawa et al., who found that its σ can be significantly enhanced by Sr-doping 26 . Follow-up studies also found that the carrier mobility, effective mass, and thermal conductivity of LaCuOSe are of anisotropic character due to its layered crystal structure 27,28 . However, study on the TE properties of LaCuOSe is rare, and there is no reliable report for ZT value of this crystal so far in whether the theoretical or the experimental studies.
Furthermore, the potential of LaCuOSe as a thermoelectric material is considered to be lower than BiCuOSe due to its wide band gap (~2.8 eV) 29 and higher lattice thermal conductivity (~2.1 W m −1 K −1 ) 26 than BiCuOSe. Recently, some wide band gap materials (1.5-3.5 eV) have been found to exhibit excellent TE performance due to high motilities and low lattice thermal conductivities [30][31][32] . For example, the σ of LaCuOSe has been reported to be 2560 S m −1 at 300 K 33 , which is more than 20 times larger than that of BiCuOSe. To make it clear, although BiCuOSe has higher S and lower κ 1 , LaCuOSe exhibits larger σ, thus leading to an open question that which one shows better overall TE performance? Existing literatures are mainly focused on the comparison of other properties, such as elastic 27 , thermal 34 and optoelectronic 18 properties. Direct comparison of electronic and thermal transport properties between BiCuOSe and LaCuOSe have not been reported yet. Hence, in order to reevaluate the potential of layered LaCuOSe for TE applications, it is necessary to understand the electronic structure and intrinsic TE transport properties of the two crystals.
In this work, the TE performance of LaCuOSe and BiCuOSe are systematically investigated by means of the first-principles calculations combined with Boltzmann transport theory. The calculated results show that as compared with BiCuOSe, La elements in LaCuOSe can effectively induce conduction band degeneracy and further enhance the Seebeck coefficient. As a result, the optimal ZT value of n-type LaCuOSe can reach to 1.46 at 900 K, which is larger than that of 1.2 for n-type BiCuOSe, demonstrating wide band gap semiconductor LaCuOSe exhibits better TE performance than BiCuOSe. This study unveils the critical role of La element in optimizing the electronic structure and further promoting thermoelectric properties of layered oxyselenides, which will provide useful information for further experimental and theoretical studies on these and other related layered TE materials.

RESULTS AND DISCUSSION
Crystal structure and electronic structure Both LaCuOSe and BiCuOSe crystals belong to tetragonal structure with P4/nmm space group (No.129) 29,35 . As shown in Fig. 2a, b, the anti-fluorite (Cu 2 Se 2 ) 2− chalcogenide layers and fluorite-like (La 2 O 2 ) 2+ /(Bi 2 O 2 ) 2+ oxide layers are stacked perpendicularly along c-axis. The optimized lattice parameters of LaCuOSe and BiCuOSe are in good agreement with the experimental results (see Table  1) 29,35 . On the other hand, the bond lengths and interlayer spacing for layered compounds provide information about the physical nature of the chemical bonds 27 . Our results show that the bond lengths and interlayer spacing for BiCuOSe are slightly shorter than those of LaCuOSe, which are consistent with previous reports 27, 34 . See Table 1 for more detailed comparison with other reported values.
The calculated electronic band structures and DOS distributions near the Fermi level of LaCuOSe and BiCuOSe are plotted in Fig. 2c-f. It shows that the difference of A-site cations affects the electronic structures of layered oxyselenides markedly. The band gap of BiCuOSe exhibits indirect character, i.e., the CBM is positioned at Z point and the VBM is on the M-G line, while LaCuOSe is exhibits direct gap with both CBM and VBM locate at Gamma (G) point. The band gaps of LaCuOSe and BiCuOSe are 2.59 eV and 0.87 eV without vdW, respectively, which are consistent with the previous experimental data of 2.8 eV for LaCuOSe and 0.8 eV for BiCuOSe 11,26 . It is found that the CBM of LaCuOSe is mostly contributed by the 5d orbitals of La atoms, and the 6p orbitals of Bi atoms contribute to the CBM of BiCuOSe (see Fig. 2e, f). It also shows that the VBMs of both LaCuOSe and BiCuOSe are mostly contributed by the 3d orbitals of Cu atoms and 4p orbitals of Se atoms.
It is noted that the energy difference between the lowest two conduction bands (i.e., E CB1 −E CBM ) of LaCuOSe at Gamma point is 0.1 eV, which is smaller than that of 0.19 eV for BiCuOSe at Z point (see Fig. 2c, d). As shown in Fig. 1, smaller energy difference between the nearest two bands around the Fermi level results in higher band degeneracy, leading to higher effective DOS for LaCuOSe than BiCuOSe in the n-type doing regime. In contrast, the energy difference between the two valence band valleys (one is at the Gamma point, and the other is on the M-G line) of BiCuOSe is almost 0 eV (i.e., E VBM −E VB1 ), which is much smaller than that of 0.15 eV for LaCuOSe, suggesting higher effective DOS for BiCuOSe than LaCuOSe in the p-type doing regime. The strategy of improving effective DOS by band degeneracy have also been reported in other TE materials with multi-band and/or multi-band valleys 10,36,37 .

Relaxation time and electronic transport property
Based on the optimized lattice structures of LaCuOSe and BiCuOSe, their electrical conductivities σ and relaxation time τ are studied by the Boltzmann transport equation with the BoltzTraP2 code. According to the deformation potential (DP) theory, the energy-dependent relaxation time is defined as 38 where r is related to specific scatting mechanism. Since the couplings of carriers and acoustic phonons for LaCuOSe and BiCuOSe are the predominant scattering, the r is set as −1/2 10 . The h, C ii , k B , m* and E 1 are the Planck constant, elastic constant, Boltzmann constant, effective mass and DP constant, respectively. The calculated C ii , m* and E 1 are listed in Table 2. Inserting these coefficients to Eq. (2), we obtain the relaxation times of BiCuOSe and LaCuOSe as a function of cattier concentration, as shown in Fig. 3a, b. The results show that, firstly, the relaxation time of BiCuOSe is always larger than that of LaCuOSe at given carrier type, carrier concentration, and temperature because of the smaller average effective mass. Secondly, the relaxation times of the two crystals are dependent on carrier types, i.e., the hole relaxation times are always smaller than the electron relaxation times at given temperature and carrier concentration due to the larger hole effective masses. Thirdly, the relaxation times for BiCuOSe and LaCuOSe are dependent on the temperature regardless of carrier types, i.e., the relaxation times decrease as temperature rises due to the more frequent scattering of carriers at the higher temperature, as is the case with layered GeAs 2 10 and GeSe 39 . For example, the relaxation times for n-type BiCuOSe and LaCuOSe vary over large range from 8.2 to 2.8 fs and 7.3 to 2.4 fs at the temperature of 300-900 K, respectively. The calculated hole relaxation time of BiCuOSe is close to the value reported by Yang et al. (1.26 fs at 780 K) 17 . However, there is no reliable report for relaxation time of LaCuOSe.
Inserting the relaxation time to σ/τ, the electrical conductivities for BiCuOSe and LaCuOSe are obtained, as shown in Fig. 3c, d. Generally, in all cases, the σ of both materials increase with the increase of carrier concentration. On one hand, the electrical conductivities of n-type BiCuOSe and LaCuOSe are always higher than those of p-type BiCuOSe and LaCuOSe because of lower scattering caused by longer electron relaxation times. We further plot the band decomposed charge densities around the CBM and VBM for both materials in Fig. 2g-j to explain the different transport behaviors of electrons and holes. Here, the isosurface level is set at 0.0005. For n-type BiCuOSe and LaCuOSe, all of the atoms are connected by charge densities, which can provide many conduction paths of electrons and effectively improve the electronic transport properties. For p-type BiCuOSe and LaCuOSe, there is no conduction channel can be found between atoms due to the strong localization of charge densities, which leads to weak charge transfer between atoms and thus reduces the electrical conductivities. Therefore, for BiCuOSe and LaCuOSe, the σ of n-type systems are most likely larger than those of p-type systems. On the other hand, regardless of carrier types, the electrical conductivity of BiCuOSe is higher than that of LaCuOSe at same temperature and carrier concentration due to the lower scattering caused by longer relaxation time. We note, however, experimental room temperature conductivity of BiCuOSe (112 S m −1 ) 13 is smaller than that of LaCuOSe (2560 S m −1 ) 33 , due to the lower intrinsic carrier concentration of BiCuOSe (~1 × 10 18 cm −3 ) than that of LaCuOSe (~2 × 10 19 cm −3 ). This can be understood from electronic structure point of view. Hole carrier concentrations of both materials can be calculated based on the Fermi distribution theory, i.e.,

Seebeck coefficient and power factor
The calculated Seebeck coefficients S are given in Fig. 4a, b. For a specific temperature, the absolute value of S for p-type (n-type) BiCuOSe is always larger (smaller) than LaCuOSe at the same carrier concentration, apart from the case of the S for p-type BiCuOSe at 900 K with the concentration below 1 × 10 19 cm −3 . For example, the Seebeck coefficients of p-type (n-type) LaCuOSe and BiCuOSe are 563 and 621 μV K −1 (−582 and −495 μV K −1 ) at 300 K with the carrier concentration of 1 × 10 18 cm −3 , respectively. In order to explain the difference in the Seebeck coefficients between LaCuOSe and BiCuOSe, the internal relationship between S and electronic structure is explored 7 , i.e., in which, the N, n, r, e, and k B are the effective DOS near Fermi level, carrier concentration, scattering mechanism parameter, electron charge and Boltzmann constant, respectively. According to the above equations, we can find that at given carrier concentration, the S is determined by N, which is affect by band degeneracy, as discussed in Fig. 1. This is also consistent with the charge density and DOS analysis. For example, Fig. 2e, f show that the slope of DOS around the CBM (VBM) for LaCuOSe is larger (smaller) than that of BiCuOSe, thus the absolute of S for n-type (ptype) LaCuOSe is larger (smaller) than that of BiCuOSe. Here, the large slope of DOS and Seebeck coefficient for p-type BiCuOSe are attributed to the band valley degeneracy around the CBM, and the band degeneracy around the VBM contributes to the large Seebeck coefficient of n-type LaCuOSe.
On the other band, the abnormal S-coefficient (greed dashed lines) for p-type BiCuOSe at 900 K with the carrier concentration form 1 × 10 18 cm −3 to 1 × 10 19 cm −3 can be attributed to bipolar conduction effect. The clear signature of bipolar conduction effect is that at high temperatures, the absolute value of Seebeck coefficient for narrow gap semiconductor increases first and then  The experimental data and theoretical values are also listed. (X = La and Bi). d is the interlayer spacing.  Fig. 4a. Power factor PF is an important parameter to measure the electronic transport property, which determines the coupling effects between S and σ, i.e., the σ increases as the carrier concentration rises, while the S is inversely proportional to carrier concentration. Hence, one can fit a balanced value of carrier concentration between S and σ to obtain the optimal power factor PF. The calculated power factors of LaCuOSe and BiCuOSe are plotted Fig. 4c, d. For both BiCuOSe and LaCuOSe, the optimal PF values for n-type systems are always higher than those for p-type systems due to the larger electrical conductivities. Regardless of carrier types, the optimal PF value decreases as the temperature rises because of the inverse correlation between electrical conductivity and temperature. More important is that the highest PF value of n-type LaCuOSe is 2.58 mW m −1 K −2 at 300 K, which is 1.5 times larger than that of 1.71 mW m −1 K −2 for n-type BiCuOSe. This value is also found larger than the optimal PF values of many layered TE materials at room temperature, such as Bi 2 O 2 Se (~0.4 mW m −1 K −2 ) 43 and Bi 2 O 2 Te (~0.55 mW m −1 K −2 ) 43 . Higher power factor of n-type LaCuOSe is mainly because the larger Seebeck coefficient caused by conduction band degeneracy and higher optimal carrier concentration (which corresponds to a higher electrical conductivity). For p-type systems, BiCuOSe shows higher optimal PF value than LaCuOSe under each considered temperature, which can be attribute to its larger Seebeck coefficient and higher electrical conductivity.

Thermal transport property
The electronic thermal conductivities κ e and lattice thermal conductivities κ 1 of LaCuOSe and BiCuOSe are further investigated. In this study, the electronic thermal conductivities for the two compounds are obtained by [44][45][46] κ e ¼ κ 0 À TσS 2 ¼ LσT (4) in which, κ 0 is the electronic thermal conductivity calculated under closed circuit condition, and L is Lorenz constant. The calculated results are given in Fig. 5a, b. It shows that the electronic thermal conductivities of n-type systems are always larger than p-type systems for both BiCuOSe and LaCuOSe. Besides, regardless of carrier types, the κ e of BiCuOSe is higher than that of LaCuOSe at same temperature and carrier concentration. The thermal transport behavior of carriers can be explained by that the κ e is linearly related to σ, and thus the variation of κ e follows the same trend of σ with the carrier concentration and temperature (see Fig. 3c, d).
Based on second and third-order interatomic force constants (IFCs), the lattice thermal conductivities of LaCuOSe and BiCuOSe are obtained, as shown in Fig. 5c, d. Obviously, the average lattice thermal conductivity of LaCuOSe is larger than that of BiCuOSe under each considered temperature. For example, the average lattice thermal conductivities of LaCuOSe and BiCuOSe are 1.73 W m −1 K −1 and 1.00 W m −1 K −1 at 300 K, respectively, which are close to the experiment values of 2.1 W m −1 K −1 for LaCuOSe 26 and 0.88 W m −1 K −1 for BiCuOSe 12 . Such low lattice thermal conductivities for LaCuOSe and BiCuOSe are comparable to those of many well-known TE materials, like PbTe and Bi 2 O 2 Se, for which the lattice thermal conductivities have been reported to bẽ 2.1 W m −1 K −1 47 and~1.1 W m −1 K −1 43 , respectively. Simultaneously, the lattice thermal conductivities of the two crystals are of anisotropic character. At room temperature, the lattice thermal conductivities along the in-plane and out-of-plane directions (κ 1a and κ 1c ) are 2.51 W m −1 K −1 and 0.18 W m −1 K −1 for LaCuOSe, 1.36 W m −1 K −1 and 0.32 W m −1 K −1 for BiCuOSe, respectively. The ratios κ 1a /κ 1c of LaCuOSe and BiCuOSe are 13.64 and 4.24 at 300 K, respectively, which are larger than those of 6.55 for GeAs 2 10 and 3.44 for phosphorene 48 .
The phonon spectra of LaCuOSe and BiCuOSe are plotted in Fig.  6a, b. It shows that the frequency gaps between the high and low frequency optical branches of LaCuOSe and BiCuOSe are 2.01 and 1.92 THz, respectively. The high frequency optical modes of LaCuOSe (BiCuOSe) are governed by O atoms, and the low frequency optical branches and acoustic branches are mainly dominated by La (Bi), Se and Cu atoms. Since the vibrational frequencies of Bi atoms are lower than those of La atoms, all of phonon modes of BiCuOSe move toward lower frequencies compared to LaCuOSe, implying that the phonon branches of BiCuOSe are more localized than LaCuOSe. According to the phonon spectra, the larger difference in the lattice thermal conductivities between LaCuOSe and BiCuOSe can be understand by the definition of lattice thermal conductivity 49 : in which, the C V , V g and l are the lattice heat capacity, phonon group velocity, and phonon mean-free path (PMFP), respectively. Within the harmonic approximation, the heat capacity can be given by 50 where, the C ph ½ω λ ðqÞ is the contribution of phonon modes for q points in the first Brillouin zone to the heat capacity. The λ is the index of the phonon modes. The calculated temperaturedependent heat capacities of LaCuOSe and BiCuOSe are given in Fig. 6c. It shows that the heat capacities of the two compounds show almost identical trends and values with temperature. Besides, their heat capacities exhibit the expected T 3 law at the low temperatures, and reach the Dulong-Petit classical limit at high temperatures (99.1 J mol −1 K −1 for BiCuOSe and 98.9 J mol −1 K −1 for LaCuOSe). Therefore, the larger difference in the lattice thermal conductivities between LaCuOSe and BiCuOSe is independent to the heat capacities of the two compounds. Our calculated result of the heat capacity for BiCuOSe is in good agreement with the result reported by Liu et al. 50 The phonon group velocities of acoustic phonon modes for LaCuOSe and BiCuOSe are given in Table 3. It can be seen that BiCuOSe exhibits lower acoustic phonon group velocity V g at given branch and orientation as compared with LaCuOSe. This is because the acoustic phonon group velocity is proportional to the slope of phonon modes at Gamma point, and the phonon branches of BiCuOSe are more localized than LaCuOSe. On the other hand, the cumulative lattice thermal conductivities as a function of PMFPs for the two crystals are plotted in Fig. 6d. It shows that the κ 1 increases rapidly with the increase of PMFP and then keeps stable when reach the upper limit of PMFP. The upper limit of PMFP for BiCuOSe (~80 nm) is lower than that of LaCuOSe (~150 nm) at room temperature. As compared with LaCuOSe, the lower V g and PMFP of BiCuOSe indicates more frequent scattering of phonons and lower heat transport property, which gives a reasonable explanation for the lower κ 1 of BiCuOSe than LaCuOSe. In contrast, the upper limits of PMFP are~300 nm for PbTe 47 and 20 nm for Bi 2 O 2 Te 43 , corresponding to the lattice thermal conductivities of~2.1 W m −1 K −1 and~0.57 W m −1 K −1 , respectively. We also note that the lattice thermal conductivities of LaCuOSe and BiCuOSe at 300 K will decrease 50% when the size of nanostructures below 4.1 nm and 2.0 nm, respectively. Hence, the low-dimensional LaCuOSe and BiCuOSe may exhibit better TE performance than bulk phase 49,51 .
The anisotropic lattice thermal conductivities for LaCuOSe and BiCuOSe are further explained. Generally, the thermal transport properties of layered crystals often are of anisotropic character due to the strong in-plane bonding and weak out-of-plane bonding, which can be reflected in the phonon group velocity V g and Debye temperature Θ D , i.e., κ 1 is positively correlated with Θ D and V g for a specific direction 10 . Table 3 shows large difference of Θ D and V g for BiCuOSe and LaCuOSe along the different directions, implying that lattice thermal conductivities of the Fig. 6 Analysis of lattice thermal transport property. The phonon spectra (a, b) for BiCuOSe and LaCuOSe. c The heat capacities C V for BiCuOSe and LaCuOSe as a function of temperature. d The normalized integration of lattice thermal conductivities with respect to phonon mean free paths (PMFP) at 300 K for BiCuOSe and LaCuOSe. Dimensionless figure of merit According to the calculated electronic and thermal transport coefficients, the ZT values of the two crystals are obtained, as shown in Fig. 7. Firstly, for both BiCuOSe and LaCuOSe, the optimal ZT values of n-type systems are always higher than those of p-type systems at a given temperature because of the higher electrical conductivities. As compared with the p-type doping approach, the maximum ZT values at 900 K could be improved 60% and 450% by using n-type doping for BiCuOSe and LaCuOSe, respectively. Secondly, for p-type systems, BiCuOSe exhibits higher optimal ZT value than LaCuOSe regardless of temperatures due to the larger Seebeck coefficient, higher electrical conductivity and lower lattice thermal conductivity. For instance, the optimal ZT value of 0.75 for BiCuOSe is 2.3 times than that of 0.32 for LaCuOSe at the temperature of 900 K. Thirdly, for n-type systems, the optimal ZT value of LaCuOSe is always larger than that of BiCuOSe at same temperature because of the larger Seebeck coefficient and higher optimal carrier concentration. For example, the optimal ZT value of LaCuOSe is 1.46 with the carrier concentration of~2.5 × 10 20 cm −3 at 900 K, which is about 22% larger than that of 1.2 with the carrier concentration of~1 × 10 20 cm −3 for BiCuOSe. The highest ZT value of LaCuOSe is also superior to that of many typical TE materials such as Bi 2 O 2 Se 52 and LaOBiSSe 53 , for which the ZT values have been reported to bẽ 0.95 at 800 K and~0.54 at 900 K, respectively. Therefore, the TE performance of LaCuOSe is better than that of BiCuOSe, and ntype doping can improve the ZT values of BiCuOSe and LaCuOSe remarkably.
We reiterate that for practical thermoelectric applications, both n-type and p-type components are very important for realization of complete thermoelectric modules toward large area, portable, high power and energy efficient applications. However, stable ntype doping for layered oxyselenides, such as LaCuOSe and BiCuOSe, are known to be difficulty due to the high volatility of more than one cations (Bi and Cu) under high-temperature materials synthesis. Therefore, previous literatures are primarily focused on hole doping because p-type BiCuOSe is commonly believed to show best thermoelectric performance among these materials. However, our work demonstrates that the n-type version of oxyselenides, although not easy to be synthesized, show unexpectedly high performance than the p-type, thus highlight the significances of n-type oxyselenides. We note that ntype doping of Cu-based oxyselenides have been recently experimentally synthesized by cation and anion codoping strategy 54 , thus providing illumination for n-type LaCuOSe. Additionally, from energy efficiency and cost-effective perspectives, the ZT values of materials should exceed 1.5 for practical application 55 . Thus, to fully realize the potential of LaCuOSe for thermoelectric energy harvesting, the working temperature is suggested to be on the order of 900 K.
Finally, we discuss the validity of models and approximations used in this calculation. (1) The DP theory considers the scattering between longitudinal acoustic phonons and carriers, which is proven as predominant scattering mechanism in BiCuOSe and LaCuOSe 33,56 , thus DP theory can provide relatively good approximation for their relaxation times and electrical conductivities. However, the phonon anharmonicity, i.e., the scattering among various phonons, is not taken into consideration in DP theory, which will become stronger as temperature increases. Thus, the relaxation times and ZT values of BiCuOSe and LaCuOSe may be underestimated at high temperatures, due to the weakening of phonon-carrier scattering caused by strong phonon-phonon interactions 57 . Unfortunately, there is no reliable theory to explicitly evaluate the effect of anharmonicity on carrier scattering (the difficulty of many-body interactions) 58,59 , thus DP theory is generally widely adopted to approximately evaluate the relaxation times of TE materials (including BiCuOSe-type materials) at high temperatures 10,60,61 . Meanwhile, we also note that the contribution of phonon anharmonicity of LaCuOSe and BiCuOSe are comparable because of similar Grüneisen parameters, which are 1.93 and 1.84 for BiCuOSe and LaCuOSe, respectively 27 . Hence, the neglection of phonon anharmonicity effect is not expected to significantly change the comparison of the thermoelectric performance of LaCuOSe and BiCuOSe. (2) The carrier concentration can also impact scattering mechanism. We analyze carrier transport based on single parabolic band (SPB) model, in which the m* is calculated straightforwardly by band curvatures at CBM/ VBM. While the effects of inter-pocket scattering are generally neglected, because of their relatively weak contributions at normal carrier concentration levels. Under this approximation, DP theory with SPB model is shown to yield a good prediction for relaxation time and thermoelectric properties in such conditions 55,62 and it is thus also widely used for materials with multi- band and/or multi-band valleys 10,63 . Strictly speaking, for materials with multi-band and/or multi-band valley behaviors such as LaCuOSe and BiCuOSe, electrons (holes) will first occupy the energy pocket right at the band edge, e.g. VBM or CBM. Then for the case of higher carrier concentrations, the extra carriers will continue to fill the orbitals next to the level of VBM/CBM. The m* of carriers in different energy pockets is different, and carriers transition between different energy pockets will cause inter-picket scattering, thus the average m* and carrier mobility of carriers will be changed at high doping level 28 . Therefore, SPB model may overestimate the ZT values of LaCuOSe and BiCuOSe at high carrier concentrations. However, in our work, the discussed carrier concentration is not so high, therefore, it is still valid to use SPB model to compare the thermoelectric properties of LaCuOSe and BiCuOSe. (3) The maximum ZT value of p-type BiCuOSe in this work is less than the experimentally reported value of 1.5 at 873 K 64 . The difference between our calculated results and experimental values comes from the following two aspects: on one hand, electrical conductivity is calculated by the semi-classical Boltzmann transport theory, which only considers doped carriers with associated scattering and neglects the effects of remaining ions in the host lattice (ionization scattering) on the crystal and electronic structures. On the other hand, element doping induced phonon anharmonicity can lead to extra phonon scattering, thus screening lattice thermal conductivity. While this can be easily measured experimentally, it is computationally prohibitive to calculate by current transport theory under Boltzmann framework. Thus, theoretical simulations based on ideal crystal with carrier doping tends to overestimate κ 1 65 . Our obtained results are comparable to other theoretical reported ZT value 17,52 .
In conclusion, we systematically investigate and compare the intrinsic TE performance of LaCuOSe and BiCuOSe by utilizing firstprinciples calculations in this paper. The calculated results show that because of the different A-site cations between BiCuOSe and LaCuOSe, their electronic structures show significant difference. In addition to the wider band gap of LaCuOSe than BiCuOSe, the most important thing is that band valley degeneracy of BiCuOSe is found around the VBM, while LaCuOSe shows band degeneracy behavior around the CBM. Band degeneracy can significant enhance the Seebeck coefficient, which leads to that the Scoefficient of p-type (n-type) BiCuOSe is larger (smaller) than that of LaCuOSe at same temperature and carrier concentration. Due to the longer electron relaxation times, the electrical conductivities of both BiCuOSe and LaCuOSe for n-type systems are higher than those for p-type systems. A large power factor of 2.58 mW m −1 K −2 for n-type LaCuOSe can be achieved at 900 K because of the larger Seebeck coefficient and higher optimal carrier concentration, which is 1.5 times larger than that of 1.71 mW m −1 K −2 for n-type BiCuOSe. Besides, the average lattice thermal conductivity of BiOCuSe (1.00 W m −1 K −1 at 300 K) is lower than that of LaCuOSe (1.73 W m −1 K −1 at 300 K) due to lower phonon group velocity and phonon mean-free path. Under the co-effect of those transport coefficients, the optimal ZT value of 1.46 for LaCuOSe is obtained with the carrier concentration of~2.5 × 10 20 cm −3 at 900 K, which is about 22% larger than that of 1.2 for BiCuOSe with the carrier concentration of~1 × 10 20 cm −3 . The present study thus demonstrates the critical role of La element in promoting TE performance of layered oxyselenides, which will provide beneficial information for further related studies.

First principle calculations
Our calculations are performed within DFT by using the Vienna Ab Initio Simulation Package (VASP) [66][67][68] . The generalized gradient approximations (GGA) 69 with the Perdew-Burke-Ernzerhof (PBE) parameterization is applied as the exchange-correlation potential. The Heyd-Scuseria-Ernzerhof (HSE06) hybrid density functional 70 is also employed to obtain accurate band structures and transport properties. The Brillouin zones of the primitive cells are sampled with Monkhorst-Pack k-meshes of 12 × 12 × 6 for relaxation and self-consistent calculation. The energy and the force convergence criterion are chosen to be 10 −4 eV and 10 −5 eV Å −1 , respectively. The cutoff energy of the plane-wave is set at 400 eV. The electronic transport properties are implemented in the BoltzTraP2 code 71 by utilizing Boltzmann transport theory and relaxation time approximation (RTA). The phonon structure is obtained by utilizing the Phonopy package 72 . The lattice thermal conductivity is calculated by solving the Boltzmann transport equation as implemented in ShengBTE code 65 . The second and third-order interatomic force constants (IFCs) are obtained using 4 × 4 × 1 and 3 × 3 × 1 supercells as inputs of Phonopy and thirdorder.py script 65 . The 10th nearest neighbors are chosen as the interaction cutoff for the third-order IFCs. The effects of SOC are considered for all calculations.

Van der Waals interactions test
For van der Waals (vdW) systems, the semi-empirical schemes like DFT-D3 method have been widely used to consider vdW interaction [73][74][75] . Besides the semi-empirical schemes, the vdW treatments with the exchangecorrelation functional correction like optB86b-vdW are also effective methods [76][77][78] . To investigate the influences of vdW interaction on the calculated results, we compare the lattice parameters, band structure, and phonon spectra for LaCuOSe and BiCuOSe obtained by standard DFT, DFT-D3, and optB86b-vdW methods in Table 1 and Figs. 8, 9. For LaCuOSe, it is shown that the lattice parameters obtained by DFT method with and without vdW correction are all in reasonable agreement with experimental values 29 , as shown in Table 1. Particularly, the band structure (and phonon spectra) show similar character to each other for standard DFT, DFT-D3 and optB86b-vdW calculations (see Figs. 8,9). In the case of BiCuOSe, we also find that consideration of vdW correction has small influences on the band structure and phonon spectra. The vdW correction, thus, is not taken into account in the subsequent calculations.

DATA AVAILABILITY
The data that support the findings of this study are available from corresponding author upon reasonable request. Fig. 9 Effects of vdW correction on phonon structure. The calculated phonon spectra for BiCuOSe (a-c) and LaCuOSe (d-f) with and without vdW correction.