BiOCuSe is formed by ionic layers of composition [Bi2O2]2+ and [Cu2Se2]2−, and adopts a tetragonal cell with space group *P4/nmm* (Fig. 1). The optimized lattice constants in our calculations (a = 3.926 Å, c = 8.957 Å) show good agreement with experimental values (a = 3.930 Å, c = 8.927 Å).9 The electronic band structures and density of states (DOS) calculated with the PBE + U + SOC functional are shown in Fig. 1a. The calculated band gap is 0.37 eV. Previous HSE + SOC calculations show a band gap of 0.81 eV, which is closer to the experimental measure of 0.80 eV.9,23 However, our band structure compares well with that calculated with the more expensive HSE + SOC functional after a rigid band shift of 0.44 eV.23,29 The top of the valence band of BiOCuSe is composed of Se and Cu orbitals and the bottom of the conduction band is constituted of Bi orbitals. The valence band maximum (VBM) is located along the M-Γ direction. The charge density at the VBM is shown in Fig. 1c. There are three valence band maxima, one along Γ-M, one at Z, and one along Z-R with only a minor energy difference (15 meV between VBM and the valence band maxima at Z, 51 meV between VBM and the valence band maxima along Z-R). The hole effective mass around the three valence band maxima is calculated to be 0.21 *m*0, 0.68 *m*0, and 4.60 *m*0, (where *m*0 is the electron mass), respectively, which agree with previous theoretical results (Γ-M: 0.23, Z-Γ: 0.47, and R-Z: 4.53).23

The phonon band structure and phonon DOS of BiOCuSe are shown in Fig. 1b. The phonons with frequencies above 200 cm− 1 are mainly projected to oxygen and those below 200 cm− 1 are mostly due to the heavier atoms. The Γ point has D4h point group symmetry, with phonons belonging to the following irreducible representations: Γ = 4A2u + 2A1g + 2B1g + 4Eu + 4Eg. For the doubly degenerate Eu and Eg modes, atoms vibrate in plane, while for A2u, A1g, B1g, and B2u modes, atoms vibrate out of plane. Figure 1e shows the six PO modes of BiOCuSe. The calculated frequencies of the PO modes compare well with literature (Table S1).30 The lower four PO phonon modes with frequencies below 200 cm− 1 are mainly due to the atomic displacements in the Cu-Se layer, while the two PO modes with frequencies above 200 cm− 1 are due to oxygen vibrations.

In order to discriminate the relative importance of LA and PO phonon scattering on charge transport in BiOCuSe, the LA phonon scattering is calculated by using deformation potential (DP) theory31 and the PO phonon scattering is calculated by accounting for the Fröhlich interaction.27 The DP constants are calculated by monitoring the variation of the energy of VBM with strain (Fig. 1d). Although this method is approximate,32 it has previously been successfully applied to describe the LA phonon scattering of CH3NH3PbI3, and the calculated DP constants (4.3 eV for electrons and 1.5 eV for holes)33 are very close to experimental results (2.93 for electrons, 2.2 eV for holes).34 The DP constants for holes in BiOCuSe are calculated to be 4.5 and 4.1 eV in the *a* and *c* directions, respectively (Table 1). We see no reason why the DP constant of holes in BiOCuSe should be as high as 24 eV,23,24 except for fitting experimental results. The elastic constants of BiOCuSe are calculated to be 155.1 and 103.7 GPa in a and c directions, respectively. The Young’s modulus is 104.1 and 63.5 GPa in the a and **c** directions, respectively, which agree well with the experimental value of 76.5 GPa measured on polycrystalline samples.35

Table 1

Deformation potential *E*1, elastic constant *C*ii, Young’s Modulus, LA phonon scattering limited relaxation time *τ*LA, and mobility *µ*LA, PO phonon scattering limited relaxation time *τ*PO and mobility *µ*PO, and the total carrier mobility limited by LA and PO phonon scattering *µ*Total for holes in BiOCuSe at 300 K.

| This work | Experiment |

direction | a | c |

*E*1 (eV) | 4.5 | 4.1 | |

*C**ii* (GPa) | 155.1 | 103.7 | |

Young’s Modulus (GPa) | 104.1 | 63.5 | 76.535 |

*τ*LA (fs) | 80 | |

*τ*PO (fs) | 15.0 | |

*µ*LA (cm2 V− 1 s− 1) | 168 | 64.7 | |

*µ*PO (cm2 V− 1 s− 1) | 45.1 | 17.7 | |

*µ*Total (cm2 V− 1 s− 1) | 31.6 | 12.4 | 2015 |

Fig. 2a and 2b show the temperature dependence of hole mobility calculated considering both LA and PO phonon scattering, and the limiting contributions when one only of the two terms is accounted for. The calculated room temperature mobilities in the *a *and *c* directions are 31.6 cm2 V-1 s-1 and 12.4 cm2 V-1 s-1, close to the experimental value (~ 20 cm2 V-1 s-1) for pristine BiOCuSe measured on polycrystalline samples which is thus reported as isotropic.15 We attribute the mobility anisotropy between in-plane and out-of-plane directions to the anisotropy of effective mass (0.21 *m*0, in plane vs. 0.68* m*0 out-of-plane), clearly linked with the layered structure of the compound. Our results indicate that the PO dominates over the LA phonon scattering for both in-plane and out-of-plane directions, in particular at low temperature; only for temperatures above 600K does the LA scattering become the predominant mobility-limiting mechanism. The combined scattering of LA and PO phonons results in the power law of *T*-1.5 and *T*-1.6 in *a* and *c* directions respectively, in good agreement with the power low from experiment.15 The LA phonon scattering limited mobility (i.e. the mobility obtained neglecting PO scattering) follows a power law of *T*-2.4. The assumption made until now when investigating mobility in BiOCuSe that LA phonon scattering follows a power law of ~ *T*-1.5,15,20,24 is therefore incorrect. This can be attributed to the multivalley valence bands of BiOCuSe, as occurs for InTe, where LA phonon scattering coupled with intervalley scattering results in a power law of *T*-2.01.36 A valley degeneracy of eight is found in previous theoretical work and confirmed by our band structure calculations shown in Figure 1.23 In our analysis, the multiple valley scattering effects are included by calculating the scattering between electron states with the same energy in different valleys.

We further separated the effects of PO modes above and below 200 cm− 1 to discriminate their relative contribution to hole mobility. As shown in Fig. 2c and 2d, the heavier atoms in the Cu-Se layer (modes below 200 cm− 1) play a dominant role at *T* below ~ 500 K, while oxygen vibrations become more important on increasing *T* above ~ 500 K. It is interesting to note that increased carrier mobility has been reported experimentally upon Te doping into the Cu-Se layer, which weakens the electron-phonon scattering through higher bond covalency.22 A temperature power law of *T*− 1.9 is obtained for the Te doped samples.22 This is consistent with our simulation, because when the polar optical phonon scattering from the Cu-Se layer is weakened, the power law will increase towards that of LA phonons.

Fig. 3 shows the dependence of the thermoelectric coefficients of BiOCuSe on carrier concentration at 300 K for the in-plane direction. For simplicity, we only discuss results in the in-plane direction in the following. Both LA and PO scattering mechanisms are included. The Seebeck coefficient *S* is strongly affected by the band gap size due to the possible bipolar effect. Therefore, we applied a rigid band shift (0.44 eV) on the conduction bands to match the gap of 0.81eV from experiment and HSE06+SOC calculations. Fig. 3a shows that *S* decreases with increasing carrier concentration and has excellent quantitative match with experimental results, especially at higher carrier concentrations. Fig. 3b and 3c show that both electrical conductivity *s* and electronic thermal conductivity *k*e increase linearly with carrier concentration, and again agree exceptionally with experiments at low carrier concentrations < 1019 cm-3. However, *s* and *k*e are overestimated at elevated carrier concentration due to the exclusion of extrinsic scattering mechanisms such as defect and boundary scattering in our simulations. The Wiedemann-Franz law connects σ and *k*e by *k*e = *LσT*, where L is the Lorenz number. The Lorenz number ranges between 1.53 × 10-8 and 2.38 × 10-8 W Ω K-2 when carrier concentration changes from 1018 cm-3 to 2 × 1021 cm-3, and compares well with the Mg doped BiOCuSe at low carrier concentration (Fig. 3d). It becomes closer to the theoretical value for metals (2.44 × 10-8 W Ω K-2) at high carrier concentrations. As shown in Fig. 3e, the theoretical limit of the power factor is 21.3 μW cm-1 K-2 at a carrier concentration of 2.0 × 1020 cm-3, which is approximately double the values currently achieved experimentally in various doped systems. The experimental lattice thermal conductivity of polycrystalline BiOCuSe samples (Fig. S1) is used to calculate the thermoelectric figure of merit of BiOCuSe.18 The optimal room temperature *ZT* is calculated to be 0.48 at 1.0 × 1020 cm-3. As shown in Fig. 3f, most of the experimental *ZT* values are below the theoretical values presented in this study. Some however become larger than the theoretical value at higher carrier concentrations. The possible reason is that microstructural effects, such as creation of nanodots19 and all-scale structural optimization21 are employed in experiments to reduce the lattice thermal conductivity below the polycrystalline value used here.18