# Continuum model study of optical absorption by hybridized moiré excitons in
transition metal dichalcogenide heterobilayers

Yao-Wen Chang  Physics Division, National Center for Theoretical Sciences,
National Taiwan University, Taipei 10617, Taiwan

###### Abstract

We propose a continuum model for the theoretical study of hybridized moiré
excitons in transition metal dichalcogenides heterobilayers, and we use a
variational method to solve the exciton wavefunction and calculate the optical
absorption spectrum. The exciton continuum model is built by the charge
continuum model for electrons and holes in moiré superlattices, thereby
preserving the moiré periodicity and lattice symmetry from the charge
continuum model. The momentum-space shift of interlayer electron-hole
distribution is included, and thus the indirect nature of interlayer excitons
is described. The spin and valley degrees of freedom and related interactions
are omitted, except for the spin-orbit energy splitting of A and B excitons.
This continuum model is applied to the simulation of optical absorption by
hybridized moiré excitons in $\text{WSe}_{2}$/$\text{WS}_{2}$ and
$\text{MoSe}_{2}$/$\text{WS}_{2}$ heterobilayers. Twist-angle and electric-
field dependences of absorption spectra are studied. Calculated spectra are
compared with experimental observations in the literature, and correspondences
of signatures are found. The deficiency and the potential of the present model
are discussed.

##  I Introduction

Two-dimensional (2D) semiconducting materials, such as transition metal
dichalcogenides (TMDCs), are promising to be applied in future photonic and
optoelectronic deviceswang2012electronics ; xia2014two ; mak2016photonics .
Excitonic effects are critical for deciding the optical properties of 2D
semiconductors due to the reduction of dielectric screening and the
enhancement of quantum confinementxiao2017excitons ; berkelbach2017optical .
Recently, moiré heterostructures stacked by 2D semiconducting materials with
twisted angles or lattice mismatches have been predicted theoretically and
synthesized experimentally. Emergent physical phenomena related to the
additional moiré degree of freedom have been observed and
discussedtran2020moire ; kennes2021moire ; mak2022semiconductor ;
zhang2023controlling ; du2023moire . Exciton physics in this system is also
studied. Excitons in moiré heterostructures are modulated by interlayer
charge-transfer couplings and moiré periodic potentials, which are formed by
local atomic registry changeskunstmann2018momentum ; rivera2018interlayer ;
alexeev2019resonantly ; jin2019observation ; jin2019identification ;
tran2019evidence ; deilmann2020light ; shimazaki2020strongly ; zhang2020twist
; tang2021tuning ; jiang2021interlayer ; liu2021moire ; wilson2021excitons ;
huang2022excitons ; lin2023remarkably . Optical properties of excitons in 2D
heterostructures can be tuned by this additional moiré degree of freedom. The
applications of 2D layered materials are widened by this new tunability.

An exciton is an electron-hole pair bound by Coulomb
attractionxiao2017excitons ; berkelbach2017optical . Excitonic signatures can
be found in direct-gap semiconductors and detected by optical spectroscopy.
For multilayer materials, excitons can be sorted as intralayer excitons, where
the electron and the hole are at the same layer, and interlayer excitons,
where the electron and the hole are at different layers. In TMDC
heterobilayers with type-II band alignment, the lowest-energy excitation is
contributed from the interlayer exciton, which has a small electron-hole
wavefunction overlaprivera2018interlayer ; deilmann2020light . Therefore,
interlayer valence-to-conduction band transition has a low oscillator
strength, and thus interlayer-exciton signatures are difficult to be observed
in absorption spectra. While there is no direct interlayer exciton transition
in TMDC heterobilayers, electron transfer between conduction bands on each
layer or hole transfer between valence bands on each layer could occur. The
charge transfer contributes to the formation of interlayer excitons. In some
situations, such as by applying a perpendicularly electric field to tune the
band gaps, the difference between the optical transition energy of interlayer
excitons and the transition energy of intralayer excitons can be at the same
scale as the interlayer charge-transfer coupling. Hybridization between
interlayer excitons and intralayer excitons forms hybrid excitons. These
hybrid excitons borrow the oscillator strength from the intralayer exciton,
and thus they can be observed in absorption spectrachaves2018electrical ;
hsu2019tailoring ; hagel2021exciton .

With the additional moiré degree of freedom in TMDC heterobilayers, excitons
become trapped by the moiré potential, and the translational symmetry of the
center-of-mass (COM) motion of excitons becomes broken. The interlayer charge-
transfer coupling should also follow the same moiré periodicity with the
potential. If the physical properties of an exciton are decided or strongly
affected by this moiré degree of freedom, the exciton can be called a moiré
excitonliu2021moire ; wilson2021excitons ; huang2022excitons . Various
theoretical methods have been proposed and applied to the research of moiré
excitonsyu2015anomalous ; yu2017moire ; wu2017topological ; wu2018theory ;
danovich2018localized ; ruiz2019interlayer ; ruiz2020theory ;
brem2020hybridized ; brem2020tunable ; guo2020shedding ; choi2021twist ;
hichri2021resonance ; yu2021luminescence ; julku2022nonlocal ;
naik2022intralayer ; hagel2022electrical . An effective continuum model which
describes the COM motion of an exciton being modulated by a periodic potential
as the moiré potential has been widely used to find exciton band structures
and exciton wavefunctionsyu2017moire ; wu2017topological ; wu2018theory . For
more generalized continuum modelsruiz2019interlayer ; choi2021twist ;
yu2021luminescence ; julku2022nonlocal , the interlayer charge-transfer
coupling is included. However, these models only describe the COM motion of
moiré excitons, such that the effect of moiré periodicity on the internal
motion of excitons has not been discussed. A more sophisticated method is to
apply atomistic calculation with density functional theory (DFT) to simulate
the hybridized moiré excitonsbrem2020hybridized ; brem2020tunable ;
guo2020shedding ; naik2022intralayer . However, since the unit cell for moiré
superlattices is usually quite large, the atomistic approach with DFT is
expensive computationally. A method to model impartially and simulate
efficiently hybridized moiré excitons is still lacking. It is the gap we
intend to fill in this work.

In this article, we provide an improved version of the exciton continuum model
to study the hybridized moiré excitons in TMDC heterobilayers with lattice
mismatches and small twist angles. The exciton continuum model is written down
based on four conditions. First, the exciton continuum model is built by the
charge continuum model for electrons and holes in moiré superlattices, thereby
preserving the moiré periodicity and lattice symmetry from the charge
continuum model. Second, the momentum-space shift of interlayer electron-hole
distribution is included, and thus the indirect nature of interlayer excitons
is described. Third, the spin and valley degrees of freedom and related
interactions are omitted, except for the spin-orbit energy splitting of A and
B excitons. Forth, in long moiré-wavelength and zero charge-transfer-coupling
limits, the exciton model and the optical absorption formula can be reduced to
the counterparts of an isolated exciton. A variational method using Slater-
type orbitals (STOs) as the basis function is applied to the numerical
solution of the exciton continuum model. Optical absorption by hybridized
moiré excitons is studied with varying twist angles in small degrees and
varying out-of-plane electric fields to tune the band-edge energy. Twist-
angle-dependent and electric-field-dependent absorption spectra of
$\text{WSe}_{2}$/$\text{WS}_{2}$ and $\text{MoSe}_{2}$/$\text{WS}_{2}$
heterobilayers are simulated. A good correspondence between the calculated
spectra from this work and the experimentally observed spectra from the
literature is achieved. In Sec. II, the exciton continuum model is derived
based on the moiré periodicity and the charge continuum model. The methods to
solve the exciton wavefunction and the optical transition amplitude are also
introduced. In Sec. III, the model is applied to the simulation of twist-
angle-dependent and electric-field-dependent absorption spectra of
$\text{WSe}_{2}$/$\text{WS}_{2}$ and $\text{MoSe}_{2}$/$\text{WS}_{2}$
heterobilayers. In Sec. IV, the applications and improvements of this model in
the future study are discussed, and the conclusion is given. In Appendix B,
matrix elements of the exciton continuum model spanned by STO and plane-wave
basis functions are given.

##  II Theory

In this section, the continuum model for hybridized moiré excitons in TMDC
heterobilayers is derived, and the method to calculate the exciton
wavefunction and the optical spectrum is given. In Sec. II.1, the geometry of
moiré superlattices and reciprocal lattices is illustrated, and the concept of
moiré Brillouin zone (MBZ) is introduced. In Sec. II.2, the charge continuum
model and the many-body formulation are introduced. The expression of screened
Coulomb interaction in the many-body Hamiltonian is also discussed. In Sec.
II.3, the exciton continuum model is derived through the many-body
formulation. A transformation from electron and hole position coordinates to
exciton COM and internal coordinates is applied to the continuum model. In
Sec. II.4, the variational method is introduced to solve the exciton
wavefunction. In Sec. II.5, the method to calculate the optical absorption
spectrum is discussed.

###  II.1 Moiré superlattice and MBZ

Figure 1: (a) Schematic plot of the R-stacked heterobilayer with lattice
mismatch $\delta=0.1$ and a twisted angle $\theta=3^{\circ}$. (b) Schematic
plot of the Brillouin zone of the top monolayer (inside blue dash-dot
hexagonal) and the Brillouin zone of the bottom monolayer (inside orange dot
hexagonal) of the R-stacked heterobilayer. The filled zone is the MBZ of
heterobilayers. (c) Energy-level diagram of band edges at the top monolayer
$\mathbf{K}_{1}$, $\mathbf{K}^{\prime}_{1}$ and band edges at the bottom
monolayer $\mathbf{K}_{2}$, $\mathbf{K}^{\prime}_{2}$ in R-stacked and
H-stacked heterobilayers. (d) A closer look at the MBZ of the heterobilayer
and high-symmetry points.

There are two types of stacking for TMDC heterobilayers: R-stacking and
H-stacking. R-stacked TMDC heterobilayers are stacking of two triangular
lattices with near $0^{\circ}$ rotation and H-stacked TMDC heterobilayers are
stacking with near $60^{\circ}$ rotationjiang2021interlayer . Small twist-
angle variants straying from $0^{\circ}$ and $60^{\circ}$ rotations can still
be considered as R-stacked and H-stacked heterobilayers. Schematic plots of
the moiré superlattice and MBZ of the R-stacked heterobilayer are illustrated
in Fig. 1 (a), (b). The lattice constant of the moiré superlattice
$a_{\text{M}}$ as a function of twist angle $\theta$ and the lattice mismatch
$\delta=|a^{\prime}_{0}-a_{0}|/a_{0}$, with $a_{0}$, $a^{\prime}_{0}$ the
lattice constants for the atomistic lattices, is given bytran2020moire ;
zhang2023controlling

$\displaystyle
a_{\text{M}}=\frac{(1+\delta)a_{0}}{\sqrt{2(1+\delta)(1-\cos\theta)+\delta^{2}}}.$
(1)

The MBZ is built by the difference between reciprocal vectors in the Brillouin
zone of the top monolayer (inside blue dash-dot hexagonal) and the Brillouin
zone of the bottom monolayer (inside orange dot hexagonal) as illustrated in
Fig. 1 (b). The moiré superlattice can be considered a triangular lattice with
the unit cell containing local geometry. The primitive vectors of the moiré
superlattice can be defined as
$\bm{a}_{1}=a_{\text{M}}[(\sqrt{3}/2)\mathbf{e}_{x}+(1/2)\mathbf{e}_{y}]$,
$\bm{a}_{2}=a_{\text{M}}[(-\sqrt{3}/2)\mathbf{e}_{x}+(1/2)\mathbf{e}_{y}]$.
For the moiré reciprocal lattice, the reciprocal primitive vectors can be
defined by $\bm{g}_{i}\cdot\bm{a}_{j}=2\pi\delta_{ij}$ for $i,j=1,2$ and found
to be
$\bm{g}_{1}=\sqrt{3}k_{\text{M}}[(1/2)\mathbf{e}_{x`}+(\sqrt{3}/2)\mathbf{e}_{y}]$,
$\bm{g}_{2}=\sqrt{3}k_{\text{M}}[(-1/2)\mathbf{e}_{x}+(\sqrt{3}/2)\mathbf{e}_{y}]$,
with $k_{\text{M}}=4\pi/(3a_{\text{M}})$. A set of reciprocal primitive
vectors can be defined as
$\bm{g}_{j}=\sqrt{3}k_{\text{M}}\left[\mathbf{e}_{x}\cos({j\pi}/{3})+\mathbf{e}_{y}\sin({j\pi}/{3})\right]$,
with $j=1,2,\cdots,6$. The vector connects between $\kappa$ and $\gamma$ is
given by $\bm{\kappa}_{1}=\left(2\bm{g}_{1}-\bm{g}_{2}\right)/3$ and the
vector connects between $\kappa^{\prime}$ and $\gamma$ is given by
$\bm{\kappa}_{2}=\left(\bm{g}_{1}-2\bm{g}_{2}\right)/3$. These two vectors
$\bm{\kappa}_{1}$ and $\bm{\kappa}_{2}$ can be assigned as reciprocal basis
vectors. A set of reciprocal basis vectors can be defined as
$\bm{\kappa}_{j}=k_{\text{M}}\left[\mathbf{e}_{x}\cos(\pi/6-{j\pi}/{3})+\mathbf{e}_{y}\sin(\pi/6-{j\pi}/{3})\right]$,
with $j=1,2,\cdots,6$. Part of the reciprocal primitive and basis vectors are
shown in Fig. 1 (d).

Energy-level diagrams of band edges on different layers with the same spin
index on R-stacked and H-stacked heterobilayers are illustrated in Fig. 1
(c)rivera2018interlayer . Spin-orbit coupling (SOC) induced energy level
splittings occur in both conduction bands and valence bands of both layers.
Due to the intrinsic spin-valley locking of TMDCs, the $\mathbf{K}$ and
$\mathbf{K}^{\prime}$ valleys with the same spin index have different band-
edge energies. As illustrated in Fig. 1 (b), interlayer tunneling only occurs
among energy levels with the same types of lines (solid line or dash line).
Optical transition is only allowed between valence and conduction bands at the
same valley ($\mathbf{K}_{1}$, $\mathbf{K}^{\prime}_{1}$, $\mathbf{K}_{2}$,
$\mathbf{K}^{\prime}_{2}$). Via this diagram and effective-mass-model
parameters for TMDC monolayers, band-edge energies for different intralayer or
interlayer excitons in R-stacked or H-stacked heterobilayers can be evaluated.

###  II.2 Charge continuum model

The many-body Hamiltonian for electrons and holes in moiré superlattices
including the electron-hole attraction is given by

$\displaystyle\hat{\mathcal{H}}$ $\displaystyle=$
$\displaystyle\sum_{\mathbf{k}\mathbf{k}^{\prime},ll^{\prime}}\Big{[}\hat{c}^{\dagger}_{l,\mathbf{k}}\tilde{H}_{\text{e},ll^{\prime}}(\mathbf{k},\mathbf{k}^{\prime})\hat{c}_{l^{\prime},\mathbf{k}^{\prime}}+\hat{d}^{\dagger}_{l,\mathbf{k}}\tilde{H}_{\text{h},ll^{\prime}}(\mathbf{k},\mathbf{k}^{\prime})\hat{d}_{l^{\prime},\mathbf{k}^{\prime}}\Big{]}$
(2)
$\displaystyle+\frac{1}{2S}\sum_{ll^{\prime},\mathbf{q}}\tilde{W}_{ll^{\prime}}(\mathbf{q})\Big{(}\hat{\varrho}_{\text{e},l,\mathbf{q}}\hat{\varrho}_{\text{e},l^{\prime},-\mathbf{q}}+\hat{\varrho}_{\text{h},l,\mathbf{q}}\hat{\varrho}_{\text{h},l^{\prime},-\mathbf{q}}\Big{)}$
$\displaystyle-\frac{1}{S}\sum_{ll^{\prime},\mathbf{q}}\tilde{W}_{ll^{\prime}}(\mathbf{q})\hat{\varrho}_{\text{e},l,\mathbf{q}}\hat{\varrho}_{\text{h},l^{\prime},-\mathbf{q}},$

where $\hat{c}^{\dagger}_{l,\mathbf{k}}$/$\hat{c}_{l,\mathbf{k}}$ is the
electron creation/annihilation operators at layer $l$ with wavevector
$\mathbf{k}$, $\hat{d}^{\dagger}_{l,\mathbf{k}}$/$\hat{d}_{l,\mathbf{k}}$ is
the hole creation/annihilation operators,
$\tilde{H}_{\text{c},ll^{\prime}}(\mathbf{k},\mathbf{k}^{\prime})$ is the
charge Hamiltonian with $\text{c}=\\{\text{e},\text{h}\\}$ the charge index
indicating electrons and holes,
$\hat{\varrho}_{\text{e},l,\mathbf{q}}=\sum_{\mathbf{k}}\hat{c}^{\dagger}_{l,\mathbf{k}+\mathbf{q}}\hat{c}_{l,\mathbf{k}}$
is the electron density operator,
$\hat{\varrho}_{\text{h},l,\mathbf{q}}=\sum_{\mathbf{k}}\hat{d}^{\dagger}_{l,\mathbf{k}+\mathbf{q}}\hat{d}_{l,\mathbf{k}}$
is the hole density operator, $S$ is the area of the lattice, and
$\tilde{W}_{ll^{\prime}}(\mathbf{q})=\int
e^{-\mathtt{i}\mathbf{q}\cdot\mathbf{r}}{W}_{ll^{\prime}}(\mathbf{r})\text{d}^{2}r$
is the screened Coulomb interaction. Note that we have ignored the spin and
valley degrees of freedom in this many-body Hamiltonian.

The charge Hamiltonian in real space is related to the charge Hamiltonian in
$k$-space through the Fourier transform
$\tilde{H}_{\text{c},ll^{\prime}}(\mathbf{k},\mathbf{k}^{\prime})=\int
e^{-\mathtt{i}(\mathbf{k}-\mathbf{k}^{\prime})\cdot\mathbf{r}}{H}_{\text{c},ll^{\prime}}(\mathbf{r})\text{d}^{2}r$
and is given bywu2018hubbard

$\displaystyle
H_{\text{c}}(\mathbf{r})=\begin{pmatrix}h_{\text{c},1}(\mathbf{r})&t_{\text{c}}(\mathbf{r})\\\
t^{*}_{\text{c}}(\mathbf{r})&h_{\text{c},2}(\mathbf{r})\end{pmatrix},$ (3)

where $h_{\text{c},l}(\mathbf{r})$ is the layer-charge Hamiltonian of the
$l$-th layer and $t_{\text{e}}(\mathbf{r})$, $t_{\text{h}}(\mathbf{r})$ is the
interlayer tunneling coupling. The layer-charge Hamiltonian is given by

$\displaystyle h_{\text{c},l}(\mathbf{r})$ $\displaystyle=$
$\displaystyle\epsilon_{\text{c},l}-z_{\text{c}}U_{l}(\mathbf{r})+\frac{|\mathbf{p}+z_{\text{c}}\bm{\kappa}_{l}|^{2}}{2m_{\text{c},l}},$
(4)

where $\mathbf{p}=-\mathtt{i}\bm{\nabla}$ is the momentum operator of the
charge particle, $z_{\text{e}}=-$ and $z_{\text{h}}=+$ are the charge parity,
$\epsilon_{\text{c},l}$ is the band-edge energy, $m_{\text{c},l}$ is the
charge mass at layer $l$, $U_{l}(\mathbf{r})$ is the moiré potential. The
moiré potential is given by the harmonic-function form

$\displaystyle
U_{l}(\mathbf{r})=2V_{l}\sum_{j=1,3,5}\cos\left(\bm{g}_{j}\cdot\mathbf{r}+\psi_{l}\right),$
(5)

and the interlayer tunneling couplings are given by

$\displaystyle
t_{\text{e}}(\mathbf{r})=w_{\text{e}}\left(1+e^{\mathtt{i}\bm{g}_{1}\cdot\mathbf{r}}+e^{\mathtt{i}\bm{g}_{2}\cdot\mathbf{r}}\right),$
(6) $\displaystyle
t_{\text{h}}(\mathbf{r})=w_{\text{h}}\left(1+e^{-\mathtt{i}\bm{g}_{1}\cdot\mathbf{r}}+e^{-\mathtt{i}\bm{g}_{2}\cdot\mathbf{r}}\right),$
(7)

where $V_{l}$ is the moiré potential depth at layer $l$, $\psi_{l}$ is the
phase angle, and $w_{\text{e}}$/$w_{\text{h}}$ are the electron/hole
interlayer-tunneling coupling strength. In this work, the phase angle is
always chosen to be $\psi_{l}=-(-1)^{l}\pi/2$, such that the moiré potential
can be rewritten as
$U_{l}(\mathbf{r})=2V_{l}(-1)^{l}\sum_{j=1,3,5}\sin\left(\bm{g}_{j}\cdot\mathbf{r}\right)$.

The screened Coulomb potential is described by a modified Rytova-Keldysh
potentialrytova ; keldysh  which is written as

$\displaystyle\tilde{W}_{ll^{\prime}}(\mathbf{q})=\frac{2\pi}{\epsilon_{ll^{\prime}}({q})q},$
(8)

where $\epsilon_{ll^{\prime}}({q})$ is the bilayer dielectric function
determined by intrinsic properties of heterobilayers. While there are many
formulations for the bilayer dielectric functiondanovich2018localized ;
ruiz2020theory ; thygesen2017calculating ; cavalcante2018electrostatics ;
van2018interlayer ; florez2020effects , here we use a simple form derived in
Appendix A. The intralayer and interlayer dielectric functions are given by

$\displaystyle\epsilon_{11}(q)$ $\displaystyle=$
$\displaystyle\frac{\kappa_{0}\epsilon_{12}(q)}{\left(\frac{\kappa_{1}+\kappa_{0}}{2}+r_{2}q\right)e^{q{d}}-\left(\frac{\kappa_{1}-\kappa_{0}}{2}+r_{2}q\right)e^{-q{d}}},$
(9) $\displaystyle\epsilon_{22}(q)$ $\displaystyle=$
$\displaystyle\frac{\kappa_{0}\epsilon_{12}(q)}{\left(\frac{\kappa_{1}+\kappa_{0}}{2}+r_{1}q\right)e^{q{d}}-\left(\frac{\kappa_{1}-\kappa_{0}}{2}+r_{1}q\right)e^{-q{d}}},\hskip
4.30554pt$ (10) $\displaystyle\epsilon_{12}(q)$ $\displaystyle=$
$\displaystyle\left(\frac{\kappa_{1}+\kappa_{0}}{2}+r_{1}q\right)\left(\frac{\kappa_{1}+\kappa_{0}}{2}+r_{2}q\right)\frac{e^{q{d}}}{\kappa_{0}}$
$\displaystyle-\left(\frac{\kappa_{1}-\kappa_{0}}{2}+r_{1}q\right)\left(\frac{\kappa_{1}-\kappa_{0}}{2}+r_{2}q\right)\frac{e^{-q{d}}}{\kappa_{0}},$

and $\epsilon_{21}(q)=\epsilon_{12}(q)$, where $\kappa_{1}$ is the dielectric
constant outside the bilayer, $\kappa_{0}$ is the dielectric constant inside
the bilayer, $r_{l}$ is the screening length on the $l$-th layer, and $d$ is
the interlayer distance.

###  II.3 Exciton continuum model

An exciton state is assumed to be written as
$\ket{\text{X}_{I}}=\text{X}^{\dagger}_{I}\ket{0}$, where

$\displaystyle\text{X}^{\dagger}_{I}$ $\displaystyle=$
$\displaystyle\sum_{l_{\text{e}},l_{\text{h}}}\sum_{\mathbf{k}_{\text{e}},\mathbf{k}_{\text{h}}}\tilde{\Psi}_{l_{\text{e}}l_{\text{h}},I}(\mathbf{k}_{\text{e}},\mathbf{k}_{\text{h}})\hat{c}^{\dagger}_{l_{\text{e}},\mathbf{k}_{\text{e}}}\hat{d}^{\dagger}_{l_{\text{h}},\mathbf{k}_{\text{h}}}$
(12)

is the exciton creation operator, $\ket{0}$ is the ground-state ket, and
$\tilde{\Psi}_{I}(\mathbf{k}_{\text{e}},\mathbf{k}_{\text{h}})$ is the exciton
wavefunction. By variation method,
$\delta[\braket{\text{X}_{I}}{\hat{\mathcal{H}}}{\text{X}_{I}}-\lambda(\braket{\text{X}_{I}}{\text{X}_{I}}-1)]=0$,
an eigenvalue equation can be derived as

$\displaystyle\sum_{\mathbf{k}^{\prime}_{\text{e}},\mathbf{k}^{\prime}_{\text{h}}}\tilde{\mathcal{H}}_{\text{X}}(\mathbf{k}_{\text{e}},\mathbf{k}_{\text{h}};\mathbf{k}^{\prime}_{\text{e}},\mathbf{k}^{\prime}_{\text{h}})\tilde{\Psi}_{I}(\mathbf{k}^{\prime}_{\text{e}},\mathbf{k}^{\prime}_{\text{h}})=\varepsilon_{\text{X},I}\tilde{\Psi}_{I}(\mathbf{k}_{\text{e}},\mathbf{k}_{\text{h}}),$

where
$\tilde{\mathcal{H}}_{\text{X}}(\mathbf{k}_{\text{e}},\mathbf{k}_{\text{h}};\mathbf{k}^{\prime}_{\text{e}},\mathbf{k}^{\prime}_{\text{h}})=\bra{0}\hat{d}_{l_{\text{h}},\mathbf{k}_{\text{h}}}\hat{c}_{l_{\text{e}},\mathbf{k}_{\text{e}}}\hat{\mathcal{H}}\hat{c}^{\dagger}_{l^{\prime}_{\text{e}},\mathbf{k}^{\prime}_{\text{e}}}\hat{d}^{\dagger}_{l^{\prime}_{\text{h}},\mathbf{k}^{\prime}_{\text{h}}}\ket{0}$
is the momentum-space exciton Hamiltonian, and $\varepsilon_{\text{X},I}$ is
the exciton eigenenergy. The continuum model to describe excitons in 2D moiré
superlattices can be derived from the Fourier transform of the momentum-space
exciton Hamiltonian

$\displaystyle\mathcal{H}_{\text{X}}(\mathbf{r}_{\text{e}},\mathbf{r}_{\text{h}})$
$\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{\Delta\mathbf{k}_{\text{e}},\Delta\mathbf{k}_{\text{h}}}\exp\left({\mathtt{i}\Delta\mathbf{k}_{\text{e}}\cdot\mathbf{r}_{\text{e}}}\right)\exp\left({\mathtt{i}\Delta\mathbf{k}_{\text{h}}\cdot\mathbf{r}_{\text{h}}}\right)$
(14)
$\displaystyle\times\tilde{\mathcal{H}}_{\text{X}}(\mathbf{k}_{\text{e}}+\Delta\mathbf{k}_{\text{e}},\mathbf{k}_{\text{h}}+\Delta\mathbf{k}_{\text{h}};\mathbf{k}_{\text{e}},\mathbf{k}_{\text{h}}),$

and it is written as

$\displaystyle\mathcal{H}_{\text{X}}(\mathbf{r}_{\text{e}},\mathbf{r}_{\text{h}})=$
$\displaystyle\begin{pmatrix}\mathcal{H}_{11}(\mathbf{r}_{\text{e}},\mathbf{r}_{\text{h}})&0&\mathcal{T}_{\text{h}}(\mathbf{r}_{\text{e}},\mathbf{r}_{\text{h}})&\mathcal{T}_{\text{e}}(\mathbf{r}_{\text{e}},\mathbf{r}_{\text{h}})\\\
0&\mathcal{H}_{22}(\mathbf{r}_{\text{e}},\mathbf{r}_{\text{h}})&\mathcal{T}^{*}_{\text{e}}(\mathbf{r}_{\text{e}},\mathbf{r}_{\text{h}})&\mathcal{T}^{*}_{\text{h}}(\mathbf{r}_{\text{e}},\mathbf{r}_{\text{h}})\\\
\mathcal{T}^{*}_{\text{h}}(\mathbf{r}_{\text{e}},\mathbf{r}_{\text{h}})&\mathcal{T}_{\text{e}}(\mathbf{r}_{\text{e}},\mathbf{r}_{\text{h}})&\mathcal{H}_{12}(\mathbf{r}_{\text{e}},\mathbf{r}_{\text{h}})&0\\\
\mathcal{T}^{*}_{\text{e}}(\mathbf{r}_{\text{e}},\mathbf{r}_{\text{h}})&\mathcal{T}_{\text{h}}(\mathbf{r}_{\text{e}},\mathbf{r}_{\text{h}})&0&\mathcal{H}_{21}(\mathbf{r}_{\text{e}},\mathbf{r}_{\text{h}})\end{pmatrix},$

where
$\mathcal{T}_{\text{e}}(\mathbf{r}_{\text{e}},\mathbf{r}_{\text{h}})=t_{\text{e}}(\mathbf{r}_{\text{e}})$,
$\mathcal{T}_{\text{h}}(\mathbf{r}_{\text{e}},\mathbf{r}_{\text{h}})=t_{\text{h}}(\mathbf{r}_{\text{e}})$
are interlayer charge-transfer couplings,

$\displaystyle\mathcal{H}_{l_{\text{e}}l_{\text{h}}}(\mathbf{r}_{\text{e}},\mathbf{r}_{\text{h}})$
$\displaystyle=$
$\displaystyle\Delta_{l_{\text{e}}l_{\text{h}}}+\frac{|\mathbf{p}_{\text{e}}-\bm{\kappa}_{l_{\text{e}}}|^{2}}{2m_{\text{e},l_{\text{e}}}}+\frac{|\mathbf{p}_{\text{h}}+\bm{\kappa}_{l_{\text{h}}}|^{2}}{2m_{\text{h},l_{\text{h}}}}$
(16)
$\displaystyle+U_{l_{\text{e}}}(\mathbf{r}_{\text{e}})-U_{l_{\text{h}}}(\mathbf{r}_{\text{h}})-W_{l_{\text{e}}l_{\text{h}}}(\mathbf{r}_{\text{eh}}),$

is the layer-exciton Hamiltonian,
$\Delta_{l_{\text{e}}l_{\text{h}}}=\epsilon_{\text{e},l_{\text{e}}}+\epsilon_{\text{h},l_{\text{h}}}$
is the band-gap energy, and $W_{l_{\text{e}}l_{\text{h}}}(\mathbf{r})=\int
e^{\mathtt{i}\mathbf{k}\cdot\mathbf{r}}\tilde{W}(\mathbf{k})\text{d}^{2}k/(2\pi)^{2}$
the screened Coulomb potential. The exciton continuum model can be
reformulated by the coordinate transformation from electron position and hole
position coordinates to exciton COM and internal coordinates,
$\mathbf{R}=\gamma_{\text{e},l_{\text{e}}}\mathbf{r}_{\text{e}}+\gamma_{\text{h},l_{\text{h}}}\mathbf{r}_{\text{h}}$,
$\mathbf{r}=\mathbf{r}_{\text{e}}-\mathbf{r}_{\text{h}}$, with
$\gamma_{\text{e},l_{\text{e}}}=m_{\text{e},l_{\text{e}}}/(m_{\text{e},l_{\text{e}}}+m_{\text{h},l_{\text{h}}})$,
$\gamma_{\text{h},l_{\text{h}}}=m_{\text{h},l_{\text{h}}}/(m_{\text{e},l_{\text{e}}}+m_{\text{h},l_{\text{h}}})$,
which leads to the following transformation for momentums
$\mathbf{p}_{\text{e}}=\gamma_{\text{e},l_{\text{e}}}\mathbf{P}+\mathbf{p}$,
$\mathbf{p}_{\text{h}}=\gamma_{\text{h},l_{\text{h}}}\mathbf{P}-\mathbf{p}$,
with $\mathbf{P}$ the COM momentum and $\mathbf{p}$ the internal momentum. The
same method has been used in Ref. hannachi2023moire, . The reciprocal basis
vectors are also applied by the coordinate transformation

$\displaystyle\bm{\mathcal{K}}_{l_{\text{e}}l_{\text{h}}}=\bm{\kappa}_{l_{\text{e}}}-\bm{\kappa}_{l_{\text{h}}},\hskip
8.61108pt\bm{\kappa}_{l_{\text{e}}l_{\text{h}}}=\gamma_{\text{h},l_{\text{h}}}\bm{\kappa}_{l_{\text{e}}}+\gamma_{\text{e},l_{\text{e}}}\bm{\kappa}_{l_{\text{h}}},$
(17)

with $\bm{\mathcal{K}}_{l_{\text{e}}l_{\text{h}}}$ the COM reciprocal basis
vector and $\bm{\kappa}_{l_{\text{e}}l_{\text{h}}}$ the internal reciprocal
basis vector. The exciton Hamiltonian can be rewritten as

$\displaystyle\mathcal{H}_{\text{X}}(\mathbf{R},\mathbf{r})$
$\displaystyle=\begin{pmatrix}\mathcal{H}_{11}(\mathbf{R},\mathbf{r})&0&\mathcal{T}_{\text{h}}(\mathbf{R},\mathbf{r})&\mathcal{T}_{\text{e}}(\mathbf{R},\mathbf{r})\\\
0&\mathcal{H}_{22}(\mathbf{R},\mathbf{r})&\mathcal{T}^{*}_{\text{e}}(\mathbf{R},\mathbf{r})&\mathcal{T}^{*}_{\text{h}}(\mathbf{R},\mathbf{r})\\\
\mathcal{T}^{*}_{\text{h}}(\mathbf{R},\mathbf{r})&\mathcal{T}_{\text{e}}(\mathbf{R},\mathbf{r})&\mathcal{H}_{12}(\mathbf{R},\mathbf{r})&0\\\
\mathcal{T}^{*}_{\text{e}}(\mathbf{R},\mathbf{r})&\mathcal{T}_{\text{h}}(\mathbf{R},\mathbf{r})&0&\mathcal{H}_{21}(\mathbf{R},\mathbf{r})\end{pmatrix},$

where

$\displaystyle\mathcal{T}_{\text{e}}(\mathbf{R},\mathbf{r})$ $\displaystyle=$
$\displaystyle
w_{\text{e}}\big{[}1+\exp\left({\mathtt{i}\bm{g}_{1}\cdot\left(\mathbf{R}+\gamma_{\text{h},l_{\text{h}}}\mathbf{r}\right)}\right)$
(19)
$\displaystyle+\exp\left({\mathtt{i}\bm{g}_{2}\cdot\left(\mathbf{R}+\gamma_{\text{h},l_{\text{h}}}\mathbf{r}\right)}\right)\big{]},$
$\displaystyle\mathcal{T}_{\text{h}}(\mathbf{R},\mathbf{r})$ $\displaystyle=$
$\displaystyle
w_{\text{h}}\big{[}1+\exp\left(-{\mathtt{i}\bm{g}_{1}\cdot\left(\mathbf{R}-\gamma_{\text{e},l_{\text{e}}}\mathbf{r}\right)}\right)$
(20)
$\displaystyle+\exp\left(-{\mathtt{i}\bm{g}_{2}\cdot\left(\mathbf{R}-\gamma_{\text{e},l_{\text{e}}}\mathbf{r}\right)}\right)\big{]},$
$\displaystyle\mathcal{H}_{l_{\text{e}}l_{\text{h}}}(\mathbf{R},\mathbf{r})$
$\displaystyle=$
$\displaystyle\Delta_{l_{\text{e}}l_{\text{h}}}+\frac{|\mathbf{P}-\bm{\mathcal{K}}_{l_{\text{e}}l_{\text{h}}}|^{2}}{2m_{\text{X},l_{\text{e}}l_{\text{h}}}}+\frac{|\mathbf{p}-\bm{\kappa}_{l_{\text{e}}l_{\text{h}}}|^{2}}{2\mu_{\text{X},l_{\text{e}}l_{\text{h}}}}$
(21)
$\displaystyle+\mathcal{V}_{l_{\text{e}}l_{\text{h}}}(\mathbf{R},\mathbf{r})-W_{l_{\text{e}}l_{\text{h}}}({r}),$

with $m_{\text{X},l_{\text{e}}l_{\text{h}}}$ the exciton mass,
$\mu_{\text{X},l_{\text{e}}l_{\text{h}}}$ the reduced mass, and
$\mathcal{V}_{l_{\text{e}}l_{\text{h}}}(\mathbf{R},\mathbf{r})$ the exciton
moiré potential. The exciton mass and the reduced mass are given by

$\displaystyle
m_{\text{X},l_{\text{e}}l_{\text{h}}}=m_{\text{e},l_{\text{e}}}+m_{\text{h},l_{\text{h}}},\hskip
8.61108pt\mu_{\text{X},l_{\text{e}}l_{\text{h}}}=\frac{m_{\text{e},l_{\text{e}}}m_{\text{h},l_{\text{h}}}}{m_{\text{X},l_{\text{e}}l_{\text{h}}}}.$
(22)

The exciton moiré potential is a combination of the electron moiré potential
and the hole moiré potential. The exciton moiré potential is given by

$\displaystyle\mathcal{V}_{l_{\text{e}}l_{\text{h}}}(\mathbf{R},\mathbf{r})$
$\displaystyle=$ $\displaystyle
U_{l_{\text{e}}}(\mathbf{R}+\gamma_{\text{h},l_{\text{h}}}\mathbf{r})-U_{l_{\text{h}}}(\mathbf{R}-\gamma_{\text{e},l_{\text{e}}}\mathbf{r}).$
(23)

Note that the exciton band minimum is centered at
$\bm{\mathcal{K}}_{l_{\text{e}}l_{\text{h}}}=\bm{\kappa}_{l_{\text{e}}}-\bm{\kappa}_{l_{\text{h}}}$.
For intralayer excitons $l_{\text{e}}=l_{\text{h}}=l$ thus
$\bm{\mathcal{K}}_{l_{\text{e}}l_{\text{h}}}=0$, the exciton band minimum
locates at $\gamma$ point in MBZ. For interlayer exciton, $l_{\text{e}}=1$,
$l_{\text{h}}=2$ or $l_{\text{e}}=2$, $l_{\text{h}}=l$ thus
$\bm{\mathcal{K}}_{l_{\text{e}}l_{\text{h}}}=\pm(\bm{\kappa}_{2}-\bm{\kappa}_{1})$,
the exciton band minimum locates at one of $\kappa$, $\kappa^{\prime}$ points
in the MBZ. Based on Fig. 1 (b), the momentum shift of interlayer excitons
comes from the momentum difference between $\mathbf{K}_{1}$
($\mathbf{K}^{\prime}_{1}$) and $\mathbf{K}_{2}$ ($\mathbf{K}^{\prime}_{2}$)
points on the Brillouin zones of the top and bottom layers. Therefore, the
intralayer exciton is assigned as a direct exciton and the interlayer exciton
is assigned as an indirect exciton. Optical transitions only occur for
intralayer excitons located at $\mathbf{K}=\mathbf{0}$, which is $\gamma$
point on the MBZ of excitons.

It is worth a mention that, by Taylor expansion of the internal coordinate,
the exciton moiré potential can be approximated as

$\displaystyle\mathcal{V}_{l_{\text{e}}l_{\text{h}}}(\mathbf{R},\mathbf{r})$
$\displaystyle\simeq$ $\displaystyle
U_{l_{\text{e}}}(\mathbf{R})+\gamma_{\text{h},l_{\text{h}}}\mathbf{r}\cdot\bm{\nabla}_{\mathbf{R}}U_{l_{\text{e}}}(\mathbf{R})$
(24) $\displaystyle-
U_{l_{\text{h}}}(\mathbf{R})+\gamma_{\text{e},l_{\text{e}}}\mathbf{r}\cdot\bm{\nabla}_{\mathbf{R}}U_{l_{\text{h}}}(\mathbf{R}).$

For intralayer excitons, $l_{\text{e}}=l_{\text{h}}=l$ and
$\psi_{l}=-(-1)^{l}\pi/2$ are used, and the exciton moiré potential becomes

$\displaystyle\mathcal{V}_{ll}(\mathbf{R},\mathbf{r})$ $\displaystyle\simeq$
$\displaystyle
2{V}_{l}(-1)^{l}\sum_{j=1,3,5}\left(\bm{g}_{j}\cdot\mathbf{r}\right)\cos(\bm{g}_{j}\cdot\mathbf{R}).$
(25)

On the other hand, for interlayer excitons, $l_{\text{e}}=2$ and
$l_{\text{h}}=1$, the exciton moiré potential becomes

$\displaystyle\mathcal{V}_{21}(\mathbf{R},\mathbf{r})$ $\displaystyle\simeq$
$\displaystyle 2({V}_{1}+V_{2})\sum_{j=1,3,5}\sin(\bm{g}_{j}\cdot\mathbf{R}).$
(26)

Therefore, exciton moiré potentials for intralayer excitons and interlayer
excitons are different. Additionally, the exciton moiré potential jointly
couples the COM motion and the internal motion. The COM wavefunction and the
internal wavefunction of an intralayer moiré exciton are entangled for any
moiré potential depths.

###  II.4 Exciton wavefunction

The exciton wavefunction and exciton band structure can be solved by the
eigenvalue problem

$\displaystyle\mathcal{H}_{\text{X}}(\mathbf{R},\mathbf{r})\Psi_{\text{X},I\mathbf{K}}(\mathbf{R},\mathbf{r})=\varepsilon_{\text{X},I\mathbf{K}}\Psi_{\text{X},I\mathbf{K}}(\mathbf{R},\mathbf{r}),$
(27)

where $\Psi_{\text{X},I\mathbf{K}}(\mathbf{R},\mathbf{r})$ is the exciton
wavefunction and $\varepsilon_{\text{X},I\mathbf{K}}$ is the exciton
eigenenergy. A variational exciton wavefunction method is used to solve the
exciton Hamiltonian. To include the entanglement between the COM motion and
the internal motion, the variational exciton wavefunction should contain a
direct product of the COM-coordinate function and the internal-coordinate
function. Such variational exciton wavefunction is written as

$\displaystyle\Psi_{\text{X},I\mathbf{K}}(\mathbf{R},\mathbf{r})$
$\displaystyle=$
$\displaystyle\frac{e^{\mathtt{i}\mathbf{K}\cdot\mathbf{R}}}{2\pi}\sum_{l_{\text{e}}l_{\text{h}},a,\mathbf{G}}C_{(l_{\text{e}}l_{\text{h}},{a},\mathbf{G}),I\mathbf{K}}\;\chi_{\mathbf{G}}(\mathbf{R})$
(28)
$\displaystyle\times\exp\left({\mathtt{i}\bm{\kappa}_{l_{\text{e}}l_{\text{h}}}\cdot\mathbf{r}}\right)\phi_{l_{\text{e}}l_{\text{h}},a}(\mathbf{r})\mathbb{E}_{l_{\text{e}}l_{\text{h}}},$

where $C_{(l_{\text{e}}l_{\text{h}},{a},\mathbf{G}),I\mathbf{K}}$ is the
variational wavefunction coefficient, $\chi_{\mathbf{G}}(\mathbf{R})$ is the
COM-coordinate basis function, $\phi_{l_{\text{e}}l_{\text{h}},a}(\mathbf{r})$
is the internal-coordinate basis function, and
$\mathbb{E}_{l_{\text{e}}l_{\text{h}}}$ is the exciton Hamiltonian-matrix
basis vector with

$\displaystyle\mathbb{E}_{11}=\begin{pmatrix}1\\\ 0\\\ 0\\\
0\end{pmatrix},\hskip 4.30554pt\mathbb{E}_{22}=\begin{pmatrix}0\\\ 1\\\ 0\\\
0\end{pmatrix},\hskip 4.30554pt\mathbb{E}_{12}=\begin{pmatrix}0\\\ 0\\\ 1\\\
0\end{pmatrix},\hskip 4.30554pt\mathbb{E}_{21}=\begin{pmatrix}0\\\ 0\\\ 0\\\
1\end{pmatrix}.$

The exciton COM degree of freedom can be solved as the problem of a particle
in a periodic potential. Therefore, the COM-coordinate basis function is given
by the plane-wave function

$\displaystyle\chi_{\mathbf{G}}(\mathbf{R})=e^{-\mathtt{i}\mathbf{G}\cdot\mathbf{R}}$
(30)

where $\mathbf{G}$ is reciprocal lattice vector. The exciton internal degree
of freedom can be treated as the eigenspectrum problem of an isolated exciton.
It is known that STOs can be used as the basis function to solve exciton
internal wavefunctions accurately in two dimensionszhang2019two ;
wu2019exciton ; mypaper0 ; mypaper1 . The internal-coordinate basis function
is given by STOs, which have the form

$\displaystyle\phi_{l_{\text{e}}l_{\text{h}},a}(\mathbf{r})=\left({e^{\mathtt{i}L_{a}\varphi}}/{\sqrt{2\pi}}\right)r^{N_{a}-1}\exp\left({-\mathcal{Z}_{l_{\text{e}}l_{\text{h}},a}r}\right)$
(31)

with $N_{a}$ the shell number, $L_{a}$ the angular momentum, and
$\mathcal{Z}_{l_{\text{e}}l_{\text{h}},a}$ the shielding constant. Note that
the phase factor
$\exp\left(\mathtt{i}\bm{\kappa}_{l_{\text{e}}l_{\text{h}}}\cdot\mathbf{r}\right)$
has been added in Eq. (28). This phase factor is utilized to shift the
internal momentum to center the origin point,

$\displaystyle\exp\left(-{\mathtt{i}\bm{\kappa}_{l_{\text{e}}l_{\text{h}}}\cdot\mathbf{r}}\right)\left(\mathbf{p}-\bm{\kappa}_{l_{\text{e}}l_{\text{h}}}\right)\exp\left({\mathtt{i}\bm{\kappa}_{l_{\text{e}}l_{\text{h}}}\cdot\mathbf{r}}\right)=\mathbf{p},$
(32)

such that the matrix element can be written in a simpler form. Matrix elements
of the exciton Hamiltonian matrix and the overlap matrix span by the present
basis are shown in Appendix B.

By using the exciton wavefunction, the radius and angular momentum of the
hybridized moiré exciton can be found. The exciton radius operator is defined
by the matrix representation

$\displaystyle r_{\text{X}}(\mathbf{R},\mathbf{r})\equiv r\mathbb{I},$ (33)

with

$\displaystyle\mathbb{I}=\begin{pmatrix}{1}&0&0&0\\\ 0&{1}&0&0\\\ 0&0&{1}&0\\\
0&0&0&{1}\end{pmatrix}$ (34)

a four-dimensional identity matrix. The exciton radius can be solved from

$\displaystyle
a_{\text{X},I}=\int\Psi^{*}_{\text{X},I\mathbf{0}}(\mathbf{R},\mathbf{r})r_{\text{X}}(\mathbf{R},\mathbf{r})\Psi_{\text{X},I\mathbf{0}}(\mathbf{R},\mathbf{r})\text{d}^{2}R\text{d}^{2}r.\hskip
8.61108pt$ (35)

The exciton angular momentum operator is defined as
$\mathcal{L}_{\text{X}}(\mathbf{R},\mathbf{r})\equiv\mathcal{L}^{\prime}_{\text{X}}(\mathbf{R})+\mathcal{L}^{\prime\prime}_{\text{X}}(\mathbf{r})$,
where

$\displaystyle\mathcal{L}^{\prime}_{\text{X}}(\mathbf{R})=\left[\mathbf{R}\times(\mathbf{P}-\bm{\mathcal{K}}_{l_{\text{e}}l_{\text{h}}})\right]\mathbb{I}$
(36)

and

$\displaystyle\mathcal{L}^{\prime\prime}_{\text{X}}(\mathbf{r})=\left[\mathbf{r}\times(\mathbf{p}-\bm{\kappa}_{l_{\text{e}}l_{\text{h}}})\right]\mathbb{I}$
(37)

are the matrix representations of the exciton COM angular momentum operator
and the exciton internal angular momentum operator. Since the internal angular
momentum of optically-active excitons should be zero, the root-mean-square
exciton angular momentum is used to measure the contribution of non-zero-
angular-momentum Rydberg excitons on moiré exciton wavefunction. The root-
mean-square exciton angular momentum $\sqrt{\braket{\mathcal{L}^{\prime\prime
2}_{\text{X}}}}$ is given by the formula

$\displaystyle\braket{\mathcal{L}^{\prime\prime
2}_{\text{X}}}_{I}=\int\Psi^{*}_{\text{X},I\mathbf{0}}(\mathbf{R},\mathbf{r})\mathcal{L}^{\prime\prime
2}_{\text{X}}(\mathbf{r})\Psi_{\text{X},I\mathbf{0}}(\mathbf{R},\mathbf{r})\text{d}^{2}R\text{d}^{2}r.\hskip
8.61108pt$ (38)

Therefore, if the exciton wavefunction contains the STO
$\phi_{l_{\text{e}}l_{\text{h}},a}(\mathbf{r})$ with $L_{a}\neq{0}$, a non-
zero value of root-mean-square exciton angular momentum can be found.

###  II.5 Optical absorption

Based on Fermi’s golden rule, the optical absorption spectrum can be given by
the formula

$\displaystyle\Gamma(\omega)$ $\displaystyle=$
$\displaystyle\frac{2\pi}{\omega}\sum_{I,\mu}|j^{\mu}_{0I}|^{2}\delta(\omega-
E_{I}+E_{0}),$ (39)

where ${E}_{I}$ the $I$-th state eigenenergy,
$\bm{j}_{II^{\prime}}=\braket{I}{\hat{\bm{j}}}{I^{\prime}}$ is the transition
amplitude, with $\ket{I}$ the $I$-th state ket. The transition between a one-
exciton excited state and the ground state is known as the one-exciton
transition. By assuming the one-exciton excited states being excited states
$\ket{I}=\ket{\text{X}_{I\mathbf{K}}}$, the one-exciton transition amplitudes
is given by

$\displaystyle\bm{j}_{0,I\mathbf{K}}=\int\bm{\mathcal{J}}_{\text{X}}(\mathbf{r})\Psi_{\text{X},I\mathbf{K}}(\mathbf{R},\mathbf{r})\text{d}^{2}r\text{d}^{2}R,$
(40)

where
$\bm{\mathcal{J}}_{\text{X}}(\mathbf{r})=\begin{pmatrix}\bm{\mathcal{J}}_{11}(\mathbf{r})&\bm{\mathcal{J}}_{22}(\mathbf{r})&0&0\end{pmatrix}$,
with $\bm{\mathcal{J}}_{ll}(\mathbf{r})$ the momentum matrix element. The
momentum matrix element is given by

$\displaystyle{\mathcal{J}}^{\pm}_{ll}(\mathbf{r})$ $\displaystyle=$
$\displaystyle
e\big{[}{\mathcal{J}}^{x}_{ll}(\mathbf{r})\pm\mathtt{i}{\mathcal{J}}^{y}_{ll}(\mathbf{r})\big{]}/\sqrt{2}=\delta_{\tau,\pm}\delta(\mathbf{r})ev_{\text{F},l},$
(41)

where $v_{\text{F},l}$ is the Fermi velocity for each layer $l$. Note that
$\tau$ is the valley index and the matrix element is circular-polarization
selected. However, this selection can be ignored since we have not included
the valley degree of freedom in this model. By using the exciton wavefunction,
the one-exciton transition amplitude is given by

$\displaystyle j^{\pm}_{0,I\mathbf{K}}$ $\displaystyle=$ $\displaystyle
2\pi\sum_{l,a}\sum_{\mathbf{G}}ev_{\text{F},l}\delta_{\mathbf{K},\mathbf{G}}C_{(ll,a,\mathbf{G}),I\mathbf{G}}\phi_{ll,a}(\mathbf{0}).$
(42)

Since only the first Brillouin zone is considered in the calculation, the
selection for the COM wavevector is given by $\mathbf{K}=\mathbf{0}$ for
optical absorption. The ground-state energy can be assigned as $E_{0}=0$ and
the excited-state energy for every one-exciton excited state is given by
$E_{I}=\varepsilon_{\text{X},I\mathbf{0}}$. Combining all the information, the
optical absorption spectrum can be rewritten as

$\displaystyle\Gamma(\omega)$ $\displaystyle=$
$\displaystyle\sum_{I,\tau=\pm}\frac{|j^{\tau}_{0,I\mathbf{0}}|^{2}}{\varepsilon_{\text{X},I\mathbf{0}}}\frac{2\eta}{(\omega-\varepsilon_{\text{X},I\mathbf{0}})^{2}+\eta^{2}},$
(43)

where $\eta$ is a line-broadening factor. Eq. (43) is the working formula for
us to calculate the optical absorption spectrum of moiré heterobilayers. Note
that if
$C_{(ll,a,\mathbf{G}),I\mathbf{G}}=\delta_{l,1}C_{(11,a,\mathbf{G}),I\mathbf{G}}$
in Eq. (42), the formula in Eq. (43) can be reduced to the optical absorption
formula for intralayer exciton transitions in the layer $l=1$. The interlayer
exciton transition is only introduced through the wavefunction coefficient
$C_{(ll,a,\mathbf{G}),I\mathbf{G}}$ of the $I$-th exciton state, where each
exciton state is a linear combination of intralayer and interlayer exciton
wavefunctions.

##  III Applications

Materials |  $E_{\text{gap}}$(eV) |  $E^{\text{v}}_{\text{edge}}$(eV) |  $E^{\text{c}}_{\text{edge}}$(eV) |  $E^{\text{v}}_{\text{SO}}$(meV) |  $E^{\text{c}}_{\text{SO}}$(meV) |  $a$(Å) | $m_{\text{e}}/m_{0}$ | $m_{\text{h}}/m_{0}$ |  $r_{0}$(Å) |  $E_{\text{X}}$(meV) | ${10}^{3}v_{\text{F}}/c_{0}$
---|---|---|---|---|---|---|---|---|---|---|---
$\text{MoSe}_{2}$ | $1.874$ | $-5.750$ | $-3.876$ | $-184$ | $20$ | $3.288$ | $0.70$ | $0.70$ | $39$ | $232$ | $1.62$
$\text{WS}_{2}$ | $2.238$ | $-6.190$ | $-3.952$ | $-425$ | $-31$ | $3.154$ | $0.35$ | $0.35$ | $34$ | $177$ | $2.50$
$\text{WSe}_{2}$ | $1.890$ | $-5.490$ | $-3.600$ | $-462$ | $-37$ | $3.286$ | $0.40$ | $0.40$ | $45$ | $165$ | $2.15$
Table 1: Effective-mass-model parameters for TMDC monolayers. The dielectric
constants for the modified Rytova-Keldysh potential are given by
$\kappa_{1}=4.4$ and $\kappa_{0}=2.0$. $E_{\text{gap}}$ is the band-gap
energy, $E^{\text{v}}_{\text{edge}}$ is the valence band-edge energy,
$E^{\text{c}}_{\text{edge}}$ is the conduction band-edge energy,
$E^{\text{v}}_{\text{SO}}$ is the valence band spin-orbit splitting energy,
$E^{\text{c}}_{\text{SO}}$ is the conduction band spin-orbit splitting energy,
$a$ is the lattice constant of the TMDC monolayer, $m_{\text{e}}$ and
$m_{\text{h}}$ are effective electron and hole masses, $m_{0}$ is the bare
electron mass, $r_{0}$ is the screening length, $E_{\text{X}}$ is the exciton
binding energy for the $1s$-orbital exciton, $v_{\text{F}}$ is the Fermi
velocity, and $c_{0}$ is the light velocity. The exciton binding energy is
calculated by variationally solving the effective mass model with Rytova-
Keldysh potential as the screened Coulomb potential. The Fermi velocity is
estimated by $E_{\text{gap}}=(m_{\text{e}}+m_{\text{h}})v_{\text{F}}$. The
effective masses and screening lengths are acquired from Ref.
goryca2019revealing, . The lattice constants are acquired from Ref.
wilson1969transition, . The band-edge energies are acquired from Ref.
zhang2016systematic, . The spin-orbit splitting energies are acquired from
Ref. kormanyos2015k, .

In this section, the exciton continuum model is applied to the simulation of
optical absorption spectra of $\text{WSe}_{2}$/$\text{WS}_{2}$ heterobilayers
and $\text{MoSe}_{2}$/$\text{WS}_{2}$ heterobilayers. Twist-angle and
electric-field dependences of optical absorption spectra are calculated and
discussed. Parameters for the effective-mass band and the modified Rytova-
Keldysh potential of heterobilayers embedded in hexagonal boron nitride
substrates are given in Table 1. All other parameters, including moiré
potential depths and charge-transfer coupling strengths, are chosen to fit
simulated optical spectra with experimental observations.

###  III.1 Excitons in $\text{WSe}_{2}$/$\text{WS}_{2}$ heterobilayers

Figure 2: (a) Calculated absorption spectra of $\text{WSe}_{2}$ and $\text{WS}_{2}$ monolayers. (b) Calculated absorption spectra of R-stacked and H-stacked $\text{WSe}_{2}$/$\text{WS}_{2}$ heterobilayers, with parameters $w_{\text{e}}=w_{\text{h}}=20$ meV, $V_{\text{WSe}_{2}}=30$ meV, $V_{\text{WS}_{2}}=5$ meV, $\eta=5$ meV, and $d=7.0$ Å. Excitons |  $E_{exc}$ (eV) |  $E_{\text{X}}$ (meV) |  $|j_{\text{X}}|^{2}$ (a.u.) |  $a_{\text{X}}$ (Å) | $\sqrt{\braket{\mathcal{L}^{\prime\prime 2}_{\text{X}}}}$
---|---|---|---|---|---
$\text{A-}{1s}$ | $1.725$ | $165$ | $1.000$ | $14.0$ | $0$
$\text{A-}{2s}$ | $1.851$ | $39$ | $0.107$ | $66.8$ | $0$
$\text{A-}{2p}$ | $1.830$ | $50$ | $0$ | $43.0$ | $1$
I | $1.690$ | $195$ | $0.763$ | $19.9$ | $0.395$
II | $1.739$ | $144$ | $0.088$ | $16.5$ | $0.261$
III | $1.789$ | $98$ | $0.337$ | $39.0$ | $0.603$
$\alpha$ | $1.898$ | $119$ | $0.025$ | $38.0$ | $0.660$
$\beta$ | $1.947$ | $69$ | $0.017$ | $56.9$ | $1.001$
Table 2: Calculated properties of A excitons in $\text{WSe}_{2}$ monolayer
(top three rows) and moiré excitons in H-stacked
$\text{WSe}_{2}$/$\text{WS}_{2}$ heterobilayers (remaining rows). $E_{exc}$ is
the exciton transition energy, $E_{\text{X}}$ is the exciton binding energy,
$|j_{\text{X}}|^{2}$ is the transition amplitude, $a_{\text{X}}$ is the
exciton radius, and $\sqrt{\braket{\mathcal{L}^{\prime\prime 2}_{\text{X}}}}$
is the exciton root-mean-square angular momentum. The properties of intralayer
moiré excitons are calculated without considering the interlayer tunneling.
Figure 3: (a) Calculated twist-angle dependence and (b) electric-field
dependence of optical absorption spectra of H-stacked
$\text{WSe}_{2}$/$\text{WS}_{2}$ heterobilayers.

Calculated absorption spectra of $\text{WSe}_{2}$ monolayer, $\text{WS}_{2}$
monolayer, and $\text{WSe}_{2}$/$\text{WS}_{2}$ heterobilayers are shown in
Fig. 2 (a) and (b). The contributions from A and B excitons in
$\text{WSe}_{2}$ monolayer and A excitons in $\text{WS}_{2}$ monolayer have
been included in Fig. 2 (a). Signatures of $2s$, $3s$ Rydberg excitons in
$\text{WSe}_{2}$ monolayer are also shown around $1.85$ eV for A excitons and
around $2.25$ eV for B excitons as the satellite peaks of the dominated $1s$
exciton transition peak. By including moiré potentials and charge-transfer
couplings in the exciton calculation for $\text{WSe}_{2}$/$\text{WS}_{2}$
heterobilayers, signatures of intralayer and interlayer moiré excitons are
shown. The three prominent peaks around $1.7$-$1.8$ eV, labeled I, II, and III
peaks as in Ref. jin2019observation, , can be attributed to the intralayer
moiré excitons extending from A excitons of $\text{WSe}_{2}$. The inset
figures in Fig. 2 (b) show enlarged peaks from optical absorption by
interlayer excitons. The signature around $1.4$-$1.5$ eV in the spectrum of
the R-stacked heterobilayer can be attributed to the interlayer moiré excitons
with the electron locating at $\mathbf{K}$ valley of $\text{WS}_{2}$ monolayer
and the hole locating at $\mathbf{K}$ valley of $\text{WSe}_{2}$ monolayer.
The signature around $1.4$-$1.5$ eV in the spectrum of the H-stacked
heterobilayer can be attributed to the interlayer moiré excitons with the
electron locating at $\mathbf{K}^{\prime}$ valley of $\text{WS}_{2}$ monolayer
and the hole locating at $\mathbf{K}$ valley of $\text{WSe}_{2}$ monolayer.
The signature around $1.85$-$1.95$ eV in the spectrum of the H-stacked
heterobilayer can be attributed to the interlayer moiré excitons with the
electron locating at $\mathbf{K}$ valley of $\text{WS}_{2}$ monolayer and the
hole locating at $\mathbf{K}^{\prime}$ valley of $\text{WSe}_{2}$ monolayer.
Interlayer-exciton transition peaks are much less obvious in comparison with
intralayer-exciton transition peaks because of the large energy difference
between interlayer excitons and intralayer excitons. These signatures are
overall coincident with experimental observations in the
literaturejin2019identification ; tang2021tuning .

In Table 2, exciton transition energies, exciton binding energies, transition
amplitudes, exciton radii, and exciton root-mean-square angular momentums of
exciton transitions in $\text{WSe}_{2}$ monolayer and H-stacked
$\text{WSe}_{2}$/$\text{WS}_{2}$ heterobilayers are shown. Excitons $\alpha$,
$\beta$ are two interlayer-exciton transitions of H-stacked
$\text{WSe}_{2}$/$\text{WS}_{2}$ heterobilayers, as indicated in Fig. 2 (b).
Because the internal degree of freedom is included in the exciton calculation,
intralayer excitons and interlayer excitons can have different radii and
angular momentum. According to the variational exciton wavefunction in Eq.
(28), an exciton in heterobilayers can be seen as the hybridization of various
Rydberg excitons in monolayers. Therefore, hybridized excitons listed in Table
2 show mixture properties from different Rydberg excitons in $\text{WSe}_{2}$
monolayer.

In Fig. 3 (a), (b), calculated twist-angle and electric-field dependences of
optical absorption spectra in H-stacked $\text{WSe}_{2}$/$\text{WS}_{2}$
heterobilayers are shown. For the twist-angle dependence, the three
intralayer-exciton peaks I, II, and III show spectral shifts with increasing
twist angles. The peak I converges to the A-$1s$ exciton of $\text{WSe}_{2}$
monolayer, the peak III converges to the A-$2s$ exciton of $\text{WSe}_{2}$
monolayer, and the peak II diminishes. For the electric-field dependence, an
out-of-plane electric field is included in the model by modulating the band-
edge energy of $\text{WS}_{2}$ monolayer as

$\displaystyle\tilde{\epsilon}_{\text{e},\text{WS}_{2}}=\epsilon_{\text{e},\text{WS}_{2}}-\xi_{z}F_{z},\hskip
8.61108pt\tilde{\epsilon}_{\text{h},\text{WS}_{2}}=\epsilon_{\text{h},\text{WS}_{2}}-\xi_{z}F_{z},\hskip
8.61108pt$ (44)

where $\epsilon_{\text{e},\text{WS}_{2}}$ and
$\epsilon_{\text{h},\text{WS}_{2}}$ are band-edge energies for the conduction
band and the valence band, $\xi_{z}$ is the out-of-plane electric dipole
moment, and $F_{z}$ is the out-of-plane electric field strength. The electric
dipole moment is assigned as $\xi_{z}=0.4\;e\cdot$nm, where $e$ denotes the
elementary charge. Signatures of hybridization can be found in both the twist-
angle dependence and the electric-field dependence, particularly near the
A-exciton transition energy of $\text{WS}_{2}$ monolayer. As shown in Fig. 3
(a), signatures of avoiding crossing are found at $7^{\circ}$ twist angle and
around $2.1$ eV excitation energy. In Fig. 3 (b), signatures of avoiding
crossing are found around $-0.3\sim-0.5$ V nm-1 electric field and around
$2.1$ eV excitation energy. These signatures indicate the hybridization
between intralayer excitons and interlayer excitons through charge-transfer
couplings as the excitation energies of intralayer excitons and interlayer
excitons are close.

Note that, for real-world experiments, the moiré heterobilayer is embedded in
a dielectric substrate, which is normally composed of multiple layers of
hexagonal boron nitrides (hBNs) or silicon dioxide (SiO2), and the applied
external electric field can be screened by the substrate. Therefore, the out-
of-plane electric field mentioned above is the out-of-plane displacement field
and the simulated electric-field dependence might be different from the
experimental observation due to different dielectric screening for different
experimental conditions.

###  III.2 Excitons in $\text{MoSe}_{2}$/$\text{WS}_{2}$ heterobilayers

Figure 4: (a) Calculated absorption spectra of $\text{MoSe}_{2}$ and
$\text{WS}_{2}$ monolayers. (b) Calculated absorption spectra of R-stacked and
H-stacked $\text{MoSe}_{2}$/$\text{WS}_{2}$ heterobilayers, with parameters
$w_{\text{e}}=w_{\text{h}}=10$ meV, $V_{\text{MoSe}_{2}}=35$ meV,
$V_{\text{WS}_{2}}=5$ meV, $\eta=5$ meV, and $d=7.0$ Å. Figure 5: Calculated
twist-angle-dependent absorption spectra of (a) R-stacked
$\text{MoSe}_{2}$/$\text{WS}_{2}$ heterobilayers and (b) H-stacked
$\text{MoSe}_{2}$/$\text{WS}_{2}$ heterobilayers. Figure 6: Calculated
electric-field-dependent absorption spectra of (a) R-stacked
$\text{MoSe}_{2}$/$\text{WS}_{2}$ heterobilayers and (b) H-stacked
$\text{MoSe}_{2}$/$\text{WS}_{2}$ heterobilayers.

Calculated absorption spectra of $\text{MoSe}_{2}$ monolayer, $\text{WS}_{2}$
monolayer, R-stacked and H-stacked $\text{MoSe}_{2}$/$\text{WS}_{2}$
heterobilayers are shown in Fig. 4 (a) and (b). Again, optical absorption
signatures of A excitons, B excitons, and $2s$, $3s$ Rydberg excitons of
$\text{MoSe}_{2}$ monolayer and $\text{WS}_{2}$ monolayer are shown in Fig. 4
(a). The interlayer-exciton transition energy of the R-stacked heterobilayer
is close to the A-exciton transition energy of $\text{MoSe}_{2}$ monolayer and
the interlayer-exciton transition energy of the H-stacked heterobilayer is
close to the B-exciton transition energy of $\text{MoSe}_{2}$ monolayer. The
hybridization between the intralayer excitons and the interlayer excitons
induces small changes in the shape of peaks around the A-exciton transition
energy of $\text{MoSe}_{2}$ monolayer ($1.6$-$1.7$ eV) for the R-stacked
heterobilayer and around the B-exciton transition energy of $\text{MoSe}_{2}$
monolayer ($1.8$-$1.9$ eV) for the H-stacked heterobilayer.

In Fig. 5, calculated twist-angle-dependent absorption spectra of R-stacked
$\text{MoSe}_{2}$/$\text{WS}_{2}$ heterobilayers and H-stacked
$\text{MoSe}_{2}$/$\text{WS}_{2}$ heterobilayers are shown, and in Fig. 6,
calculated electric-field-dependent absorption spectra of R-stacked
$\text{MoSe}_{2}$/$\text{WS}_{2}$ heterobilayers and H-stacked
$\text{MoSe}_{2}$/$\text{WS}_{2}$ heterobilayers are shown. The out-of-plane
electric field is included in the model by the same band-edge modulation of
$\text{WS}_{2}$ monolayer in Eq. (44). Signatures of avoiding crossing due to
exciton hybridization can be found in both simulated spectra. Experimental
twist-angle-dependent absorption spectra of R-stacked and H-stacked
$\text{MoSe}_{2}$/$\text{WS}_{2}$ heterobilayers can be found in Ref.
alexeev2019resonantly, . Experimental electric-field-dependent absorption
spectra of R-stacked $\text{MoSe}_{2}$/$\text{WS}_{2}$ heterobilayers can be
found in Ref. tang2021tuning, . The simulation and the observation share
resemblances but also noticeable differences. Those differences might be
attributed to the simplicity of the present model. Particularly, intervalley
exciton exchange, spin-flip interlayer charge transfer, and electron-phonon
coupling are ignored. Some exciton transition signatures could be lost due to
the omission. However, since it is important to test the limit of the proposed
model and the context of this article is already quite ample, we will leave
those topics to future works.

##  IV Discussions and conclusion

A critical issue we have not discussed is the properness of parameters chosen
in the exciton continuum model. Apart from the parameters also used in
monolayers, four parameters decide the optical spectra of heterobilayers: two
moiré potential depths $V_{1}$, $V_{2}$ and two charge-transfer coupling
strength $w_{\text{e}}$, $w_{\text{h}}$. These parameters can be obtained by
fitting observed optical spectra as what we have done in this work, or
calculated by atomistic simulation. However, there is an inconsistency between
the moiré potential depths acquired from atomistic simulation and fitting the
spectrum. By DFT calculation, moiré potential depths on $\text{MoSe}_{2}$
layer and $\text{WSe}_{2}$ layer of different heterobilayers are estimated to
be around ${8}\sim{10}$ meV (peak-to-peak energy difference around
$80\sim{100}$ meV)naik2022intralayer . Via various simulations including this
work by exciton continuum models, observed moiré-exciton signatures can only
be explained by larger moiré potential depths about $15\sim{35}$
meVjin2019observation ; lin2023remarkably . The discrepancy may be attributed
to the oversimplification of the effective moiré potential, the insufficiency
of DFT calculation to simulate atom-atom potential without sufficient good van
der Waals density functionalsrydberg2003van , or the underestimation of
geometry relaxation and strain in moiré heterobilayerlin2023remarkably . A
theoretical work eliminates the oversimplification by using DFT calculation to
simulate moiré excitons atomisticallynaik2022intralayer . They found that one
of the moiré-exciton signatures can be attributed to the charge-transfer
exciton across different moiré unit cells, and the low moiré potential depth
($\sim{9}$ meV) is sufficient to explain all intralayer moiré-exciton
signatures. However, we do not reach the same conclusion since all radii of
intralayer moiré excitons calculated in this work and listed in Table 2 are
shorter than the moiré lattice constant ($a_{\text{M}}\sim 80$ Å) of
$\text{WSe}_{2}$/$\text{WS}_{2}$ heterobilayers. A scanning tunneling
microscopy experiment on twisted TMDC heterobilayers also indicates that DFT
calculations might underestimate the moiré potential depthsshabani2021deep .
Further investigation on this issue is required in future works.

Note that the exciton continuum model in this study is only applied to the
simulation of the optical spectrum of $\text{WSe}_{2}$/$\text{WS}_{2}$ and
$\text{MoSe}_{2}$/$\text{WS}_{2}$ heterobilayers with small twist angles. This
is because the present model is still too simple to accurately simulate
twisted TMDC homobilayers or twisted TMDC heterobilayers with large twist
angles. Several critical elements are not considered in this model. One is the
valley degree of freedom of TMDCs, which is related to the optical selection
rule of circular-polarized light. Another is the mirror-symmetry-breaking
interaction, which should exist in chiral or helical materials. Without these
elements, the chirality and optical activity of twisted TMDC bilayers can not
be described correctly, and thus the electronic structure of twisted TMDC
bilayers may not be simulated accurately by this model. Additionally, as
mentioned in Sec. III.2, intervalley exciton exchange and spin-flip interlayer
charge transfer are ignored in this model. These effects are also related to
the valley degree of freedom and this ignorance might also contribute to
inaccuracy. Since this article is already lengthy, we decided to keep the
simplicity of this model and omit these elements. In future works, we intend
to include the valley degree of freedom and the mirror-symmetry-breaking
interaction in an improved model and study a broader range of twisted TMDC
bilayers.

Although this exciton continuum model might suffer from oversimplification, it
can be beneficial with its simplicity for future extended studies. For
example, in this work, only lattice mismatch and twist angle on the moiré
pattern are considered, and the rotational symmetry of the system is not
broken. However, in reality, the threefold rotational symmetry of moiré
heterobilayers could be broken in the presence of the uniaxial strain, and the
optical properties of interlayer and hybrid excitons would be affected by the
symmetry breaking. The modification could be done by including a pseudo-gauge
potential in terms of strain tensor componentszheng2021twist ;
zheng2022anomalous . Another example is that exciton complexes such as trions
and biexcitons are known to be computationally exhausting to simulate. It is
easier to use effective-mass models and variational wavefunctions to study
exciton complexesmypaper0 ; mypaper1 . It is a natural extension of our theory
to build the theoretical framework of using continuum models and variational
wavefunctions to study hybridized moiré exciton complexes in multilayer
materials. Additionally, with the future inclusion of intervalley exchange and
electron-phonon coupling into this model, charge and exciton dynamics in moiré
materials could also be studied without too exhausted computation.

In conclusion, we propose a continuum model to study hybridized moiré excitons
in TMDC heterobilayers. A variational wavefunction method is used to solve
this model. We use the present theory to study optical absorption spectra of
$\text{WSe}_{2}$/$\text{WS}_{2}$ and $\text{MoSe}_{2}$/$\text{WS}_{2}$
heterobilayers. A good correspondence between theoretical simulation and
experimental observation of moiré excitons is found. This model could set the
foundation of a theoretical framework to study the physical properties of
moiré excitons in different moiré materials.

###### Acknowledgements.

This work was supported by the National Science and Technology Council, the
Ministry of Education (Higher Education Sprout Project NTU-112L104022), and
the National Center for Theoretical Sciences of Taiwan.

##  Appendix A Modified Rytova-Keldysh potential

The modified Rytova-Keldysh potential for parallel bilayers is derived in this
section. The notation and formulation of the present derivation are taken from
Ref. schwinger1998classical, .

###  A.1 Reduced Green’s function

The Poisson equation for electrostatics in three-dimensional (3D) free space
is written as

$\displaystyle\nabla^{2}\Phi(\mathbf{r})=-{\varrho(\mathbf{r})}/{\kappa},$
(45)

where $\Phi(\mathbf{r})$ is the electrostatic potential, $\varrho(\mathbf{r})$
is the charge density, and $\kappa$ is the dielectric constant. The charge-
free electrostatic potential can be solved by the Laplace equation,
$\nabla^{2}\Phi(\mathbf{r})=0$. The formal solution of the Poisson equation in
three-dimensional free space is given by

$\displaystyle\Phi(\mathbf{r})=\frac{1}{\kappa}\int\frac{\varrho(\mathbf{r}^{\prime})}{|\mathbf{r}-\mathbf{r}^{\prime}|}\text{d}^{3}r.$
(46)

The Green’s function at a point $\mathbf{r}$ interacting with the point charge
$e$ at $\mathbf{r}^{\prime}$ in 3D space is given by

$\displaystyle
G(\mathbf{r},\mathbf{r}^{\prime})=\frac{1}{\kappa|\mathbf{r}-\mathbf{r}^{\prime}|}.$
(47)

Since $\nabla^{2}(1/r)=-4\pi\delta(\mathbf{r})$ and thus

$\displaystyle\kappa\nabla^{2}G(\mathbf{r},\mathbf{r}^{\prime})=-4\pi\delta(\mathbf{r}-\mathbf{r}^{\prime}),$
(48)

Green’s function is a solution of the Poisson equation. The Green’s function
can also be derived from solving Eq. (48). By using the property of the delta
function,

$\displaystyle\delta(\mathbf{r}-\mathbf{r}^{\prime})=\int\frac{\text{d}^{3}k}{(2\pi)^{3}}e^{\mathtt{i}\mathbf{k}\cdot(\mathbf{r}-\mathbf{r}^{\prime})},$
(49)

and the derivative
$\bm{\nabla}e^{\mathtt{i}\mathbf{k}\cdot(\mathbf{r}-\mathbf{r}^{\prime})}=\mathtt{i}\mathbf{k}e^{\mathtt{i}\mathbf{k}\cdot(\mathbf{r}-\mathbf{r}^{\prime})}$,
the Fourier transform of Green’s function can be found as

$\displaystyle
G(\mathbf{r},\mathbf{r}^{\prime})=\int\frac{\text{d}^{3}k}{(2\pi)^{3}}\frac{4\pi}{\kappa{k}^{2}}e^{\mathtt{i}\mathbf{k}\cdot(\mathbf{r}-\mathbf{r}^{\prime})}.$
(50)

Singling out the $z$ direction in a 3D free space, Eq. (50) can be rewritten
as

$\displaystyle G(\mathbf{r},\mathbf{r}^{\prime})$ $\displaystyle=$
$\displaystyle\frac{4\pi}{\kappa}\int\frac{\text{d}k_{x}\text{d}k_{y}}{(2\pi)^{2}}e^{\mathtt{i}\left[k_{x}(x-x^{\prime})+k_{y}(y-y^{\prime})\right]}$
(51)
$\displaystyle\times\int\frac{\text{d}k_{z}}{2\pi}\frac{e^{\mathtt{i}k_{z}(z-z^{\prime})}}{k^{2}_{x}+k^{2}_{y}+k^{2}_{z}}.$

By using

$\displaystyle\int^{\infty}_{-\infty}\frac{\text{d}k_{z}}{2\pi}\frac{e^{\mathtt{i}k_{z}(z-z^{\prime})}}{{k}_{\parallel}^{2}+k^{2}_{z}}=\frac{1}{2{k}_{\parallel}}e^{-{k}_{\parallel}|z-z^{\prime}|},$
(52)

with ${k}_{\parallel}^{2}=k^{2}_{x}+k^{2}_{y}$, the Green’s function becomes

$\displaystyle G(\mathbf{r},\mathbf{r}^{\prime})$ $\displaystyle=$
$\displaystyle{4\pi}\int\frac{\text{d}^{2}{k}_{\parallel}}{(2\pi)^{2}}e^{\mathtt{i}\mathbf{k}_{\parallel}\cdot(\mathbf{r}_{\parallel}-\mathbf{r}^{\prime}_{\parallel})}g(z,z^{\prime};{k}_{\parallel}),$
(53)

where $\mathbf{r}_{\parallel}=(x,y)$ is the 2D coordinates,
$\mathbf{k}_{\parallel}=(k_{x},k_{y})$ is the 2D $k$ vector, and

$\displaystyle
g(z,z^{\prime};{k}_{\parallel})=\frac{1}{2\kappa{k}_{\parallel}}e^{-{k}_{\parallel}|z-z^{\prime}|}$
(54)

is the reduced Green’s function. The differential equation for the reduced
Green’s function can be found as

$\displaystyle\kappa\left[{k}_{\parallel}^{2}-\frac{\partial^{2}}{\partial{z}^{2}}\right]g(z,z^{\prime};{k}_{\parallel})=\delta(z-z^{\prime}).$
(55)

By these two integrals over the infinitesimal segment near $z=z^{\prime}$

$\displaystyle\int^{z^{\prime}+0}_{z^{\prime}-0}\text{d}z\left[{k}_{\parallel}^{2}-\frac{\partial^{2}}{\partial{z}^{2}}\right]g(z,z^{\prime};{k}_{\parallel})=\int^{z^{\prime}+0}_{z^{\prime}-0}\text{d}z\frac{\delta(z-z^{\prime})}{\kappa},$
$\displaystyle\int^{z^{\prime}+0}_{z^{\prime}-0}\text{d}z\;z\left[{k}_{\parallel}^{2}-\frac{\partial^{2}}{\partial{z}^{2}}\right]g(z,z^{\prime};{k}_{\parallel})$
$\displaystyle=\int^{z^{\prime}+0}_{z^{\prime}-0}\text{d}z\;z\frac{\delta(z-z^{\prime})}{\kappa},$
(57)

the first equation gives

$\displaystyle-\frac{\partial{g}}{\partial{z}}\Big{|}_{z=z^{\prime}+0}+\frac{\partial{g}}{\partial{z}}\Big{|}_{z=z^{\prime}-0}=\frac{1}{\kappa},$
(58)

with $0$ indicates an infinitesimal value, and the second equation becomes

$\displaystyle-\int^{z^{\prime}+0}_{z^{\prime}-0}\text{d}z\;z\frac{\partial^{2}}{\partial{z}^{2}}g(z,z^{\prime};{k}_{\parallel})=\frac{z^{\prime}}{\kappa},$
(59)

and then becomes

$\displaystyle-z\frac{\partial{g}}{\partial{z}}\Big{|}^{z^{\prime}+0}_{z^{\prime}-0}+\int^{z^{\prime}+0}_{z^{\prime}-0}\text{d}z\frac{\partial{g}}{\partial{z}}=\frac{z^{\prime}}{\kappa}.$
(60)

By using Eq. (58), it is found that

$\displaystyle g\big{|}_{z=z^{\prime}+0}=g\big{|}_{z=z^{\prime}-0}.$ (61)

The equation indicates the continuity condition for the solution of the
Poisson equation. Eq. (58) and Eq. (61) indicate that the $g$ is continuous
and $\partial{g}/\partial{z}$ is discontinuous at $z=z^{\prime}$.

###  A.2 Rytova-Keldysh potential in parallel bilayers

The parallel bilayers separate three regimes. Two are uniform dielectric
materials $z<0$, $z>d$. The regime in the middle $0<z<d$ is an empty space.
The Green’s function is written as

$\displaystyle z\leq{0}:\hskip
2.15277pt\kappa_{1}\nabla^{2}G(\mathbf{r},\mathbf{r}^{\prime})=-4\pi[\delta(\mathbf{r}-\mathbf{r}^{\prime})+\sigma_{1}\delta(z)],$
(62) $\displaystyle 0<z<d:\hskip
2.15277pt\kappa_{0}\nabla^{2}G(\mathbf{r},\mathbf{r}^{\prime})=-4\pi\delta(\mathbf{r}-\mathbf{r}^{\prime}),$
(63) $\displaystyle z\geq{d}:\hskip
2.15277pt\kappa_{2}\nabla^{2}G(\mathbf{r},\mathbf{r}^{\prime})=-4\pi[\delta(\mathbf{r}-\mathbf{r}^{\prime})+\sigma_{2}\delta(z-d)],$

where $\sigma_{1}$, $\sigma_{2}$ are the surface bound charges on the layers,
$\kappa_{1}$, $\kappa_{2}$ are the dielectric constants of the uniform
dielectric materials, and $\kappa_{0}$ is the dielectric constant for the
empty space. The surface-bound charge can be related to the polarization of
the layer by

$\displaystyle\bm{\nabla}_{\parallel}\cdot\mathbf{P}_{\parallel}=-\sigma,$
(65)

where $\bm{\nabla}_{\parallel}=(\partial/\partial{x},\partial/\partial{y})$
and $\mathbf{P}_{\parallel}=(P_{x},P_{y})$ are the divergence and the
polarization on the 2D space. The polarization can be related to the electric
field by

$\displaystyle\mathbf{P}_{\parallel}=\chi_{\text{2D}}\mathbf{E}_{\parallel},$
(66)

with $\mathbf{E}_{\parallel}=(E_{x},E_{y})$ the electric field on the two-
dimensional space and $\chi_{\text{2D}}$ is the electric susceptibility (or
polarizability) of the layer. Since
$\mathbf{E}_{\parallel}=-\bm{\nabla}_{\parallel}\Phi$, we can derive the 2D
Poisson equation,

$\displaystyle\chi_{1}\nabla^{2}_{\parallel}\Phi(\mathbf{r})\big{|}_{z=0}=\sigma_{1},\hskip
8.61108pt\chi_{2}\nabla^{2}_{\parallel}\Phi(\mathbf{r})\big{|}_{z=d}=\sigma_{2},$
(67)

where $\chi_{1}$ and $\chi_{2}$ are the electric susceptibilities of the
corresponding layers, the Green’s function can be rewritten as

$\displaystyle z\leq{0}:$
$\displaystyle\kappa_{1}\left[\nabla^{2}+\frac{4\pi}{\kappa_{1}}\chi_{1}\delta(z)\nabla^{2}_{\parallel}\right]G(\mathbf{r},\mathbf{r}^{\prime})=-4\pi\delta(\mathbf{r}-\mathbf{r}^{\prime}),$
(68) $\displaystyle 0<z<d:\hskip
4.30554pt\kappa_{0}\nabla^{2}G(\mathbf{r},\mathbf{r}^{\prime})=-4\pi\delta(\mathbf{r}-\mathbf{r}^{\prime}),$
(69) $\displaystyle z\geq{d}:$
$\displaystyle\kappa_{2}\left[\nabla^{2}+\frac{4\pi}{\kappa_{2}}\chi_{2}\delta(z-d)\nabla^{2}_{\parallel}\right]G(\mathbf{r},\mathbf{r}^{\prime})=-4\pi\delta(\mathbf{r}-\mathbf{r}^{\prime}).$

By using Eq. (53), the Poisson equation for the reduced Green’s function can
be written as

$\displaystyle z\leq{0}:\hskip
4.30554pt\kappa_{1}\left[{k}_{\parallel}^{2}-\frac{\partial^{2}}{\partial{z}^{2}}+\frac{4\pi\chi_{1}}{\kappa_{1}}{k}_{\parallel}^{2}\delta(z)\right]g=\delta(z-z^{\prime}),$
(71) $\displaystyle 0<z<d:\hskip
4.30554pt\kappa_{0}\left[{k}_{\parallel}^{2}-\frac{\partial^{2}}{\partial{z}^{2}}\right]g=\delta(z-z^{\prime}),$
(72) $\displaystyle z\geq{d}:$
$\displaystyle\kappa_{2}\left[{k}_{\parallel}^{2}-\frac{\partial^{2}}{\partial{z}^{2}}+\frac{4\pi\chi_{2}}{\kappa_{2}}{k}_{\parallel}^{2}\delta(z-d)\right]g=\delta(z-z^{\prime}).$
(73)

Assuming that the test charge is at $z^{\prime}=0$, the boundary conditions
can be derived

$\displaystyle-\kappa_{0}\frac{\partial}{\partial{z}}g\Big{|}_{z=+0}+\kappa_{1}\frac{\partial}{\partial{z}}g\Big{|}_{z=-0}=1-4\pi\chi_{1}{k}_{\parallel}^{2}g\big{|}_{z=0},$
$\displaystyle-\kappa_{2}\frac{\partial}{\partial{z}}g\Big{|}_{z=d+0}+\kappa_{0}\frac{\partial}{\partial{z}}g\Big{|}_{z=d-0}=-4\pi\chi_{2}{k}_{\parallel}^{2}g\big{|}_{z=d}.$

The general solution of Green’s function is assumed as

$\displaystyle z\leq{0}:$ $\displaystyle g=Ae^{{k}_{\parallel}z},$ (76)
$\displaystyle 0<z<d:$ $\displaystyle
g=Be^{{k}_{\parallel}z}+Ce^{-{k}_{\parallel}z},$ (77) $\displaystyle
z\geq{d}:$ $\displaystyle g=De^{-{k}_{\parallel}z},$ (78)

By using the continuity condition in Eq. (61) and the boundary contions, the
four coefficients can be connected to four coupled linear equations,

$\displaystyle A=B+C,\hskip
8.61108ptBe^{{k}_{\parallel}d}+Ce^{-{k}_{\parallel}d}=De^{-{k}_{\parallel}d}.$
(79)
$\displaystyle\left(\kappa_{1}{k}_{\parallel}+2\xi_{1}{k}_{\parallel}^{2}\right)A-\kappa_{0}{k}_{\parallel}B+\kappa_{0}{k}_{\parallel}C=1,$
(80)
$\displaystyle\kappa_{0}{k}_{\parallel}\left(Be^{{k}_{\parallel}{d}}-Ce^{-{k}_{\parallel}{d}}\right)$
$\displaystyle+\left(\kappa_{2}{k}_{\parallel}+4\pi\chi_{2}{k}_{\parallel}^{2}\right)De^{-{k}_{\parallel}{d}}=0.$
(81)

By using the first equation $A=B+C$ in other linear equations to replace the
coefficient $A$, we can reduce one of the four variables. We get the three
coupled equations,

$\displaystyle
Be^{{k}_{\parallel}d}+Ce^{-{k}_{\parallel}d}-De^{-{k}_{\parallel}d}=0.$ (82)
$\displaystyle(\lambda_{1}-\kappa_{0})B+(\lambda_{1}+\kappa_{0})C=1/{k}_{\parallel},$
(83)
$\displaystyle\kappa_{0}Be^{{k}_{\parallel}{d}}-\kappa_{0}Ce^{-{k}_{\parallel}{d}}+\lambda_{2}De^{-{k}_{\parallel}{d}}=0.$
(84)

with $\lambda_{1}=\kappa_{1}+4\pi\chi_{1}{k}_{\parallel}$ and
$\lambda_{2}=\kappa_{2}+4\pi\chi_{2}{k}_{\parallel}$. These equations can be
rewritten as the matrix form,

$\displaystyle\begin{bmatrix}\lambda_{1}-\kappa_{0}&\lambda_{1}+\kappa_{0}&0\\\
e^{{k}_{\parallel}d}&e^{-{k}_{\parallel}d}&-e^{-{k}_{\parallel}d}\\\
\kappa_{0}e^{{k}_{\parallel}{d}}&-\kappa_{0}e^{-{k}_{\parallel}{d}}&\lambda_{2}e^{-{k}_{\parallel}{d}}\end{bmatrix}\begin{bmatrix}B\\\
C\\\ D\end{bmatrix}=\begin{bmatrix}1/{k}_{\parallel}\\\ 0\\\
0\end{bmatrix},\hskip 12.91663pt$ (85)

By Cramer’s rule, these coefficients can be solved as

$\displaystyle B$ $\displaystyle=$
$\displaystyle\frac{1}{\textit{Det}}\begin{vmatrix}1/{k}_{\parallel}&\lambda_{1}+\kappa_{0}&0\\\
0&e^{-{k}_{\parallel}d}&-e^{-{k}_{\parallel}d}\\\
0&-\kappa_{0}e^{-{k}_{\parallel}{d}}&\lambda_{2}e^{-{k}_{\parallel}{d}}\end{vmatrix}$
(86) $\displaystyle=$
$\displaystyle\frac{(\lambda_{2}-\kappa_{0})e^{-2{k}_{\parallel}{d}}}{{k}_{\parallel}\textit{Det}},$
$\displaystyle C$ $\displaystyle=$
$\displaystyle\frac{1}{\textit{Det}}\begin{vmatrix}\lambda_{1}-\kappa_{0}&1/{k}_{\parallel}&0\\\
e^{{k}_{\parallel}d}&0&-e^{-{k}_{\parallel}d}\\\
\kappa_{0}e^{{k}_{\parallel}{d}}&0&\lambda_{2}e^{-{k}_{\parallel}{d}}\end{vmatrix}=-\frac{\lambda_{2}+\kappa_{0}}{{k}_{\parallel}\textit{Det}},\hskip
12.91663pt$ (87) $\displaystyle D$ $\displaystyle=$
$\displaystyle\frac{1}{\textit{Det}}\begin{vmatrix}\lambda_{1}-\kappa_{0}&\lambda_{1}+\kappa_{0}&1/{k}_{\parallel}\\\
e^{{k}_{\parallel}d}&e^{-{k}_{\parallel}d}&0\\\
\kappa_{0}e^{{k}_{\parallel}{d}}&-\kappa_{0}e^{-{k}_{\parallel}{d}}&0\end{vmatrix}=-\frac{2\kappa_{0}}{{k}_{\parallel}\textit{Det}},\hskip
12.91663pt$ (88) Det $\displaystyle=$
$\displaystyle\begin{vmatrix}\lambda_{1}-\kappa_{0}&\lambda_{1}+\kappa_{0}&0\\\
e^{{k}_{\parallel}d}&e^{-{k}_{\parallel}d}&-e^{-{k}_{\parallel}d}\\\
\kappa_{0}e^{{k}_{\parallel}{d}}&-\kappa_{0}e^{-{k}_{\parallel}{d}}&\lambda_{2}e^{-{k}_{\parallel}{d}}\end{vmatrix}$
(89) $\displaystyle=$
$\displaystyle(\lambda_{1}-\kappa_{0})(\lambda_{2}-\kappa_{0})e^{-2{k}_{\parallel}{d}}$
$\displaystyle-(\lambda_{1}+\kappa_{0})(\lambda_{2}+\kappa_{0}).$

Therefore, the intralayer and interlayer Rytova-Keldysh potentials can be
found as

$\displaystyle W_{11}({k}_{\parallel})=4\pi g\big{|}_{z=0,\;z^{\prime}=0}$
$\displaystyle=\frac{(4\pi/{k}_{\parallel})\left[(\lambda_{2}+\kappa_{0})e^{{k}_{\parallel}{d}}-(\lambda_{2}-\kappa_{0})e^{-{k}_{\parallel}{d}}\right]}{(\lambda_{1}+\kappa_{0})(\lambda_{2}+\kappa_{0})e^{{k}_{\parallel}{d}}-(\lambda_{1}-\kappa_{0})(\lambda_{2}-\kappa_{0})e^{-{k}_{\parallel}{d}}},$
$\displaystyle W_{21}({k}_{\parallel})=4\pi g\big{|}_{z=d,\;z^{\prime}=0}$
$\displaystyle=\frac{8\pi\kappa_{0}/{k}_{\parallel}}{(\lambda_{1}+\kappa_{0})(\lambda_{2}+\kappa_{0})e^{{k}_{\parallel}{d}}-(\lambda_{1}-\kappa_{0})(\lambda_{2}-\kappa_{0})e^{-{k}_{\parallel}{d}}}.$

By replacing $q=k_{\parallel}$ and using dielectric functions, these
potentials can be rewritten as

$\displaystyle W_{11}(q)=\frac{2\pi}{\epsilon_{11}(q)q},\hskip
4.30554ptW_{21}(q)=\frac{2\pi}{\epsilon_{21}(q)q},$ (92)

where the intralayer and interlayer dielectric functions are written as

$\displaystyle\epsilon_{11}(q)$ $\displaystyle=$
$\displaystyle\frac{\kappa_{0}\epsilon_{21}(q)}{\left(\frac{\kappa_{2}+\kappa_{0}}{2}+r_{2}q\right)e^{q{d}}-\left(\frac{\kappa_{2}-\kappa_{0}}{2}+r_{2}q\right)e^{-q{d}}},$
$\displaystyle\epsilon_{21}(q)$ $\displaystyle=$
$\displaystyle\left(\frac{\kappa_{1}+\kappa_{0}}{2}+r_{1}q\right)\left(\frac{\kappa_{2}+\kappa_{0}}{2}+r_{2}q\right)\frac{e^{q{d}}}{\kappa_{0}}$
$\displaystyle-\left(\frac{\kappa_{1}-\kappa_{0}}{2}+r_{1}q\right)\left(\frac{\kappa_{2}-\kappa_{0}}{2}+r_{2}q\right)\frac{e^{-q{d}}}{\kappa_{0}},$

with $r_{1}=2\pi\chi_{1}$ and $r_{2}=2\pi\chi_{2}$ being the screening
lengths. By assigning $\kappa_{1}=\kappa_{2}$, the modified Rytova-Keldysh
potentials in Eq. (9), Eq. (10), and Eq. (LABEL:inter) can be derived.

##  Appendix B Matrix elements of exciton Hamiltonian

Matrix elements of the exciton Hamiltonian matrix and overlap matrix expanded
by the basis functions in Sec. II.4 are derived. By using the variational
method, the exciton wavefunction can be solved from the matrix eigenvalue
equation

$\displaystyle\sum_{l^{\prime}_{\text{e}}l^{\prime}_{\text{h}},b}\bar{\mathcal{H}}_{(l_{\text{e}}l_{\text{h}},{a},\mathbf{G})(l^{\prime}_{\text{e}}l^{\prime}_{\text{h}},{b},\mathbf{G}^{\prime}),\mathbf{K}}C_{(l^{\prime}_{\text{e}}l^{\prime}_{\text{h}},{b},\mathbf{G}^{\prime}),I\mathbf{K}}=\varepsilon_{\text{X},I\mathbf{K}}\sum_{b}\mathcal{O}_{(l_{\text{e}}l_{\text{h}},a,\mathbf{G})(l_{\text{e}}l_{\text{h}},b,\mathbf{G}^{\prime})}C_{(l_{\text{e}}l_{\text{h}},{b},\mathbf{G}^{\prime}),I\mathbf{K}}.$
(95)

The exciton Hamiltonian matrix elements are given by

$\displaystyle\bar{\mathcal{H}}_{(l_{\text{e}}l_{\text{h}},{a},\mathbf{G})(l^{\prime}_{\text{e}}l^{\prime}_{\text{h}},{b},\mathbf{G}^{\prime}),\mathbf{K}}$
$\displaystyle=$
$\displaystyle\delta_{l_{\text{e}},l^{\prime}_{\text{e}}}\delta_{l_{\text{h}},l^{\prime}_{\text{h}}}\int\psi^{*}_{\mathbf{G}}(\mathbf{R})\phi^{*}_{l_{\text{e}}l_{\text{h}},a}(\mathbf{r})\tilde{\mathcal{H}}_{l_{\text{e}}l_{\text{h}}}(\mathbf{R},\mathbf{r})\psi_{\mathbf{G}^{\prime}}(\mathbf{R})\phi_{l_{\text{e}}l_{\text{h}},b}(\mathbf{r})\text{d}^{2}R\text{d}^{2}r$
(96)
$\displaystyle+(1-\delta_{l_{\text{e}},l^{\prime}_{\text{e}}})\delta_{l_{\text{h}},l^{\prime}_{\text{h}}}\int\psi^{*}_{\mathbf{G}}(\mathbf{R})\phi^{*}_{l_{\text{e}}l_{\text{h}},a}(\mathbf{r})\tilde{\mathcal{T}}_{\text{e},(l_{\text{e}}l_{\text{h}})(l^{\prime}_{\text{e}}l_{\text{h}})}(\mathbf{R},\mathbf{r})\psi_{\mathbf{G}^{\prime}}(\mathbf{R})\phi_{l^{\prime}_{\text{e}}l_{\text{h}},b}(\mathbf{r})\text{d}^{2}R\text{d}^{2}r$
$\displaystyle+\delta_{l_{\text{e}},l^{\prime}_{\text{e}}}(1-\delta_{l_{\text{h}},l^{\prime}_{\text{h}}})\int\psi^{*}_{\mathbf{G}}(\mathbf{R})\phi^{*}_{l_{\text{e}}l_{\text{h}},a}(\mathbf{r})\tilde{\mathcal{T}}_{\text{h},(l_{\text{e}}l_{\text{h}})(l_{\text{e}}l^{\prime}_{\text{h}})}(\mathbf{R},\mathbf{r})\psi_{\mathbf{G}^{\prime}}(\mathbf{R})\phi_{l_{\text{e}}l^{\prime}_{\text{h}},b}(\mathbf{r})\text{d}^{2}R\text{d}^{2}r,$

with

$\displaystyle\tilde{\mathcal{H}}_{l_{\text{e}}l_{\text{h}}}(\mathbf{R},\mathbf{r})=e^{-\mathtt{i}\mathbf{K}\cdot\mathbf{R}}\exp\left(-{\mathtt{i}\bm{\kappa}_{l_{\text{e}}l_{\text{h}}}\cdot\mathbf{r}}\right)\mathcal{H}_{l_{\text{e}}l_{\text{h}}}(\mathbf{R},\mathbf{r})\exp\left({\mathtt{i}\bm{\kappa}_{l_{\text{e}}l_{\text{h}}}\cdot\mathbf{r}}\right)e^{\mathtt{i}\mathbf{K}\cdot\mathbf{R}},$
(97)
$\displaystyle\tilde{\mathcal{T}}_{\text{e},(1l_{\text{h}})(2l_{\text{h}})}(\mathbf{R},\mathbf{r})=e^{-\mathtt{i}\mathbf{K}\cdot\mathbf{R}}\exp\left(-{\mathtt{i}\bm{\kappa}_{1l_{\text{h}}}\cdot\mathbf{r}}\right)\mathcal{T}_{\text{e}}(\mathbf{R},\mathbf{r})\exp\left({\mathtt{i}\bm{\kappa}_{2l_{\text{h}}}\cdot\mathbf{r}}\right)e^{\mathtt{i}\mathbf{K}\cdot\mathbf{R}},$
(98)
$\displaystyle\tilde{\mathcal{T}}_{\text{e},(2l_{\text{h}})(1l_{\text{h}})}(\mathbf{R},\mathbf{r})=e^{-\mathtt{i}\mathbf{K}\cdot\mathbf{R}}\exp\left(-{\mathtt{i}\bm{\kappa}_{2l_{\text{h}}}\cdot\mathbf{r}}\right)\mathcal{T}^{*}_{\text{e}}(\mathbf{R},\mathbf{r})\exp\left({\mathtt{i}\bm{\kappa}_{1l_{\text{h}}}\cdot\mathbf{r}}\right)e^{\mathtt{i}\mathbf{K}\cdot\mathbf{R}},$
(99)
$\displaystyle\tilde{\mathcal{T}}_{\text{h},(l_{\text{e}}1)(l_{\text{e}}2)}(\mathbf{R},\mathbf{r})=e^{-\mathtt{i}\mathbf{K}\cdot\mathbf{R}}\exp\left(-{\mathtt{i}\bm{\kappa}_{l_{\text{e}}1}\cdot\mathbf{r}}\right)\mathcal{T}_{\text{h}}(\mathbf{R},\mathbf{r})\exp\left({\mathtt{i}\bm{\kappa}_{l_{\text{e}}2}\cdot\mathbf{r}}\right)e^{\mathtt{i}\mathbf{K}\cdot\mathbf{R}},$
(100)
$\displaystyle\tilde{\mathcal{T}}_{\text{h},(l_{\text{e}}2)(l_{\text{e}}1)}(\mathbf{R},\mathbf{r})=e^{-\mathtt{i}\mathbf{K}\cdot\mathbf{R}}\exp\left(-{\mathtt{i}\bm{\kappa}_{l_{\text{e}}2}\cdot\mathbf{r}}\right)\mathcal{T}^{*}_{\text{h}}(\mathbf{R},\mathbf{r})\exp\left({\mathtt{i}\bm{\kappa}_{l_{\text{e}}1}\cdot\mathbf{r}}\right)e^{\mathtt{i}\mathbf{K}\cdot\mathbf{R}}.$
(101)

The overlap matrix elements are given by

$\displaystyle\bar{\mathcal{O}}_{(l_{\text{e}}l_{\text{h}},{a},\mathbf{G})(l_{\text{e}}l_{\text{h}},{b},\mathbf{G}^{\prime})}=\int\psi^{*}_{\mathbf{G}}(\mathbf{R})\phi^{*}_{l_{\text{e}}l_{\text{h}},a}(\mathbf{r})\psi_{\mathbf{G}^{\prime}}(\mathbf{R})\phi_{l_{\text{e}}l_{\text{h}},b}(\mathbf{r})\text{d}^{2}R\text{d}^{2}r.$
(102)

The diagonal Hamiltonian matrix elements are given by

$\displaystyle\bar{\mathcal{H}}_{(l_{\text{e}}l_{\text{h}},{a},\mathbf{G})(l_{\text{e}}l_{\text{h}},{b},\mathbf{G}^{\prime}),\mathbf{K}}$
$\displaystyle=$
$\displaystyle\delta_{\mathbf{G},\mathbf{G}^{\prime}}o_{ab}\left[\Delta_{l_{\text{e}}l_{\text{h}}}+\frac{|\mathbf{K}-\mathbf{G}-\bm{\mathcal{K}}_{l_{\text{e}}l_{\text{h}}}|^{2}}{2m_{\text{X},l_{\text{e}}l_{\text{h}}}}\right]+\delta_{\mathbf{G},\mathbf{G}^{\prime}}\langle{a}|\left[-\frac{\nabla^{2}}{2\mu_{\text{X},l_{\text{e}}l_{\text{h}}}}-W_{l_{\text{e}}l_{\text{h}}}(r)\right]|{b}\rangle$
(103)
$\displaystyle+\langle{a,\mathbf{G}}|\mathcal{V}_{l_{\text{e}}l_{\text{h}}}(\mathbf{R},\mathbf{r})|{b,\mathbf{G}^{\prime}}\rangle,$
$\displaystyle-\frac{1}{2\mu_{\text{X},l_{\text{e}}l_{\text{h}}}}\langle{a}|\nabla^{2}|{b}\rangle$
$\displaystyle=$
$\displaystyle-\frac{1}{2\mu_{\text{X},l_{\text{e}}l_{\text{h}}}}\frac{(N_{a}+N_{b}-1)!}{\left(\mathcal{Z}_{l_{\text{e}}l_{\text{h}},a}+\mathcal{Z}_{l_{\text{e}}l_{\text{h}},b}\right)^{N_{a}+N_{b}}}\Bigg{\\{}(1-\delta_{N_{b},1})\frac{\left[(N_{b}-1)^{2}-L_{b}^{2}\right]\left(\mathcal{Z}_{l_{\text{e}}l_{\text{h}},a}+\mathcal{Z}_{l_{\text{e}}l_{\text{h}},b}\right)^{2}}{(N_{a}+N_{b}-1)(N_{a}+N_{b}-2)}$
(104)
$\displaystyle-\frac{\left[(2N_{b}-1)\mathcal{Z}_{l_{\text{e}}l_{\text{h}},b}\right]\left(\mathcal{Z}_{l_{\text{e}}l_{\text{h}},a}+\mathcal{Z}_{l_{\text{e}}l_{\text{h}},b}\right)}{(N_{a}+N_{b}-1)}+\mathcal{Z}^{2}_{l_{\text{e}}l_{\text{h}},b}\Bigg{\\}},$
$\displaystyle\langle{a}|W_{l_{\text{e}}l_{\text{h}}}(\mathbf{r})|{b}\rangle$
$\displaystyle=$
$\displaystyle\frac{\delta_{L_{a},L_{b}}}{(2\pi)^{2}}\int^{\infty}_{0}\tilde{\mathcal{R}}_{N_{a}+N_{b}-1,0}(\mathcal{Z}_{l_{\text{e}}l_{\text{h}},a}+\mathcal{Z}_{l_{\text{e}}l_{\text{h}},b},k)\tilde{W}_{l_{\text{e}}l_{\text{h}}}(k)k\text{d}k,$
(105)
$\displaystyle\langle{a,\mathbf{G}}|\mathcal{V}_{l_{\text{e}}l_{\text{h}}}(\mathbf{R},\mathbf{r})|{b,\mathbf{G}^{\prime}}\rangle$
$\displaystyle=$
$\displaystyle\sum_{j=1,3,5}\Big{\\{}V_{l_{\text{e}}}\exp\left[\mathtt{i}(-1)^{l_{\text{e}}}\psi\right]\delta(\mathbf{G}-\mathbf{G}^{\prime}+\bm{g}_{j})\tilde{\rho}_{(l_{\text{e}}l_{\text{h}},a)(l_{\text{e}}l_{\text{h}},b)}(-\gamma_{\text{h},l_{\text{h}}}\bm{g}_{j})$
(106)
$\displaystyle+V_{l_{\text{e}}}\exp\left[-\mathtt{i}(-1)^{l_{\text{e}}}\psi\right]\delta(\mathbf{G}-\mathbf{G}^{\prime}-\bm{g}_{j})\tilde{\rho}_{(l_{\text{e}}l_{\text{h}},a)(l_{\text{e}}l_{\text{h}},b)}(\gamma_{\text{h},l_{\text{h}}}\bm{g}_{j})$
$\displaystyle-
V_{l_{\text{h}}}\exp\left[\mathtt{i}(-1)^{l_{\text{h}}}\psi\right]\delta(\mathbf{G}-\mathbf{G}^{\prime}+\bm{g}_{j})\tilde{\rho}_{(l_{\text{e}}l_{\text{h}},a)(l_{\text{e}}l_{\text{h}},b)}(\gamma_{\text{e},l_{\text{e}}}\bm{g}_{j})$
$\displaystyle-
V_{l_{\text{h}}}\exp\left[-\mathtt{i}(-1)^{l_{\text{h}}}\psi\right]\delta(\mathbf{G}-\mathbf{G}^{\prime}-\bm{g}_{j})\tilde{\rho}_{(l_{\text{e}}l_{\text{h}},a)(l_{\text{e}}l_{\text{h}},b)}(-\gamma_{\text{e},l_{\text{e}}}\bm{g}_{j})\Big{\\}},$

with

$\displaystyle\tilde{\rho}_{(l_{\text{e}}l_{\text{h}},a)(l^{\prime}_{\text{e}}l^{\prime}_{\text{h}},b)}(\mathbf{k})$
$\displaystyle=$
$\displaystyle\frac{e^{-\mathtt{i}(L_{a}-L_{b})\varphi_{\mathbf{k}}}}{2\pi}\tilde{\mathcal{R}}_{N_{a}+N_{b}-1,L_{b}-L_{a}}(\mathcal{Z}_{l_{\text{e}}l_{\text{h}},a}+\mathcal{Z}_{l^{\prime}_{\text{e}}l^{\prime}_{\text{h}},b},k),$
(107)

where the radial function in momentum space can be obtained by the generating
formula

$\displaystyle\tilde{\mathcal{R}}_{N,L}(\mathcal{Z},k)$ $\displaystyle=$
$\displaystyle\frac{2\pi(-\mathtt{i})^{N}}{k^{N+1}}\frac{\text{d}^{N}}{\text{d}z^{N}}\frac{\left(z-\mathtt{i}(L/|L|)\sqrt{1-z^{2}}\right)^{|L|}}{\sqrt{1-z^{2}}}\Bigg{|}_{z=\mathtt{i}{\mathcal{Z}}/{k}}.$
(108)

The nondiagonal Hamiltonian matrix elements can be rewritten as

$\displaystyle\bar{\mathcal{H}}_{(l_{\text{e}}l_{\text{h}},a,\mathbf{G})(l^{\prime}_{\text{e}}l^{\prime}_{\text{h}},b,\mathbf{G}^{\prime}),\mathbf{K}}$
$\displaystyle=$ $\displaystyle
w_{\text{e}}\delta_{l_{\text{h}},l^{\prime}_{\text{h}}}(1-\delta_{l_{\text{e}},l^{\prime}_{\text{e}}})\Big{[}\delta(\mathbf{G}-\mathbf{G}^{\prime})\tilde{\rho}_{(l_{\text{e}}l_{\text{h}},a)(l^{\prime}_{\text{e}}l_{\text{h}},b)}\big{(}(l_{\text{e}}-l^{\prime}_{\text{e}})(\gamma_{\text{h},l_{\text{h}}}\bm{\kappa}_{3}+\Delta_{\text{e}}\bm{\kappa}_{l_{\text{h}}})\big{)}$
(109)
$\displaystyle+\delta\big{(}\mathbf{G}-\mathbf{G}^{\prime}-(l_{\text{e}}-l^{\prime}_{\text{e}})\bm{g}_{1}\big{)}\tilde{\rho}_{(l_{\text{e}}l_{\text{h}},a)(l^{\prime}_{\text{e}}l_{\text{h}},b)}\big{(}(l_{\text{e}}-l^{\prime}_{\text{e}})(\gamma_{\text{h},l_{\text{h}}}\bm{\kappa}_{1}+\Delta_{\text{e}}\bm{\kappa}_{l_{\text{h}}})\big{)}$
$\displaystyle+\delta\big{(}\mathbf{G}-\mathbf{G}^{\prime}-(l_{\text{e}}-l^{\prime}_{\text{e}})\bm{g}_{2}\big{)}\tilde{\rho}_{(l_{\text{e}}l_{\text{h}},a)(l^{\prime}_{\text{e}}l_{\text{h}},b)}\big{(}(l_{\text{e}}-l^{\prime}_{\text{e}})(\gamma_{\text{h},l_{\text{h}}}\bm{\kappa}_{5}+\Delta_{\text{e}}\bm{\kappa}_{l_{\text{h}}})\big{)}\Big{]}$
$\displaystyle+w_{\text{h}}\delta_{l_{\text{e}},l^{\prime}_{\text{e}}}(1-\delta_{l_{\text{h}},l^{\prime}_{\text{h}}})\Big{[}\delta(\mathbf{G}-\mathbf{G}^{\prime})\tilde{\rho}_{(l_{\text{e}}l_{\text{h}},a)(l_{\text{e}}l^{\prime}_{\text{h}},b)}\big{(}(l_{\text{h}}-l^{\prime}_{\text{h}})(\gamma_{\text{e},l_{\text{e}}}\bm{\kappa}_{3}+\Delta_{\text{h}}\bm{\kappa}_{l_{\text{e}}})\big{)}$
$\displaystyle+\delta\big{(}\mathbf{G}-\mathbf{G}^{\prime}+(l_{\text{h}}-l^{\prime}_{\text{h}})\bm{g}_{1}\big{)}\tilde{\rho}_{(l_{\text{e}}l_{\text{h}},a)(l_{\text{e}}l^{\prime}_{\text{h}},b)}\big{(}(l_{\text{h}}-l^{\prime}_{\text{h}})(\gamma_{\text{e},l_{\text{e}}}\bm{\kappa}_{1}+\Delta_{\text{h}}\bm{\kappa}_{l_{\text{e}}})\big{)}$
$\displaystyle+\delta\big{(}\mathbf{G}-\mathbf{G}^{\prime}+(l_{\text{h}}-l^{\prime}_{\text{h}})\bm{g}_{2}\big{)}\tilde{\rho}_{(l_{\text{e}}l_{\text{h}},a)(l_{\text{e}}l^{\prime}_{\text{h}},b)}\big{(}(l_{\text{h}}-l^{\prime}_{\text{h}})(\gamma_{\text{e},l_{\text{e}}}\bm{\kappa}_{5}+\Delta_{\text{h}}\bm{\kappa}_{l_{\text{e}}})\big{)}\Big{]},$

with $\Delta_{\text{e}}=\gamma_{\text{e},2}-\gamma_{\text{e},1}$,
$\Delta_{\text{h}}=\gamma_{\text{h},2}-\gamma_{\text{h},1}$. The overlap
matrix element is given by

$\displaystyle\bar{\mathcal{O}}_{(l_{\text{e}}l_{\text{h}},{a},\mathbf{G})(l_{\text{e}}l_{\text{h}},{b},\mathbf{G}^{\prime})}$
$\displaystyle=$
$\displaystyle\delta(\mathbf{G}-\mathbf{G}^{\prime})\delta_{L_{a},L_{b}}\frac{(N_{a}+N_{b}-1)!}{(\mathcal{Z}_{l_{\text{e}}l_{\text{h}},a}+\mathcal{Z}_{l_{\text{e}}l_{\text{h}},b})^{N_{a}+N_{b}}}.$
(110)

A more detailed derivation of matrix elements of the exciton Hamiltonian by
using STOs as the basis function can be found in Ref. mypaper0, ; mypaper1, .

## References

  * (1) Q. H. Wang, K. Kalantar-Zadeh, A. Kis, J. N. Coleman, and M. S. Strano, Electronics and optoelectronics of two-dimensional transition metal dichalcogenides, Nat. Nanotech., 7, 699 (2012).
  * (2) F. Xia, H. Wang, D. Xiao, M. Dubey, and A. Ramasubramaniam, Two-dimensional material nanophotonics, Nat. Photon., 8, 899 (2014).
  * (3) K. F. Mak and J. Shan, Photonics and optoelectronics of 2D semiconductor transition metal dichalcogenides, Nat. Photon., 10, 216 (2016).
  * (4) J. Xiao, M. Zhao, Y. Wang, and X. Zhang, Excitons in atomically thin 2D semiconductors and their applications, Nanophotonics, 6, 1309 (2017).
  * (5) T. C. Berkelbach, and D. R. Reichman, Optical and Excitonic Properties of Atomically Thin Transition-Metal Dichalcogenides Ann. Rev. Condens. Matter Phys., 9, 379 (2018).
  * (6) K. Tran, J. Choi, and A. Singh, Moiré and beyond in transition metal dichalcogenide twisted bilayers, 2D Mater., 8, 022002 (2021).
  * (7) D. M. Kennes, M. Claassen, L. Xian, A. Georges, A. J. Millis, J. Hone, C. R. Dean, D. N. Basov, A. N. Pasupathy, and A. Rubio Moiré heterostructures as a condensed-matter quantum simulator, Nat. Phys., 17, 155 (2021).
  * (8) K. F. Mak and J. Shan, Semiconductor moiré materials, Nat. Nanotech., 17, 686 (2022).
  * (9) L. Zhang, R. Ni, and Y. Zhou, Controlling quantum phases of electrons and excitons in moiré superlattices, J. Appl. Phys., 133, 080901 (2023).
  * (10) L. Du, M. R. Molas, Z. Huang, G. Zhang, F. Wang, and Z. Sun, Moiré photonics and optoelectronics, Science, 379, eadg0014 (2023).
  * (11) J. Kunstmann, F. Mooshammer, P. Nagler, A. Chaves, F. Stein, N. Paradiso, G. Plechinger, C. Strunk, C. Sch’́uller, G. Seifert, D. R. Reichman, and T. Korn, Momentum-space indirect interlayer excitons in transition-metal dichalcogenide van der Waals heterostructures, Nat. Phys., 14, 801 (2018).
  * (12) P. Rivera, H. Yu, K. L. Seyler, N. P. Wilson, W. Yao, and X. Xu, Interlayer valley excitons in heterobilayers of transition metal dichalcogenides, Nat. Nanotech., 13, 1004 (2018).
  * (13) E. M. Alexeev, D. A. Ruiz-Tijerina, M. Danovich, M. J. Hamer, D. J. Terry, P. K. Nayak, S. Ahn, S. Pak, J. Lee, J. I. Sohn, M. R. Molas, M. Koperski, K. Watanabe, T. Taniguchi, K. S. Novoselov, R. V. Gorbachev, H. S. Shin, V. I. Fal’ko, and A. I. Tartakovskii, Resonantly hybridized excitons in moiré superlattices in van der Waals heterostructures, Nature, 567, 81 (2019).
  * (14) C. Jin, E. C. Regan, A. Yan, M. I. B. Utama, D. Wang, S. Zhao, Y. Qin, S. Yang, Z. Zheng, S. Shi, K. Watanabe, T. Taniguchi, S. Tongay, A. Zettl, and F. Wang, Observation of moiré excitons in WSe2/WS2 heterostructure superlattices, Nature, 567, 76 (2019).
  * (15) C. Jin, E. C. Regan, D. Wang, M. I. B. Utama, C.-S. Yang, J. Cain, Y. Qin, Y. Shen, Z. Zheng, K. Watanabe, T. Taniguchi, S. Tongay, A. Zettl, and F. Wang, Identification of spin, valley and moiré quasi-angular momentum of interlayer excitons, Nat. Phys., 15, 1140 (2019).
  * (16) K. Tran, G. Moody, F. Wu, X. Lu, J. Choi, K. Kim, A. Rai, D. A. Sanchez, J. Quan, A. Singh, J. Embley, A. Zepeda, M. Campbell, T. Autry, T. Taniguchi, K. Watanabe, N. Lu, S. K. Banerjee, K. L. Silverman, S. Kim, E. Tutuc, L. Yang, A. H. MacDonald, and Xiaoqin Li, Evidence for moiré excitons in van der Waals heterostructures, Nature, 567, 71 (2019).
  * (17) T. Deilmann, M. Rohlfing, and U. Wurstbauer, Light-matter interaction in van der Waals hetero-structures, J. Phys.: Condens. Matter, 32, 333002 (2020).
  * (18) Y, Shimazaki, I. Schwartz, K. Watanabe, T. Taniguchi, M. Kroner, and A. Imamoğlu, Strongly correlated electrons and hybrid excitons in a moiré heterostructure, Nature, 580, 472 (2020).
  * (19) L. Zhang, Z. Zhang, F. Wu, D. Wang, R. Gogna, S. Hou, K. Watanabe, T. Taniguchi, K. Kulkarni, T. Kuo, S. R. Forrest, and H. Deng, Twist-angle dependence of moiré excitons in WS2/MoSe2 heterobilayers, Nat. Comm., 11, 5888 (2020).
  * (20) Y. Tang, J. Gu, S. Liu, K. Watanabe, K. F. Mak, and J. Shan, Tuning layer-hybridized moiré excitons by the quantum-confined Stark effect, Nat. Nanotech., 16, 52 (2021).
  * (21) Y. Jiang, S. Chen, W. Zheng, B. Zheng, and A. Pan, Interlayer exciton formation, relaxation, and transport in TMD van der Waals heterostructures, Light: Science & Applications, 10, 72 (2021).
  * (22) Y. Liu, C. Zeng, J. Yu, J. Zhong, B. Li, Z. Zhang, Z. Liu, Z. M. Wang, A. Pan, and X. Duan, Moiré superlattices and related moiré excitons in twisted van der Waals heterostructures, Chem. Soc. Rev., 50, 6401 (2021).
  * (23) N. P. Wilson, W. Yao, J. Shan, and X. Xu, Excitons and emergent quantum phenomena in stacked 2D semiconductors, Nature, 599, 7885 (2018).
  * (24) D. Huang, J. Choi, C.-K. Shih, and X. Li, Excitons in semiconductor moiré superlattices, Nat. Nanotech., 17, 227 (2022).
  * (25) B.-H. Lin, Y.-C. Chao, I-T. Hsieh, C.-P. Chuu, C.-J. Lee, F.-H. Chu, L.-S. Lu, W.-T. Hsu, C.-W. Pao, C.-K. Shih, J.-J. Su, and W.-H. Chang, Remarkably Deep Moiré Potential for Intralayer Excitons in MoSe2/MoS2 Twisted Heterobilayers, Nano Lett., (2023).
  * (26) A. Chaves, J. G. Azadani, V. O. Özçelik, R. Grassi, and T. Low, Electrical control of excitons in van der Waals heterostructures with type-II band alignment, Phys. Rev. B, 98, 121302(R) (2018).
  * (27) W.-T. Hsu, B.-H. Lin, L.-S. Lu, M.-H. Lee, M.-W. Chu, L.-J. Li, W. Yao, W.-H. Chang, and C.-K. Shih, Tailoring excitonic states of van der Waals bilayers through stacking configuration band alignment, and valley spin, Sci. Adv., 5, eaax7407 (2019).
  * (28) J. Hagel, S. Brem, C. Linderälv, P. Erhart, and E. Malic, Exciton landscape in van der Waals heterostructures, Phys. Rev. Rearch., 3, 043217 (2021).
  * (29) H. Yu, Y. Wang, Q. Tong, X. Xu, and W. Yao, Anomalous light cones and valley optical selection rules of interlayer excitons in twisted heterobilayers, Phys. Rev. Lett., 115, 187002 (2015).
  * (30) H. Yu, G.-B. Liu, J. Tang, X. Xu, and W. Yao, Moiré excitons: From programmable quantum emitter arrays to spin-orbit–coupled artificial lattices, Science advances, 3, e1701696 (2017).
  * (31) F. Wu, T. Lovorn, and A. H. MacDonald, Topological exciton bands in moiré heterojunctions, Phys. Rev. Lett., 118, 147401 (2017).
  * (32) F. Wu, T. Lovorn, and A. H. MacDonald Theory of optical absorption by interlayer excitons in transition metal dichalcogenide heterobilayers, Phys. Rev. B, 97, 035306 (2018).
  * (33) M. Danovich, D. A. Ruiz-Tijerina, R. J. Hunt, M. Szyniszewski, N. D. Drummond, and V. I. Fal’ko, Localized interlayer complexes in heterobilayer transition metal dichalcogenides, Phys. Rev. B, 97, 195452 (2018).
  * (34) D. A. Ruiz-Tijerina and V. I. Fal’ko, Interlayer hybridization and moiré superlattice minibands for electrons and excitons in heterobilayers of transition-metal dichalcogenides, Phys. Rev. B, 99, 125424 (2019).
  * (35) D. A. Ruiz-Tijerina, I. Soltero, and F. Mireles, Theory of moiré localized excitons in transition metal dichalcogenide heterobilayers, Phys. Rev. B, 102, 195403 (2020).
  * (36) S. Brem, K.-Q. Lin, R. Gillen, J. M. Bauer, J. Maultzsch, J. M. Luptonb, and E. Malic, Hybridized intervalley moiré excitons and flat bands in twisted WSe 2 bilayers, Nanoscale, 12, 11088 (2020).
  * (37) S. Brem, C. Linderälv, and P. Erhart, and E. Malic, Tunable phases of moiré excitons in van der Waals heterostructures, Nano Lett., 20, 8534 (2020).
  * (38) H. Guo, X. Zhang, and G. Lu, Shedding light on moiré excitons: A first-principles perspective, Sci. Adv., 6, eabc5638 (2020).
  * (39) J. Choi, M. Florian, A. Steinhoff, D. Erben, K. Tran, D. S. Kim, L. Sun, J. Quan, R. Claassen, S. Majumder, J. A. Hollingsworth, T. Taniguchi, K. Watanabe, K. Ueno, A. Singh, G. Moody, F. Jahnke, and X. Li, Twist angle-dependent interlayer exciton lifetimes in van der Waals heterostructures, Phys. Rev. Lett., 126, 047401 (2021).
  * (40) A. Hichri, T. Amand, and S. Jaziri, Resonance energy transfer from moiré-trapped excitons in MoSe2/WSe2 heterobilayers to graphene: Dielectric environment effect, Phys. Rev. Mater., 5, 114002 (2021).
  * (41) H. Yu and W. Yao, Luminescence anomaly of dipolar valley excitons in homobilayer semiconductor moiré superlattices, Phy. Rev. X, 11, 021042 (2021).
  * (42) A. Julku, Nonlocal interactions and supersolidity of moiré excitons, Phys. Rev. B, 106, 035406 (2022).
  * (43) M. H. Naik, E. C. Regan, Z. Zhang, Y.-H. Chan, Z. Li, D. Wang, Y. Yoon, C. S. Ong, W. Zhao, S. Zhao, M. I. B. Utama, B. Gao, X. Wei, M. Sayyad, K. Yumigeta, K. Watanabe, T. Taniguchi, S. Tongay, F. H. da Jornada, F. Wang, and S. G. Louie, Intralayer charge-transfer moiré excitons in van der Waals superlattices, Nature, 609, 52 (2022).
  * (44) J. Hagel, S. Brem, and E. Malic, Electrical tuning of moiré excitons in MoSe2 bilayers, 2D Mater., 10, 014013 (2023).
  * (45) F. Wu, T. Lovorn, E. Tutuc, and A. H. MacDonald, Hubbard model physics in transition metal dichalcogenide moiré bands, Phys. Rev. Lett., 121, 026402 (2018).
  * (46) N. S. Rytova, Vestn. Mosk. Univ. Fiz. Astron., 3, 30 (1967).
  * (47) L. V. Keldysh, J. Exp. Theoret. Phys. Lett., 29, 658 (1979).
  * (48) K. S. Thygesen, Calculating excitons, plasmons, and quasiparticles in 2D materials and van der Waals heterostructures, 2D Mater., 4, 022004 (2017).
  * (49) L. S. R. Cavalcante, A. Chaves, B. Van Duppen, F. M. Peeters, and D. R. Reichman, Electrostatics of electron-hole interactions in van der Waals heterostructures, Phys. Rev. B, 97, 125427 (2018).
  * (50) M. Van der Donck and F. M. Peeters, Interlayer excitons in transition metal dichalcogenide heterostructures, Phys. Rev. B, 98, 115104 (2018).
  * (51) F. García Flórez, Laurens D. A. Siebbeles, and H. T. C. Stoof, Effects of material thickness and surrounding dielectric medium on Coulomb interactions and two-dimensional excitons, Phys. Rev. B, 102, 125303 (2020).
  * (52) H. E. Hannachi, D. Elmaghraoui, and S. Jaziri, Moiré interlayer exciton relative and center of mass motions coupling. Effect on 1s-np interlayer exciton THz transitions, The European Physical Journal Plus, 138, 1.
  * (53) J.-Z. Zhang and J.-Z. Ma, Two-dimensional excitons in monolayer transition metal dichalcogenides from radial equation and variational calculations, J. Phys.: Condens. Matter, 31, 105702 (2019).
  * (54) S. Wu, L. Cheng, and Q. Wang, Exciton states and absorption spectra in freestanding monolayer transition metal dichalcogenides: A variationally optimized diagonalization method, Phys. Rev. B, 100, 115430 (2019).
  * (55) Y.-W. Chang and Y.-C. Chang, Variationally optimized orbital approach to trions in two-dimensional materials, J. Chem. Phys., 155, 024110 (2021).
  * (56) Y.-W. Chang and Y.-C. Chang, Theory of magnetic-field effect on trions in two-dimensional materials, J. Chem. Phys., 157, 044104 (2022).
  * (57) M. Goryca, J. Li, A. V. Stier, T. Taniguchi, K. Watanabe, E. Courtade, S. Shree, C. Robert, B. Urbaszek, X. Marie and S. A. Crooker, Revealing exciton masses and dielectric properties of monolayer semiconductors with high magnetic fields, Nat. Comm., 10, 4172 (2019).
  * (58) J. A. Wilson, A. D. and Yoffe, The transition metal dichalcogenides discussion and interpretation of the observed optical, electrical and structural properties, Advances in Physics, 18, 193 (1969).
  * (59) C. Zhang, C. Gong, Y. Nie, K.-A. Min, C. Liang, Y. J. Oh, H. Zhang, W. Wang, S. Hong, L. Colombo, R. M. Wallace, and K. Cho, Systematic study of electronic structure and band alignment of monolayer transition metal dichalcogenides in Van der Waals heterostructures, 2D Mater., 4, 015026 (2016).
  * (60) A. Kormányos, G. Burkard, M. Gmitra, J. Fabian, V. Zólyomi, N. D. Drummond, and V. Fal’ko, k· p theory for two-dimensional transition metal dichalcogenide semiconductors, 2D Mater., 2, 022001 (2015).
  * (61) H. Rydberg, M. Dion, N. Jacobson, E. Schröder, P. Hyldgaard, S. I. Simak, D. C. Langreth, and B. I. Lundqvist, Van der Waals density functional for layered structures, Phys. Rev. Lett., 91, 126402 (2003).
  * (62) S. Shabani, D. Halbertal, W. Wu, M. Chen, S. Liu, J. Hone, W. Yao, D. N. Basov, X. Zhu, and A. N. Pasupathy, Deep moiré potentials in twisted transition metal dichalcogenide bilayers, Nat. Phys., 17, 720 (2021).
  * (63) J. Schwinger, L. L. DeRaad Jr, K. Milton, and W.-Y. Tsai, Classical electrodynamics, Westview Press (1998).
  * (64) H. Zheng, D. Zhai, and W. Yao, Twist versus heterostrain control of optical properties of moiré exciton minibands, 2D Mater., 8, 044016 (2021).
  * (65) H. Zheng, D. Zhai, and W. Yao, Anomalous magneto-optical response and chiral interface of dipolar excitons at twisted valleys, Nano Letters, 22, 5466 (2022).