Finite-difference analysis of excitonic spectra in monolayer transition-metal dichalcogenides/hBN heterostructures

We theoretically study the excitonic properties in semiconductor transition metal dichalcogenide monolayers (ML-TMDs)(WSe $_2^{}$ , WS $_2^{}$ , MoSe $_2^{}$ , MoS $_2^{}$ , and MoTe $_2^{}$ ) encapsulated in hexagonal boron nitride (hBN). Excitonic binding energies and full non-hydrogenic Rydberg series (s, p, and d states) are calculated. As a first validation, the binding energies of the 1s-5s states of WSe $_2^{}$ /hBN are reproduced to less than 2 meV ( $<1\%$ ) over the entire range of thicknesses $N\in [0,\infty [$ , and the separation ${\Delta }_{12}=127.8$  meV reproduces the experimental value without adjusted parameters. We further derive a generalized two-sided dielectric potential $V(r,N_{\text{top}},N_{\text{bot}})$ that captures the asymmetric encapsulation geometry, and compute the two-dimensional map $E_{\text{b}}(N_{\text{top}},N_{\text{bot}})$ for WSe $_2^{}$ . A new analysis of Coulomb engineering shows that E $_b^{}$ (1s) varies by up to a factor of three across experimentally accessible dielectric environments, while the ratio $E_b\left(2s\right)/E^b\left(1s\right)$ serves as a sensitive probe of the non-hydrogenic character. Moreover, Berry-curvature corrections lift the angular-momentum degeneracy of p-and d-type excitons by 2-8 meV, with M-based compounds showing larger effects. Finally, the quadratic Stark polarizability under an out-of-plane electric field is computed via exact $s-p$ coupling, with results depending sensitively on N and the material parameters. These results provide a unified quantitative framework directly applicable to studies on van der Waals heterostructures.

The generalized Rytova-Keldysh (RK) Schrödinger equation is solved using a cell-centered finite-difference scheme on a fine real-space grid (N = 2000 points, ${\xi }_{\max }$ = 50), with the calculations incorporating both Coulomb engineering effects and an applied out-of-plane electric field. Topological Berry phase corrections are subsequently introduced: the differential Berry curvature between the conduction and valence bands computed within the bulk Dirac model is averaged over the exciton momentum distribution through a Bessel-Hankel transform, yielding an energy shift that scales linearly with the azimuthal quantum number m.