## 2 Methods and Approaches  
### Calculation Details  
The ab initio study of van der Waals heterostructures was carried out within the framework of the density functional theory (DFT) implemented in the Siesta 4.1.5 software package [37; 38; 39]. To describe exchange-correlation effects, generalized gradient approximation (GGA) was used in the parametrization of Perdew, Burke and Ernzerhof (PBE) [40]. To describe the van der Waals interaction (vdW) between heterostructure layers, we used the correction scheme proposed by [PERSON], where vdW interactions are represented by a pair force field DFT + D2 [41]. When optimizing the geometry of the structure, we used the basis set of split valence orbitals DZ2P (double zeta 2 polarization), which includes polarization functions. The integration of the Brillouin zone was carried out using the Monkhorst-Pack scheme [42] with a 12 \(\times\) 6 \(\times\) 1 k-points mesh. The relaxation of the structure was carried out until the maximum interatomic force was less than 0.025 eV/A. The effective Broyden-Pulay mixing scheme was used to minimize the energy of the electronic subsystem [43]. In order to avoid interaction between neighboring structures in calculations with periodic boundary conditions, the translation vector along the z-axis of the supercell was set to be greater than 20 A. The real space mesh cutoff was chosen to be 300 Ry throughout the calculations. To study the optical properties of van der Waals heterostructures, the complex frequency-dependent dielectric function was calculated in SIESTA software package. The calculations were carried out in the framework of the first-order time-dependent perturbation theory [44]. For this purpose, the self-consistent ground state energies and eigenfunctions were first calculated, and then they were used to calculate the transition dipole moment matrix elements. The imaginary part of the dielectric function determined the optical absorption, which was calculated from the transition rate between valance and conduction band states. The equation for calculating the absorption coefficient is as follows:  
\[\alpha(\omega)=\frac{\omega}{cn(\omega)}\varepsilon_{2}(\omega), \tag{1}\]  
where \(\varepsilon_{2}(\omega)\) is the imaginary part of the complex dielectric function, \(n(\omega)\) is the refractive index, \(c\) is the speed of light. The calculation of \(\varepsilon_{2}(\omega)\) was carried out according to the following equation:  
\[\varepsilon_{2}(\omega)=\frac{e^{2}}{\pi m^{2}\omega^{2}}\sum_{\omega,c}\int_{ BZ}d\vec{k}\Big{|}\Big{\langle}\psi_{ck}\Big{|}\hat{e},\overrightarrow{p} \Big{|}\psi_{rk}\Big{\rangle}\Big{|}^{2}\delta(E_{c}(k)-E_{v}(k)-\hbar\omega), \tag{2}\]  
where the summation was carried out over each pair of states in the valence band (filled) and the conduction band (unoccupied), and integration was performed over all k-points in the Brillouin zone; the indices \(c\) and \(\
u\) refer to the states in the conduction and valence bands, respectively, \(E_{(c,\
u)}\) (\(k\)) and \(\psi_{(c,\
u),k}\) are the energy and eigenfunction of these states, respectively. The electronic dipole transition matrix element is between the occupied and unoccupied states, where \(\hat{e}\) is the polarization vector and \(\overrightarrow{p}\) is the momentum operator. The optical characteristics of the van der Waals heterostructures were calculated for two different directions of light polarization (vector \(\mathbf{E}\) is parallel to the x-axis; vector \(\mathbf{E}\) is perpendicular to the x-axis) for the energy range from 0.6 eV to 20 eV using optical broadening of 0.05 eV. The optical characteristics of the heterostructures under study were calculated by integrating the Brillouin zone with a 114 \(\times\) 65 \(\times\) 1 k-points mesh. Based on the absorption spectrum, the photocurrent spectrum was calculated.