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

CONTEXT: 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 ∈ [ 0 , ∞ [ , and the separation Δ 12 = 127.8 meV reproduces the experimental value without adjusted parameters. We further derive a generalized two-sided dielectric potential V ( r , N top , N bot ) that captures the asymmetric encapsulation geometry, and compute the two-dimensional map E b ( N top , N 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 ( 2 s ) / E b ( 1 s ) 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. METHODS: 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, ξ 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.