# The strong modulation limit of excitons and trions in moiré materials

Yongxin Zeng and Allan H. MacDonald  Department of Physics, University of
Texas at Austin, Austin, TX 78712

###### Abstract

The optical properties of weakly-doped two-dimensional materials are dominated
by strong exciton and trion absorption and luminescence features. In this
article we examine the influence of moiré patterns in semiconductor
heterobilayers on exciton and trion states in the limit of strong moiré
modulation potentials, commenting on similarities and differences compared to
the case of excitons and trions in semiconductor quantum dots. We discuss
strategies for using optical properties as quantitative probes of moiré
materials, and the prospects for exploiting moiré materials to design unique
light emitters.

##  I Introduction

Moiré superlattices form when two-dimensional (2D) crystals are stacked with a
small twist angle or lattice mismatch. When the host materials are
semiconductors or semimetals, moiré superlattices are accurately described by
continuum models with the periodicity of the moiré pattern, giving rise to
artificial two-dimensional moiré materials with controllable large lattice
constants. Extensive theoretical and experimental work on twisted bilayer
graphene [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
20, 21, 22] and twisted bilayer transition metal dichalcogenides (TMDs) [23,
24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37] over the past several
years has established moiré materials as an important platform for studies of
strongly correlated electron physics.

In addition to modifying low-energy electronic properties, the moiré
potentials in twisted bilayer TMDs also significantly modify interband
excitonic collective modes, and associated optical properties [38, 39, 40, 41,
42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53]. The moiré modulation splits
both intralayer and interlayer excitonic peaks in the optical spectrum into a
series of subpeaks and modifies optical selection rules. Since the electrons,
holes, and excitons are often strongly localized near minima of the periodic
moiré potential, moiré superlattice systems can act as a new and more
perfectly periodic realization of quantum dot arrays [39, 40], and be used to
engineer single-photon emitters and entangled photon sources [48, 43] with
potential applications in quantum optics.

Theories of the excitonic properties of moiré materials have progressed by
treating excitons as point particles moving in a periodic moiré potential [38,
39, 40]. Recent numerical and experimental studies show that moiré modulation
strengthens when strain is taking into account [54, 55, 56, 57], and that the
conduction band electrons and valence band holes that combine to form excitons
in a moiré material can be attracted to different lateral positions within the
moiré pattern [58, 59, 60, 61, 62], distorting the excitonic modes. In the
strong modulation limit, the electron and hole moiré band widths become
extremely small and the moiré material can be viewed as an accurately periodic
array of quantum dots. In this article we study how strong confinement and
lateral shifts between the electron and hole confining centers affect the
properties of excitons and trions in the quantum dot array limit of a moiré
material. In Sec. II we introduce our model for electrons and holes in a moiré
superlattice, and present results for Coulomb interaction matrix elements in
the single-particle harmonic oscillator basis that we use for numerical
calculations. In Secs. III and IV we calculate the ground state energies and
optical absorption strengths of excitons and trions respectively, by exact
diagonalization (ED) of the corresponding few-body Hamiltonians in the
harmonic oscillator basis. Finally, in Sec. V we discuss the implications of
our results, and point to future work that can be approached with our method.
We find that the optical properties of systems in which electrons and holes
are trapped at the same position within the periodic moiré pattern are very
different from those of systems in which electrons and holes are trapped at
different positions; the former case is distinguished by fewer but stronger
absorption features because of the approximate applicability of angular
momentum conservation selection rules. Magnetic fields are potentially useful
in characterizing moiré materials, but because of of the strong electron-hole
binding in two-dimensional material excitons and trions will yield substantial
corrections that can be used to identify elecron, hole, exciton, and trion
configurations only at very strong magnetic fields.

##  II The model

Figure 1: (a) A hole-doped moiré superlattice at filling factor $\nu=1/3$. The
blue and green circles represent the electron and hole confining centers
respectively, which are assumed here to be located at different high-symmetry
stacking sites. The filled and empty circles distinguish sites that are
occupied by a quasiparticle from those that are empty. The doped holes form a
Wigner crystal in the ground state. (b) The moiré lattice with some electron-
hole excitations. The yellow ellipses show three different types of electron-
hole excitations: an exciton (right), a trion (lower left), and an example of
a more widely separated electron-hole pair (top left). (c) If the center of an
electron site and its neighboring hole site is associated with an exciton or
trion excitation, as marked by the red crosses and triangles, together they
form a kagome lattice. (d) Same as (b), except that the electrons and holes
are localized at the same sites. The yellow circles show two excitons and a
trion excitation.

In a TMD heterobilayer moiré superlattice, electrons and holes tend to
localize at one of three high-symmetry sites in a moiré unit cell [54, 55, 56,
57, 58, 59, 60, 61, 62] at which the local stacking places the metal atoms of
one layer above either metal, chalcogen, or empty sites of the layer below. As
summarized in Fig. 1, there is an important distinction between the case in
which electrons and hole are localized at the same positions in the moiré
pattern or at different positions. It is not yet confidently known which among
the various TMD heterobilayer systems fall into the former category and which
fall into the latter category. When the electron and hole sites are the same,
they form a triangular lattice of sites at which excitons can be localized.
When the electron and hole sites are separated in the 2D plane, however, their
triangular lattices combine to form a hexagonal lattice like the one shown in
Fig. 1(a). In a hole-doped system some of the hole sites are occupied, and at
certain fractional fillings the ground state of the doped holes is a lattice
Wigner crystal [24, 25, 26, 36]. Fig. 1 shows an example at filling factor
$\nu=1/3$. Possible low-energy excitations of the system include electron-hole
pairs at empty neighboring sites (an exciton), or an electron-hole pair at a
filled hole site and its neighboring electron site (a trion). Excitations with
electrons and holes farther away are also possible, but have higher energies
and weaker optical absorptions. If we ignore these, and consider only near-
neighbor exciton and trion excitations, each bond between an electron site and
its neighboring hole site can be associated with a possible exciton or trion
excitation (marked by red crosses and triangles in Fig. 1(c)), depending on
whether the hole site is occupied or not. Together these excitations form a
kagome lattice. Excitons and trions on different sites are coupled by
intersite hopping of their constituents, and via Coulomb interactions. We do
not account for these complications here.

Figure 2: Schematic illustration of the model system – a few electrons and
holes confined in a 2D quantum dot. The black dots stand for the confining
potential minima for electrons and holes, and the blue and green dots
represent electrons and holes respectively. The vectors $\bm{r}_{ei}$ and
$\bm{r}_{hj}$ are defined relative to the potential minima for electrons and
holes respectively, which are separated by the vector $\bm{s}=s\hat{x}$ when
the electrons and holes are in the same layer and by $s\hat{x}+d\hat{z}$ when
the electrons and holes are in layers separated by distance $d$.

In this article we consider both two-particle exciton and three-particle trion
electron-hole complexes, ignoring the possible role of coupling to other
excitations of the system. We study a system of $N_{e}$ electrons and $N_{h}$
holes confined in a 2D quantum dot (Fig. 2) with harmonic confining potentials
for the electrons and holes, and a magnetic field $\bm{B}=B\hat{z}$ applied
perpendicular to the 2D plane. We consider the general case in which the
electrons and holes have different effective masses $m_{e}$ and $m_{h}$ and
are localized at different lateral positions in the same or different layers
by confining potentials that can have different oscillator frequencies
$\omega_{e}$ and $\omega_{h}$. For heterobilayer moiré superlattices, states
with electrons and holes in the same layer are relevant in absorption or
reflection measurements and states with electrons and holes in different
layers are relevant in luminescence experiments. The system is described by
Hamiltonian

$H=H_{e}+H_{h}+V_{ee}+V_{hh}+V_{eh}.$ (1)

The electron and hole single-particle terms are

$\displaystyle H_{e}$
$\displaystyle=\sum_{i=1}^{N_{e}}\left[\frac{(\bm{p}_{ei}+e\bm{A}_{e}(\bm{r}_{ei}))^{2}}{2m_{e}}+\frac{1}{2}m_{e}\omega_{e}^{2}r_{ei}^{2}\right],$
(2) $\displaystyle H_{h}$
$\displaystyle=\sum_{j=1}^{N_{h}}\left[\frac{(\bm{p}_{hj}-e\bm{A}_{h}(\bm{r}_{hj}))^{2}}{2m_{h}}+\frac{1}{2}m_{h}\omega_{h}^{2}r_{hj}^{2}\right],$
(3)

where $i$ and $j$ are labels for electrons and holes, and the origins of
$\bm{r}_{e}$ and $\bm{r}_{h}$ are chosen to be at the minima of the
corresponding confining potentials. We choose a gauge such that
$\bm{A}_{e}(\bm{r}_{e})=\frac{1}{2}\bm{B}\times\bm{r}_{e}$ and
$\bm{A}_{h}(\bm{r}_{h})=\frac{1}{2}\bm{B}\times\bm{r}_{h}$, i.e., the
symmetric gauge, but with different origins for electrons and holes. The
interaction terms are

$\displaystyle V_{ee}$
$\displaystyle=\sum_{i_{1}<i_{2}}^{N_{e}}\frac{e^{2}}{\epsilon|\bm{r}_{ei_{1}}-\bm{r}_{ei_{2}}|},$
(4) $\displaystyle V_{hh}$
$\displaystyle=\sum_{j_{1}<j_{2}}^{N_{h}}\frac{e^{2}}{\epsilon|\bm{r}_{hj_{1}}-\bm{r}_{hj_{2}}|},$
(5) $\displaystyle V_{eh}$
$\displaystyle=\sum_{i=1}^{N_{e}}\sum_{j=1}^{N_{h}}\frac{-e^{2}}{\epsilon\sqrt{|\bm{r}_{ei}-\bm{r}_{hj}+\bm{s}|^{2}+d^{2}}},$
(6)

where $\bm{s}$ is the displacement between the confining potential minimum of
electrons and holes. We choose the coordinate axes such that the 2D crystals
are in the $\hat{x}$-$\hat{y}$ plane and the in-plane component of $\bm{s}$ is
along the $\hat{x}$ direction (see Fig. 2). $d$ is the distance between the
electron and hole layers ($d=0$ when the electrons and holes are in the same
layer), and $\epsilon$ is the dielectric constant of the dielectric
surrounding the system. The most commonly employed surrounding dielectric is
hexagonal boron nitride (hBN), which has an anisotropic dielectric tensor. In
this case $\epsilon=\sqrt{\epsilon_{zz}\epsilon_{\perp}}$ and the effective
layer separation $d$ is related to the physical separation $d_{0}$ by
$d=d_{0}\sqrt{\epsilon_{\perp}/\epsilon_{zz}}$, where $\epsilon_{\perp}$ and
$\epsilon_{zz}$ are respectively the perpendicular-to-plane and in-plane
components of the dielectric tensor.

We diagonalize the few-particle Hamiltonian using a configuration interaction
approach with basis state Slater determinants constructed from the eigenstates
of $H_{e}$ and $H_{h}$. Our calculations are similar to those reported in
earlier papers [63, 64, 65, 66] on quantum dots in semiconductor quantum
wells, except that we consider here the possibility that the electrons and
holes are localized at different lateral positions. The harmonic oscillator
eigenstates $|n_{+}n_{-}\rangle$ are labelled by two quantum numbers $n_{+}$
and $n_{-}$, with energy eigenvalues

$E_{n_{+}n_{-}}=(n_{+}+\frac{1}{2})\hbar\Omega_{+}+(n_{-}+\frac{1}{2})\hbar\Omega_{-},$
(7)

where $\Omega_{\pm}=\Omega\pm\omega_{c}/2$ with
$\Omega=\sqrt{\omega^{2}+(\omega_{c}/2)^{2}}$ and $\omega_{c}=eB/m$. (See
Appendix A for derivations.) To diagonalize the Hamiltonian we choose basis
states that are antisymmetrized combination of electron-hole many particle
product states

$|n_{e1+}n_{e1-},n_{e2+}n_{e2-},\dots;n_{h1+}n_{h1-},n_{h2+}n_{h2-},\dots\rangle.$
(8)

Antisymmetrization is not required between electrons or holes since these can
be considered to be distinguishable. The two-particle matrix elements for
Coulomb interactions are evaluated by expanding the interactions in terms of
ladder operators of the single-particle basis. For interactions among
electrons or holes, the result is

$\begin{split}\langle
n_{1+}^{\prime}n_{1-}^{\prime},n_{2+}^{\prime}n_{2-}^{\prime}|V_{\alpha\alpha}|n_{1+}n_{1-},n_{2+}n_{2-}\rangle&=\frac{\sqrt{\beta_{\alpha}}}{2}\frac{e^{2}}{\epsilon
l_{B}}\left(\prod_{i\lambda}n_{i\lambda}^{\prime}!n_{i\lambda}!\right)^{1/2}(-1)^{n_{1+}^{\prime}+n_{1-}+n_{2+}^{\prime}+n_{2-}}\delta_{l^{\prime}l}\\\
&\times\sum_{k_{i\lambda}=0}^{\min(n_{i\lambda}^{\prime},n_{i\lambda})}\left(\prod_{i\lambda}\frac{(-1)^{k_{i\lambda}}2^{-p}}{k_{i\lambda}!(n_{i\lambda}^{\prime}-k_{i\lambda})!(n_{i\lambda}-k_{i\lambda})!}\right)\cdot\Gamma(p+\frac{1}{2}),\end{split}$
(9)

where $\alpha=e,h$, the sums define separate factors for $\lambda=\pm$ and
particle $i=1,2$, the conserved angular momentum quantum numbers
$l=n_{1+}-n_{1-}+n_{2+}-n_{2-}$ and
$l^{\prime}=n_{1+}^{\prime}-n_{1-}^{\prime}+n_{2+}^{\prime}-n_{2-}^{\prime}$,
and $p=\sum_{i\lambda}(n_{i\lambda}^{\prime}+n_{i\lambda}-2k_{i\lambda})/2$.
In Eq. (9), the magnetic length $l_{B}=\sqrt{\hbar/eB}$ and we have defined
$\beta_{\alpha}=\sqrt{1+(2\omega_{\alpha}/\omega_{\alpha c})^{2}}$.

The two-particle matrix elements of $V_{eh}$ are similar but do not conserve
angular momentum because the electron and hole oscillators have different
centers in general:

$\begin{split}\langle
n_{e+}^{\prime}n_{e-}^{\prime},n_{h+}^{\prime}n_{h-}^{\prime}|V_{eh}|n_{e+}n_{e-},n_{h+}n_{h-}\rangle&=-\sqrt{\frac{\beta_{eh}}{2}}\frac{e^{2}}{\epsilon
l_{B}}\left(\prod_{\alpha\lambda}n_{\alpha\lambda}^{\prime}!n_{\alpha\lambda}!\right)^{1/2}(-1)^{n_{e+}^{\prime}+n_{e-}+n_{h+}^{\prime}+n_{h-}}\\\
&\times\sum_{k_{\alpha\lambda}=0}^{\min(n_{\alpha\lambda}^{\prime},n_{\alpha\lambda})}\left(\prod_{\alpha\lambda}\frac{(-1)^{k_{\alpha\lambda}}(\frac{\beta_{eh}}{\beta_{\alpha}})^{\frac{1}{2}(n_{\alpha\lambda}^{\prime}+n_{\alpha\lambda})-k_{\alpha\lambda}}}{k_{\alpha\lambda}!(n_{\alpha\lambda}^{\prime}-k_{\alpha\lambda})!(n_{\alpha\lambda}-k_{\alpha\lambda})!}\right)\cdot
I,\end{split}$ (10)

where $\beta_{eh}=\beta_{e}\beta_{h}/(\beta_{e}+\beta_{h})$, and the
$(n^{\prime},n,k)$-dependent integral

$I=\int_{0}^{\infty}d\tilde{q}^{2}\,e^{-\tilde{q}^{2}}(\tilde{q}^{2})^{p-\frac{1}{2}}e^{-\tilde{q}\tilde{d}}\cdot\int_{0}^{2\pi}\frac{d\phi}{2\pi}\,(ie^{i\phi})^{l-l^{\prime}}e^{i\tilde{q}\tilde{s}\cos\phi},$
(11)

where we introduced the dimensionless lengths
$\tilde{d}=\frac{\sqrt{2\beta_{eh}}}{l_{B}}d$,
$\tilde{s}=\frac{\sqrt{2\beta_{eh}}}{l_{B}}s$, the total angular momentum
quantum numbers $l=n_{e+}-n_{e-}-n_{h+}+n_{h-}$,
$l^{\prime}=n_{e+}^{\prime}-n_{e-}^{\prime}-n_{h+}^{\prime}+n_{h-}^{\prime}$,
and
$p=\sum_{\alpha\lambda}(n_{\alpha\lambda}^{\prime}+n_{\alpha\lambda}-2k_{\alpha\lambda})/2$.
The radial part of the integral in Eq. (11) can be computed analytically and
we obtain

$I=4^{-p}\,\Gamma(2p+1)\int_{0}^{2\pi}\frac{d\phi}{2\pi}\,(ie^{i\phi})^{l-l^{\prime}}\,U(p+\frac{1}{2},\frac{1}{2},\frac{1}{4}(\tilde{d}-i\tilde{s}\cos\phi)^{2}),$
(12)

where $U$ is the confluent hypergeometric function of the second kind. The
angular integral must be evaluated numerically for $s\neq 0$. For $s=0$,

$I(s=0)=\delta_{l^{\prime}l}\cdot
4^{-p}\,\Gamma(2p+1)\,U(p+\frac{1}{2},\frac{1}{2},\frac{\tilde{d}^{2}}{4}).$
(13)

When $d\to 0$, $I(s=0)\to\delta_{l^{\prime}l}\cdot\Gamma(p+\frac{1}{2})$. More
detail on the derivations of these formulas is provided in Appendix A. The
generalization of these formulas to properly antisymmetrized multi-particle
states basis with spin and/or valley degrees of freedom is straightforward.

##  III Excitons

We first consider the case of excitons – two-particle states with an electron
and a hole that are distinguishable. We neglect intervalley exchange
interactions [67, 68, 69, 70] so that the excitonic states are valley
independent. The $i$th eigenstate of the Hamiltonian can be written as

$|\Psi_{eh}^{(i)}\rangle=\sum_{n_{e}n_{h}}C_{n_{e}n_{h}}^{(i)}|n_{e};n_{h}\rangle,$
(14)

where $n_{e}$ is a shorthand for the two indices $(n_{e+},n_{e-})$ that
specify which electron basis state is occupied, and similarly for $n_{h}$. The
coefficients $C_{n_{e}n_{h}}^{(i)}$ are obtained by numerical diagonalization
of the Hamiltonian in the magnetic harmonic oscillator basis
$|n_{e};n_{h}\rangle$. In our model, exciton energies are given relative to
the band gap between potential extrema.

In the weak-confinement limit, the relative-motion electron-hole pair state is
very weakly distorted, and the electron-hole center-of-mass moves in a weak
modulation potential. The characteristic length scale for bound electron-hole
pairs is the effective Bohr radius $a_{B}^{*}=\epsilon\hbar^{2}/me^{2}$ where
$m=m_{e}m_{h}/(m_{e}+m_{h})$ is the reduced mass of the electron-hole pair.
The characteristic energy scale is the effective Rydberg ${\rm
Ry}^{*}=e^{2}/2\epsilon a_{B}^{*}=me^{4}/2\epsilon^{2}\hbar^{2}$. In our
calculations we set $m_{e}=0.4m_{e0}$, $m_{h}=0.5m_{e0}$ where $m_{e0}$ is the
free electron mass, and $\epsilon=5$, as approximate parameters for intralayer
excitons in bilayer TMDs [71] surrounded by hBN. Then we find
$a_{B}^{*}=$1.19\text{\,}\mathrm{n}\mathrm{m}$$ and ${\rm
Ry^{*}}=$121\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$.

In the strong-confinement and large-$s$ limit, the electron and hole are each
confined in a harmonic potential, with length scale given by the oscillator
length $l_{\alpha}=\sqrt{\hbar/m_{\alpha}\omega_{\alpha}}$ and energy scale
$\hbar\omega_{\alpha}$. The magnetic field provides another length scale, the
magnetic length $l_{B}=\sqrt{\hbar/eB}$, and another energy scale, the
cyclotron energy $\hbar\omega_{\alpha c}=\hbar eB/m_{\alpha}$. For our
explicit calculations we choose
$\hbar\omega_{e}=\hbar\omega_{h}=$50\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$.
Then the harmonic confinement lengths
$l_{e}=$1.95\text{\,}\mathrm{n}\mathrm{m}$$ and
$l_{h}=$1.75\text{\,}\mathrm{n}\mathrm{m}$$. At magnetic field
$B=$10\text{\,}\mathrm{T}$$, the magnetic length
$l_{B}=$8.11\text{\,}\mathrm{n}\mathrm{m}$$ is considerably longer and the
cyclotron energies
$\hbar\omega_{ec}=$2.9\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$,
$\hbar\omega_{hc}=$2.3\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$ considerably
smaller than other length and energy scales. In a magnetic field, the
oscillator length $L_{\alpha}=\sqrt{\hbar/m_{\alpha}\Omega_{\alpha}}$ with
$\Omega_{\alpha}=\sqrt{\omega_{\alpha}^{2}+(\omega_{\alpha c}/2)^{2}}$. We
therefore see that due to the large effective masses in the TMDs, the
oscillator length is weakly perturbed by a $B=$10\text{\,}\mathrm{T}$$
magnetic field ($L_{\alpha}\approx l_{\alpha}\ll l_{B}$). The magnetic length
and the oscillator length start to become comparable only for magnetic fields
$\sim$100\text{\,}\mathrm{T}$$; below this field scale exciton states are
weakly perturbed by magnetic fields.

Figure 3: Ground-state exciton energy $E_{eh}^{(0)}$ as a function of the in-
plane separation between the electron and hole confining potential minima $s$,
for intralayer (blue $d=0$) and interlayer (green and yellow
$d=1,$2\text{\,}\mathrm{n}\mathrm{m}$)$ excitons. These results were
calculated at magnetic field $B=$10\text{\,}\mathrm{T}$$, but differ very
little from results calculated at $B=0$. The dots and solid lines are
numerical results from ED calculations, and the dotted curves are the
large-$s$ approximation of (15). The blue cross shows the weak-confinement
approximation (16) expression which adds a zero-point confinement energy to
the free exciton energy.

Fig. 3 shows numerical results for the ground-state exciton energy at
$B=$10\text{\,}\mathrm{T}$$ as a function of the lateral electron-hole
separation $s$. When $s\gg L_{e},L_{h}$, the electron and hole are confined in
widely separated harmonic potentials, electron-hole correlations are
negligible, and interactions can be treated at lowest order in perturbation
theory. In this limit, the ground state energy is

$E_{eh}^{(0)}(s)=E_{e}^{(0)}+E_{h}^{(0)}-\frac{e^{2}}{\epsilon\sqrt{s^{2}+d^{2}}},$
(15)

where $E_{e}^{(0)}=\hbar(\Omega_{e+}+\Omega_{e-})/2=\hbar\Omega_{e}$ is the
ground state of an electron in the quantum dot (and similarly for holes) with
$\Omega_{\pm}$ and $\Omega$ defined after Eq. (7). The dotted curves in Fig. 3
show that the approximate ground state energy calculated from Eq. (15), which
entirely neglects electron-hole correlations, agrees well with the numerical
results when $s\gg L_{e},L_{h}$. For intermediate
$s\sim$3\text{\,}\mathrm{n}\mathrm{m}$$, Coulomb attraction distorts the
wavefunctions of the electron and hole, which effectively reduces the
electron-hole separation and results in a lower ground-state energy than that
given by the large-$s$ approximation (15). For small $s$, the electron and
hole wavefunctions overlap in the 2D plane and the approximation of Eq. (15),
which becomes singular in this limit when $d=0$, breaks down. The singular
behavior at small $s$ is cut off by the kinetic-energy pressure required by
the Heisenberg uncertainly relation, and this results in a ground state energy
that is higher than that given by Eq. (15). The approximate ground state
energy of Eq. (15) is in reasonable agreement with numerical results when
$s\gtrsim L_{e},L_{h}$.

We see from Fig. 3 that exciton properties depend qualitatively on the ratio
of $s$ to confinement oscillator lengths. As indicated in Fig. 1, the
experimentally relevant value of $s$ is either $s=0$, in which case electron-
hole correlation is always important, or $s=a_{M}/\sqrt{3}$, where $a_{M}$ is
the lattice constant of the triangular moiré Bravais lattice. In the former
case (see below) the large-$s$ limit is never relevant. In the latter case,
the large-$s$ limit is always achieved at large enough moiré period since the
oscillator length of electrons or holes in moiré materials varies [23, 37]
approximately as the square root of the moiré period. The large-$s$ limit is
simple not only for excitons, but for any finite-particle-number electron-hole
complex, since the interactions between electrons and holes $V_{eh}$ (Eq. (6))
simply add an electrostatic contribution to energies and do not alter
electronic wavefunctions.

If confinement is weak, the $s=0$ system can be viewed as a tightly bound
exciton of mass $M=m_{e}+m_{h}$ moving in a weak confining potential with
angular frequency
$\omega=\sqrt{(m_{e}\omega_{e}^{2}+m_{h}\omega_{h}^{2})/(m_{e}+m_{h})}$ [64].
The ground state energy is then the sum of the ground state energies of the
harmonic oscillator and the tightly bound exciton. For the intralayer exciton
case $d=0$, the system is similar to a 2D hydrogen atom, and the total energy
is given by the analytic result:

$E_{eh}^{(0)}(0)=\hbar\omega-4\,{\rm Ry^{*}}.$ (16)

This approximate result is indicated by the blue cross in Fig. 3, where we see
that the parameters chosen for this calculation are not in the weak
confinement limit – the exciton radius $\sim a_{B}^{*}$ is not small compared
with the confining length of the harmonic potential $L=\sqrt{\hbar/M\omega}$ –
and Eq. (16) is therefore inaccurate. The stronger confinement weakens the
correlations between electron and hole positions and therefore increases the
total energy.

Figure 4: Diamagnetic exciton energy shift $\Delta
E_{eh}^{(0)}(B)=E_{eh}^{(0)}(B)-E_{eh}^{(0)}(B=0)$ of the excitonic ground
state in a magnetic field $B$, for interlayer excitons at
$d=1,$2\text{\,}\mathrm{n}\mathrm{m}$$ and
$s=0,$3\text{\,}\mathrm{n}\mathrm{m}$$. The circular and triangular dots
connected by solid or dashed lines are numerical results from ED calculations,
and the dotted curve shows the energy increase obtained from the $s\to\infty$
limit (15) expression in which only the zero-point confinement energy is
shifted.

Fig. 4 shows how the exciton ground state energy shifts in a magnetic field.
In the large-$s$ limit it follows from Eq. (15) that the diamagnetic exciton
energy shift comes from the single-particle energies of electrons and holes:

$\begin{split}&\Delta E_{eh}^{(0)}(B)\equiv
E_{eh}^{(0)}(B)-E_{eh}^{(0)}(B=0)\\\
&\xrightarrow{s\to\infty}\hbar(\Omega_{e}+\Omega_{h}-\omega_{e}-\omega_{h})\\\
&\approx\frac{e^{2}B^{2}}{8}\left(\frac{l_{e}^{2}}{m_{e}}+\frac{l_{h}^{2}}{m_{h}}\right).\end{split}$
(17)

The last line of the above equation applies for a small magnetic field $B\ll
2m_{\alpha}\omega_{\alpha}/e$. Fig. 4 shows that for
$s=$3\text{\,}\mathrm{n}\mathrm{m}$$ the numerical results are already close
to the $s\to\infty$ limit. The diamagnetic exciton energy shift grows
quadratically with magnetic field both in the strong-binding case ($s=0$) and
in the strong-confinement large-$s$ limit, as expected for the ground state
with zero angular momentum. By first-order perturbation theory we see from
Eqs. (2) and (3) that the diamagnetic exciton energy shift for a zero-angular-
mementum state is

$\Delta
E_{eh}^{(0)}(B)\approx\frac{e^{2}B^{2}}{8}\left(\frac{\langle\bm{r}_{e}^{2}\rangle}{m_{e}}+\frac{\langle\bm{r}_{h}^{2}\rangle}{m_{h}}\right),$
(18)

where the angle brackets denote ground state expectation values. When $s=0$
the electron and hole orbitals shrink due to strong Coulomb attraction,
resulting in a weaker diamagnetic response. The diamagnetic response of
excitons can therefore be used as an experimental probe of the binding
strength of the electrons and holes and therefore distinguish the case in
which they localized at the same site from the case in which they are
localized at different sites.

Next we study the optical absorption of the intralayer exciton states ($d=0$).
In our model the current operator is given up to a constant factor by

$\begin{split}\hat{j}\sim\sum_{\tau}\int
d\bm{r}\,&[e_{\tau}^{\dagger}(\bm{r})h_{\tau}^{\dagger}(\bm{r})\exp(i\frac{eBs}{2\hbar}y)\\\
&+h_{\tau}(\bm{r})e_{\tau}(\bm{r})\exp(-i\frac{eBs}{2\hbar}y)],\end{split}$
(19)

where $e^{\dagger}$ and $h^{\dagger}$ are electron and hole creation
operators, $\tau=\pm$ is the valley index, and the phase factors $\exp(\pm
i\frac{eBs}{2\hbar}y)$ come from our gauge choice (see Appendix B for
details). The matrix element for a transition from the vacuum state
$|\Psi_{0}\rangle$ with no electrons or holes to an exciton state
$|\Psi_{eh}^{(i)}\rangle$ (Eq. (14)) is proportional to

$\langle\Psi_{0}|\hat{j}|\Psi_{eh}^{(i)}\rangle\sim\sum_{n_{e}n_{h}}C_{n_{e}n_{h}}^{(i)}\int
d\bm{r}\,\psi_{n_{e}}^{(e)}(\bm{r})\psi_{n_{h}}^{(h)}(\bm{r})e^{-i\frac{eBs}{2\hbar}y},$
(20)

where $\psi_{n_{e}}^{(e)}(\bm{r})$ is the wavefunction of an electron in state
$n_{e}$, and similarly $\psi_{n_{h}}^{(h)}(\bm{r})$ for holes. Then the
absorption rate of state $i$ is proportional to the dimensionless quantity

$A_{eh}^{(i)}=\Big{|}\sum_{n_{e}n_{h}}C_{n_{e}n_{h}}^{(i)}\int
d\bm{r}\,\psi_{n_{e}}^{(e)}(\bm{r})\psi_{n_{h}}^{(h)}(\bm{r})e^{-i\frac{eBs}{2\hbar}y}\Big{|}^{2},$
(21)

which is, roughly speaking, the overlap of the electron and hole
wavefunctions. The optical absorption strength of the ground state is
calculated by numerical integration of the wavefunctions and shown in Fig. 5.
The optical spectrum is a set of $\delta$-function peaks located at the energy
eigenvalues of exciton states and with heights proportional to the
corresponding optical absorption rates:

$A_{eh}(\omega)=\sum_{i}A_{eh}^{(i)}\,\delta(\hbar\omega-E_{eh}^{(i)}-E_{g}),$
(22)

where $E_{g}$ is the single-particle gap between the conduction and valence
bands.

Figure 5: Optical absorption strength $A_{eh}^{(0)}$ of intralayer excitons as
a function of the separation between the electron and hole confining potential
minima $s$. These results were calculated at magnetic field
$B=$10\text{\,}\mathrm{T}$$, but differ very little from results calculated at
$B=0$. The blue dots are the numerical results from ED calculations and the
dotted curve is from the large-$s$ approximation (24). Since the likely $s$ is
either $0$ or $a_{M}/\sqrt{3}$, based on symmetry considerations, and typical
moiré periods are longer than $5\text{\,}\mathrm{n}\mathrm{m}$, a substantial
reduction in oscillator strength is expected in the $s\neq 0$ case.

We can obtain analytic expressions for the optical matrix elements in certain
limiting cases. When $s\gg L_{e},L_{h}$ and confinement is stong the electron
and hole are in the ground states of their confining potentials and their
wavefunctions are

$\psi_{0}^{(\alpha)}(\bm{r})=\frac{1}{\sqrt{\pi}}\frac{1}{L_{\alpha}}e^{-r^{2}/2L_{\alpha}^{2}},$
(23)

where $\alpha=e,h$ and the origin $\bm{r}=0$ is chosen at the bottom of the
confining potential. The two-body wavefunction is approximately the product of
$\psi_{0}^{(e)}$ and $\psi_{0}^{(h)}$. Then we can do the integral (21) and
find

$A_{eh}^{(0)}(s)\approx\left(\frac{2L_{e}L_{h}}{L_{e}^{2}+L_{h}^{2}}\right)^{2}\exp\left[-\frac{s^{2}}{L_{e}^{2}+L_{h}^{2}}\left(1+\frac{L_{e}^{2}L_{h}^{2}}{4l_{B}^{4}}\right)\right].$
(24)

The dotted curve in Fig. 5 shows that the above approximation is in good
agreement with the numerical results at large $s$.

For $s=0$ and the weak-confinement limit, the wavefunction of the two-particle
system can be separated into center-of-mass and relative motions

$\Psi_{eh}^{(0)}(\bm{r}_{e},\bm{r}_{h})\approx\psi_{\rm cm}(\bm{R})\psi_{\rm
rel}(\bm{r}),$ (25)

where $\bm{R}=(m_{e}\bm{r}_{e}+m_{h}\bm{r}_{h})/M$ and
$\bm{r}=\bm{r}_{e}-\bm{r}_{h}$ are the center-of-mass and relative
coordinates. The center-of-mass motion is the ground state of the harmonic
oscillator:

$\psi_{\rm cm}(\bm{R})=\frac{1}{\sqrt{\pi}}\frac{1}{L}e^{-R^{2}/2L^{2}},$ (26)

where the confining length
$L=\sqrt{\hbar/M\omega}=$1.30\text{\,}\mathrm{n}\mathrm{m}$$. The relative
motion is the ground state of the 2D hydrogen atom

$\psi_{\rm
rel}(\bm{r})=\frac{1}{\sqrt{2\pi}}\frac{4}{a_{B}^{*}}e^{-2r/a_{B}^{*}}.$ (27)

With the wavefunction (25) and the current operator (19) we find

$\langle\Psi_{0}|\hat{j}|\Psi_{eh}^{(0)}\rangle=\int
d\bm{R}\,\Psi_{eh}^{(0)}(\bm{R},\bm{R})\approx 4\sqrt{2}L/a_{B}^{*},$ (28)

and therefore

$\begin{split}A_{eh}^{(0)}(0)=|\langle\Psi_{0}|\hat{j}|\Psi_{eh}^{(0)}\rangle|^{2}\approx\frac{32L^{2}}{a_{B}^{*2}}\approx
38.\end{split}$ (29)

This estimate is much larger than the result from ED calculation, as shown in
Fig. 5. The difference can be understood as follows. Confinement adds a
harmonic repuslive contribution to the relative motion which mixes the ground
state of the 2D hydrogen atom (27) with higher energy states. Since the
wavefunctions of higher-energy hydrogenic states have smaller values at
$\bm{r}=0$, the resulting optical absorption strength is reduced [65]. For the
parameters used in the calculations illustrated in Fig. 5, the confining
potential is far from weak, a lot of high energy states come into play, and
the discrepancy between the exact numerical results and the approximation of
(29) becomes large.

##  IV Trions

Next we study a system of one electron and two holes – a moiré-confined
trion.The three-particle state must be antisymmetric with respect to exchange
of the two holes. We work in a representation in which the orbital and valley
parts of the two-hole states separate, and each is symmetric or antisymmetric
with respect to the permutation. The orbital parts of the basis states are
($n_{h1}>n_{h2}$)

$\displaystyle|n_{e};n_{h1}n_{h2}\rangle_{s}=\frac{1}{\sqrt{2}}|n_{e}\rangle_{e}(|n_{h1}n_{h2}\rangle_{h}+|n_{h2}n_{h1}\rangle_{h}),$
(30)
$\displaystyle|n_{e};n_{h1}n_{h2}\rangle_{a}=\frac{1}{\sqrt{2}}|n_{e}\rangle_{e}(|n_{h1}n_{h2}\rangle_{h}-|n_{h2}n_{h1}\rangle_{h}),$
(31)
$\displaystyle|n_{e};n_{h}n_{h}\rangle_{s}=|n_{e}\rangle_{e}|n_{h}n_{h}\rangle_{h},$
(32)

where each $n$ is shorthand for the two indices $(n_{+},n_{-})$ that
distinguish two-dimensional oscillator single-particle states, the subscripts
$s$ and $a$ stand for symmetric and antisymmetric states, and the kets with
subscripts $e$ and $h$ are direct product states of electrons and holes
without (anti)symmetrization. The valley parts of the basis states are

$\displaystyle|\pm;++\rangle_{s}=|\pm\rangle_{e}|++\rangle_{h},$ (33)
$\displaystyle|\pm;--\rangle_{s}=|\pm\rangle_{e}|--\rangle_{h},$ (34)
$\displaystyle|\pm;+-\rangle_{s}=\frac{1}{\sqrt{2}}|\pm\rangle_{e}(|+-\rangle_{h}+|-+\rangle_{h}),$
(35)
$\displaystyle|\pm;+-\rangle_{a}=\frac{1}{\sqrt{2}}|\pm\rangle_{e}(|+-\rangle_{h}-|-+\rangle_{h}),$
(36)

where $+$ and $-$ stand for valley $+K$ and $-K$.

We are interested in the optical excitations of a state with a single hole to
a trion one-electron-two-hole state. The initial hole is in the ground state
$n_{h}=0$ (by which we mean $n_{h+}=n_{h-}=0$), and we assume without loss of
generality that it is in the $-K$ valley:

$|\Psi_{h}^{(0)}\rangle=|0\rangle_{h}|-\rangle_{h}.$ (37)

In order for the optical matrix elements to be non-zero, the newly added
electron and hole must be in the same valley. There are three possibilities:
$|+;+-\rangle_{a}$, $|+;+-\rangle_{s}$ and $|-;--\rangle_{s}$. Correspondingly
the three-particle state has three possible forms:

$\displaystyle|\Psi_{ehh}^{(1,i)}\rangle$
$\displaystyle=\sum_{n_{e}}\sum_{n_{h1}\geq
n_{h2}}C_{n_{e}n_{h1}n_{h2}}^{(s,i)}|n_{e};n_{h1}n_{h2}\rangle_{s}|+;+-\rangle_{a},$
(38) $\displaystyle|\Psi_{ehh}^{(2,i)}\rangle$
$\displaystyle=\sum_{n_{e}}\sum_{n_{h1}>n_{h2}}C_{n_{e}n_{h1}n_{h2}}^{(a,i)}|n_{e};n_{h1}n_{h2}\rangle_{a}|+;+-\rangle_{s},$
(39) $\displaystyle|\Psi_{ehh}^{(3,i)}\rangle$
$\displaystyle=\sum_{n_{e}}\sum_{n_{h1}>n_{h2}}C_{n_{e}n_{h1}n_{h2}}^{(a,i)}|n_{e};n_{h1}n_{h2}\rangle_{a}|-;--\rangle_{s}.$
(40)

The energy of the system depends only on the orbital part of the state, so the
second and third forms of the three-particle state (39) and (40) with
antisymmetric orbital parts are exactly degenerate in energy. Exchange
statistics suggests that the ground state takes the first form (38), since
there is then no exclusion principle that forbids the two holes from occupying
the same low-energy single-particle orbital ($n_{h1}=n_{h2}=0$). This
expectation is confirmed by our numerical results.

Using the current operator (19), we find the optical absorption strength for
intralayer trions

$\displaystyle A_{ehh}^{(1,i)}$ $\displaystyle=\Big{|}\sum_{n_{e}}\int
d\bm{r}\psi_{n_{e}}^{(e)}(\bm{r})\left[\frac{1}{\sqrt{2}}\sum_{n_{h}>0}C_{n_{e}n_{h}0}^{(s,i)}\psi_{n_{h}}^{(h)}(\bm{r})+C_{n_{e}00}^{(s,i)}\psi_{0}^{(h)}(\bm{r})\right]\exp(-i\frac{eBs}{2\hbar}y)\Big{|}^{2}\equiv
A_{ehh}^{(s,i)},$ (41) $\displaystyle A_{ehh}^{(2,i)}$
$\displaystyle=\frac{1}{2}\,\Big{|}\sum_{n_{e}}\sum_{n_{h}>0}C_{n_{e}n_{h}0}^{(a,i)}\int
d\bm{r}\psi_{n_{e}}^{(e)}(\bm{r})\psi_{n_{h}}^{(h)}(\bm{r})\exp(-i\frac{eBs}{2\hbar}y)\Big{|}^{2}\equiv
A_{ehh}^{(a,i)},$ (42) $\displaystyle A_{ehh}^{(3,i)}$
$\displaystyle=\Big{|}\sum_{n_{e}}\sum_{n_{h}>0}C_{n_{e}n_{h}0}^{(a,i)}\int
d\bm{r}\psi_{n_{e}}^{(e)}(\bm{r})\psi_{n_{h}}^{(h)}(\bm{r})\exp(-i\frac{eBs}{2\hbar}y)\Big{|}^{2}=2A_{ehh}^{(a,i)}.$
(43)

When the polarization of the incident light is random (or equivalently, when
the valley polarization of the initial hole is random), there is equal
probability that an electron-hole pair is created in the same or opposite
valley as the initial hole. Therefore the optical absorption spectrum is

$\begin{split}A_{ehh}(\omega)=&\frac{1}{2}\sum_{p=1}^{3}\sum_{i}A_{ehh}^{(p,i)}\,\delta(\hbar\omega-
E_{ehh}^{(p,i)}-E_{g})\\\
=&\frac{1}{2}\sum_{i}A_{ehh}^{(s,i)}\,\delta(\hbar\omega-
E_{eh}^{(s,i)}-E_{g})\\\
&+\frac{3}{2}\sum_{i}A_{ehh}^{(a,i)}\,\delta(\hbar\omega-
E_{ehh}^{(a,i)}-E_{g}).\end{split}$ (44)

Figure 6: Ground-state trion energy $E_{ehh}^{(0)}$ as a function of the
separation between the electron and hole confining potential minima $s$, for
intralayer ($d=0$) and interlayer ($d=1,$2\text{\,}\mathrm{n}\mathrm{m}$$)
trions represented by different colors. The magnetic field
$B=$10\text{\,}\mathrm{T}$$. The dots and solid lines are numerical results
from ED calculations, and the dotted curves show results from the large-$s$
approximation (45).

Fig. 6 shows the ground-state trion energy $E_{ehh}^{(0)}$ as a function of
the electron-hole separation $s$. Other parameters of the calculations are
chosen to be the same as for the calculations of excitons in Sec. III. The
dotted curve shows the large-$s$ approximation in which the motion of the
electron (holes) is not affected by the presence of the holes (electron) in
the other potential well, so the ground state energy is

$E_{ehh}^{(0)}(s)\approx
E_{e}^{(0)}+E_{hh}^{(0)}-\frac{2e^{2}}{\epsilon\sqrt{s^{2}+d^{2}}},$ (45)

where $E_{hh}^{(0)}$ is the ground state energy of two holes in the same
quantum dot. A similar ED calculation with only holes present yields
$E_{hh}^{(0)}=$224\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$ for the parameters
of Fig. 6. As shown in Fig. 6, the approximation (45) is in good agreement
with the numerical results for $s\gtrsim L_{e},L_{h}$.

Figure 7: Ground state trion binding energy $E_{tb}$ as a function of the
separation between the electron and hole confining potential minima $s$, for
intralayer ($d=0$) and interlayer ($d=1,$2\text{\,}\mathrm{n}\mathrm{m}$$)
trions represented by different colors. The magnetic field
$B=$10\text{\,}\mathrm{T}$$. The dots and solid lines are numerical results
from ED calculations, and the dotted curves show results from the large-$s$
approximation.

Fig. 7 shows the trion binding energy, defined as

$E_{tb}=E_{eh}^{(0)}+E_{h}^{(0)}-E_{ehh}^{(0)},$ (46)

as a function of $s$. At small $s$, an extra hole binds with the electron-hole
pair in such a way as to lower the energy. At large $s$, the repulsion between
the two holes dominates and the energy increases when an extra hole is added
to the system. The crossover from a bound trion to an unbound trion takes
place between 2 to $3\text{\,}\mathrm{n}\mathrm{m}$ for an intralayer trion,
and around $1\text{\,}\mathrm{n}\mathrm{m}$ for an interlayer trion at
interlayer distance $d=$1\text{\,}\mathrm{n}\mathrm{m}$$. An interlayer trion
at $d=$2\text{\,}\mathrm{n}\mathrm{m}$$ is always unbound even for $s=0$.

Figure 8: Optical absorption strength $A_{ehh}^{(0)}$ of intralayer trions as
a function of the separation between the electron and hole confining potential
minima $s$. The magnetic field $B=$10\text{\,}\mathrm{T}$$. The blue dots are
numerical results from ED calculations, and the dotted curve is the large-$s$
approximation (49).

Fig. 8 shows the optical absorption strength $A_{ehh}^{(0)}$ as a function of
$s$. Comparison with Fig. 5 shows that the optical absorption of an intralayer
trion is always weaker than that of an intralayer exciton. This can be
qualitatively understood by considering the $m_{e}\to\infty$ limit at $s=0$
for which the electron stays at the center of its quantum dot while the two
holes move around it. The Coulomb repulsion between the two holes pushes them
away from the electron into higher-energy states, so the overlap between the
electron and hole wavefunctions, and therefore the optical absorption
strength, decreases. For large $s\gg L_{e},L_{h}$, the optical absorption
strength is small and can be calculated from Eq. (41) with

$C^{(s,0)}_{n_{e}n_{h1}n_{h2}}=\delta_{n_{e},0}C^{(s,0)}_{n_{h1}n_{h2}},$ (47)

where $C^{(s,0)}_{n_{h1}n_{h2}}$ (with $n_{h1}>n_{h2}$) are the expansion
coefficients of the two-hole ground state in the symmetric two-hole basis. We
notice from the ED calculation that the two-hole ground state comes
predominantly from the single-particle ground state $|00,00\rangle_{s}$ and
the lowest two excited states with zero total angular momentum
$|11,00\rangle_{s}$ and $|10,01\rangle_{s}$:

$|\Psi_{hh}^{(0)}\rangle\approx
C_{0}|00,00\rangle_{s}+C_{1}|11,00\rangle_{s}+C_{2}|10,01\rangle_{s},$ (48)

with $C_{0}\approx 0.80$ and $C_{1}\approx C_{2}\approx 0.42$. Then from Eq.
(41) we find an analytic expression for the optical absorption strength in the
large-$s$ limit:

$A_{ehh}^{(0)}(s)\approx\Bigg{|}C_{0}-\frac{C_{1}}{\sqrt{2}}\left[\frac{L_{e}^{2}-L_{h}^{2}}{L_{e}^{2}+L_{h}^{2}}+\left(1-\frac{L_{e}^{4}}{4l_{B}^{4}}\right)\frac{L_{h}^{2}s^{2}}{L_{e}^{2}+L_{h}^{2}}\right]\Bigg{|}^{2}\left(\frac{2L_{e}L_{h}}{L_{e}^{2}+L_{h}^{2}}\right)^{2}\exp\left[-\frac{s^{2}}{L_{e}^{2}+L_{h}^{2}}\left(1+\frac{L_{e}^{2}L_{h}^{2}}{4l_{B}^{4}}\right)\right].$
(49)

This approximation is in good agreement with the numerical results at
$s\gtrsim$3\text{\,}\mathrm{n}\mathrm{m}$$ as shown by the dotted curve in
Fig. 8.

##  V Summary and Discussion

In this article we have numerically calculated the energies of exciton and
trion states in moiré superlattices, and the matrix elements for interband
absorption processes with exciton and trion final states. We have modeled the
moiré superlattice by harmonic confinement potentials – an approximation that
is valid when the confinement oscillator lengths are small compared to the
moiré lattice constant. Our calculations are similar to earlier ones [63, 64,
65, 66] that used harmonic confinement potentials to model quantum dots in
semiconductor quantum wells, but differ because of the relevance of the TMD
valley degree of freedom, and by the possibility in the moiré superlattice
case that electrons and holes are localized at different lattice sites.

The relative positions of electron and hole confinement locations is very
important for the optical properties of moiré superlattices. As illustrated in
Fig. 1, the two most likely possibilities are that the confinement locations
are at identical positions in the moiré unit cell, or that they are separated
by a large fraction of the moiré lattice constant. It is not yet known whether
one circumstance or the other is always realized, or whether this important
distinction is interface-dependent. Our calculations show that the two cases
should be easily observable experimentally, since the rotational symmetries
that apply approximately in the common-confinement-location case, lead to a
small number of very strong absorption peaks. In the distinct-confinement-
location case, the number of exciton and trion absorption peaks proliferates
as angular-momentum-conservation selection rules are relaxed, but the overall
absorption strength is still weakened due to electron-hole confinement.
Because the Bohr radius length scale is so short in TMD two-dimensional
semiconductors, the influence of magnetic fields on observable properties is
weak at moderate magnetic field strengths. In fact most of our illustrative
calculations were carried out at a magnetic field $B=$10\text{\,}\mathrm{T}$$,
but are nearly identical to the results we obtain at $B=0$. Magnetic fields
are potentially very useful probes in sorting out the properties of electron-
hole systems formed at particular interfaces, but will require extremely
strong magnetic fields on the $100\text{\,}\mathrm{T}$ scale.

A full theory of the optical absorption spectrum of a doped moiré system is
complicated by importance of correlations among the doped charges present in
the initial state and correlations between electrons, holes, and trions
located in different periods of the moiré superlattice [52]. As a first
approximation, we ignore these complications and assume that at filling factor
$\nu\in(0,1)$, the optical spectrum is simply a weighted average of the
exciton and trion spectra:

$A_{\nu}(\omega)=(1-\nu)A_{eh}(\omega)+\nu A_{ehh}(\omega).$ (50)

with $A_{eh}$ and $A_{ehh}$ given by Eqs. (22) and (44). The numerical results
are shown in Fig. 9, for a large $s=$3\text{\,}\mathrm{n}\mathrm{m}$$ (a) and
a small $s=$1\text{\,}\mathrm{n}\mathrm{m}$$ (b). The magnetic field
$B=$10\text{\,}\mathrm{T}$$, and all other parameters are the same as in Secs.
III and IV. The band gap $E_{g}$ is set to zero in this figure since its only
effect is a constant shift along the $\omega$-axis. The $\delta$-function is
replaced with a Lorentzian function with height 1 and width
$0.2\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$. At
$s=$1\text{\,}\mathrm{n}\mathrm{m}$$ the electrons and holes are strongly
bound and the two- or three-body eigenstates of the system have approximate
rotational symmetry, so most of these states have optical absorption that is
supressed by angular momentum conservation [49, 63]. In contrast, at
$s=$3\text{\,}\mathrm{n}\mathrm{m}$$ the electrons and holes are strongly
localized in the harmonic potentials, and Coulomb interaction distorts the
localized wavefunctions weakly. In this case the selection rules associated
with angular momentum conservation are not relevant because rotational
symmetry is strongly broken. The absorption peaks in Fig. 9(a) are therefore
much more closely spaced than in Fig. 9(b). The peaks for states with nonzero
angular momenta appear in the strong-confinement limit and are split in a
magnetic field (e.g., the exciton peaks near
$40\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ in Fig. 9(a)). The results do
demonstrate that the absorption spectrum of moiré superlattices is strongly
sensitive to electrical modulation of the carrier density.

Figure 9: The optical absorption spectra as the ground state hole moiré
filling factor $\nu$ varies from 0 to 1, assuming that excitons. holes, and
trions centered on different sties are not correlated. These results were
calculated at (a) $s=$3\text{\,}\mathrm{n}\mathrm{m}$$; (b)
$s=$1\text{\,}\mathrm{n}\mathrm{m}$$.

A more accurate description of the optical spectra of moiré materials would
take account of correlations between electrons and holes in different periods
of the moiré pattern, but is beyond the scope of the present work. These
correlations will entangle photons emitted by nearby excitons, and promise the
possibility of realizing unusual single and few photon light emitters.
Intervalley exchange interactions are crucial for the topological properties
of point-particle-like moiré excitons [38] and are expected to play an
important role for our case of spatially distorted excitons and trions as
well. Our work has focused on excitons and trions, but our method also applies
to more complex systems like biexcitons. Because the moiré confinement length
is relatively small, the transition energy between excitons and biexcitons
localized at a single moiré site is expected to be large. The method can also
be extended to study a combination of localized and itinerant particles, which
may help interpret recent experimental work [26, 72] that used excitons in a
TMD monolayer as an optical sensor of correlated insulating states in moiré
materials. We anticipate that close interaction between theory and experiment
will be necessary to sort out all the richness of the optical properties of
moiré superlattices.

###### Acknowledgements.

Y.Z. and A.H.M. acknowledge many helpful conversations with Hui Deng and Feng
Wang. This work was supported by Army Research Office (ARO) Grant #
W911NF-17-1-0312 (MURI).

##  Appendix A Derivation of Coulomb matrix elements in the harmonic
oscillator basis

###  A.1 Single-particle basis

We consider an electron or a hole moving in a two-dimensional quantum dot,
described by the Hamiltonian

$H=\frac{1}{2m}(\bm{p}\pm e\bm{A})+\frac{1}{2}m\omega^{2}\bm{r}^{2},$ (51)

where the upper and lower signs are for electrons and holes respectively. We
work in the symmetric gauge

$A_{x}=-By/2,\quad A_{y}=Bx/2.$ (52)

To simplify notations we define the dimensionless quantities

$\tilde{\bm{r}}=\bm{r}/l_{B},\quad\tilde{\omega}=\omega/\omega_{c},\quad\tilde{H}=H/\hbar\omega_{c},$
(53)

where the magnetic length $l_{B}=\sqrt{\hbar/eB}$ and the cyclotron frequency
$\quad\omega_{c}=eB/m$. The dimensionless Hamiltonian can then be written as

$\tilde{H}=\frac{1}{2}(-i\partial_{\tilde{x}}\mp\frac{\tilde{y}}{2})^{2}+\frac{1}{2}(-i\partial_{\tilde{y}}\pm\frac{\tilde{x}}{2})+\frac{1}{2}\tilde{\omega}^{2}(\tilde{x}^{2}+\tilde{y}^{2}).$
(54)

Now we transform to the complex coordinates

$z=x+iy,\quad z^{*}=x-iy,$ (55)

and their dimensionless counterparts $\tilde{z},\tilde{z}^{*}$ defined in an
obvious way. Then after some algebra the Hamiltonian can be written as

$\tilde{H}=-2\partial_{\tilde{z}^{*}}\partial_{\tilde{z}}+(\frac{1}{8}+\frac{1}{2}\tilde{\omega}^{2})\tilde{z}^{*}\tilde{z}\pm\frac{1}{2}(\tilde{z}\partial_{\tilde{z}}-\tilde{z}^{*}\partial_{\tilde{z}^{*}}).$
(56)

Now we introduce the operators

$\begin{split}b_{+}=\frac{1}{\sqrt{2}}(\frac{\tilde{z}^{*}}{2}+2\partial_{\tilde{z}}),\quad&b_{+}^{\dagger}=\frac{1}{\sqrt{2}}(\frac{\tilde{z}}{2}-2\partial_{\tilde{z}^{*}}),\\\
b_{-}=\frac{1}{\sqrt{2}}(\frac{\tilde{z}}{2}+2\partial_{\tilde{z}^{*}}),\quad&b_{-}^{\dagger}=\frac{1}{\sqrt{2}}(\frac{\tilde{z}^{*}}{2}-2\partial_{\tilde{z}}),\end{split}$
(57)

which satisfy the usual commutation relations of raising and lowering ladder
operators. The Hamiltonian can be written in terms of these operators as

$\tilde{H}=b_{\pm}^{\dagger}b_{\pm}+\frac{1}{2}+\tilde{\omega}^{2}(b_{+}^{\dagger}+b_{-})(b_{+}+b_{-}^{\dagger}),$
(58)

where as before the upper and lower signs (in the subscript) are for electrons
and holes respectively. The Hamiltonian is quadratic in the ladder operators,
and can be diagonalized by a Bogoliubov transformation

$\begin{split}a_{+}&=b_{\pm}\cosh\frac{\eta}{2}+b_{\mp}^{\dagger}\sinh\frac{\eta}{2},\\\
a_{-}&=b_{\mp}\cosh\frac{\eta}{2}+b_{\pm}^{\dagger}\sinh\frac{\eta}{2},\end{split}$
(59)

where

$\tanh\eta=\frac{2\tilde{\omega}^{2}}{2\tilde{\omega}^{2}+1}=\frac{2\omega^{2}}{2\omega^{2}+\omega_{c}^{2}}.$
(60)

The Hamiltonian then takes the diagonal form

$\tilde{H}=\tilde{\Omega}_{+}(a_{+}^{\dagger}a_{+}+\frac{1}{2})+\tilde{\Omega}_{-}(a_{-}^{\dagger}a_{-}+\frac{1}{2}),$
(61)

where

$\tilde{\Omega}_{\pm}=\frac{1}{2}(\sqrt{1+4\tilde{\omega}^{2}}\pm 1).$ (62)

Restoring the physical dimensions, we have the diagonalized Hamiltonian

$H=\hbar\Omega_{+}(a_{+}^{\dagger}a_{+}+\frac{1}{2})+\hbar\Omega_{-}(a_{-}^{\dagger}a_{-}+\frac{1}{2}),$
(63)

with

$\Omega_{\pm}=\sqrt{\omega^{2}+(\frac{\omega_{c}}{2})^{2}}\pm\frac{\omega_{c}}{2}\equiv\Omega\pm\frac{\omega_{c}}{2}.$
(64)

The eigenstates of this Hamiltonian are

$|n_{+}n_{-}\rangle=(a_{+}^{\dagger})^{n_{+}}(a_{-}^{\dagger})^{n_{-}}|0\rangle,$
(65)

with energy

$E_{n_{+}n_{-}}=(n_{+}+\frac{1}{2})\hbar\Omega_{+}+(n_{-}+\frac{1}{2})\hbar\Omega_{-}.$
(66)

For later use, the complex coordinate $z$ can be expressed in terms of the
ladder operators as

$z=\sqrt{2}l_{B}(b_{+}^{\dagger}+b_{-})=\sqrt{2}l_{B}e^{-\eta/2}(a_{\pm}^{\dagger}+a_{\mp})=L(a_{\pm}^{\dagger}+a_{\mp}),$
(67)

where $L=\sqrt{2}l_{B}e^{-\eta/2}=\sqrt{\hbar/m\Omega}$ is the total confining
length. The normalized real-space wavefunction of the energy eigenstates is

$\psi_{n_{+}n_{-}}(r,\theta)=\left[\frac{n!}{(n+|l|)!}\frac{1}{\pi}\frac{1}{L}\right]^{1/2}\left(\frac{r}{L}\right)^{|l|}e^{-r^{2}/2L^{2}}L_{n}^{|l|}(r^{2}/L^{2})\cdot
e^{il\theta},$ (68)

where $n=\min(n_{+},n_{-})$, $l=\pm(n_{+}-n_{-})$,
$\Omega=\sqrt{\omega^{2}+(\omega_{c}/2)^{2}}$, and $L_{n}^{|l|}$ is the
associated Laguerre polynomial.

###  A.2 Coulomb matrix elements

Now we calculate the Coulomb matrix elements in the above basis (generalized
to multiple particles, see Eq. (8)). We start with the most complicated one,
$V_{eh}$ (Eq. (6)), the electron-hole interaction. The calculations for
$V_{ee}$ and $V_{hh}$ are similar and simpler. We first write $V_{eh}$ in
Fourier form (consider only one electron and one hole for simplicity):

$V_{eh}=-\frac{2\pi
e^{2}}{\epsilon}\int\frac{d^{2}q}{(2\pi)^{2}}\,\frac{1}{|q|}e^{-|q|d}e^{iq_{x}s}e^{i(q^{*}z_{e}+qz_{e}^{*})/2}e^{-i(q^{*}z_{h}+qz_{h}^{*})/2},$
(69)

where $q=q_{x}+iq_{y}$ is a complex integration variable. By Eq. (67) we can
write the complex coordinates $z$ in terms of ladder operators
$a_{\pm},a_{\pm}^{\dagger}$, then the matrix elements of $V_{eh}$ can be
written as

$\begin{split}\langle
n_{e+}^{\prime}n_{e-}^{\prime},n_{h+}^{\prime}n_{h-}^{\prime}|V_{eh}|n_{e+}n_{e-},n_{h+}n_{h-}\rangle&=-\frac{2\pi
e^{2}}{\epsilon}\int\frac{d^{2}q}{(2\pi)^{2}}\,\frac{1}{|q|}e^{-|q|d}e^{iq_{x}s}\langle
n_{e+}^{\prime}|e^{iL_{e}(q^{*}a_{e+}^{\dagger}+qa_{e+})}|n_{e+}\rangle\\\
\times\langle
n_{e-}^{\prime}|e^{iL_{e}(q^{*}a_{e-}^{\dagger}+qa_{e-})}|n_{e-}\rangle&\langle
n_{h+}^{\prime}|e^{-iL_{h}(qa_{h+}^{\dagger}+q^{*}a_{h+})}|n_{h+}\rangle\langle
n_{h-}^{\prime}|e^{-iL_{h}(qa_{h-}^{\dagger}+q^{*}a_{h-})}|n_{h-}\rangle,\end{split}$
(70)

Then we need to evaluate matrix elements of the form $\langle
n^{\prime}|e^{i(fa^{\dagger}+f^{*}a)}|n\rangle$. To this end we make use of
the operator identity

$e^{i(fa^{\dagger}+f^{*}a)}=e^{-|f|^{2}/2}e^{ifa^{\dagger}}e^{if^{*}a},$ (71)

which then leads to

$\begin{split}\langle
n^{\prime}|e^{i(fa^{\dagger}+f^{*}a)}|n\rangle&=e^{-|f|^{2}/2}\langle
n^{\prime}|e^{ifa^{\dagger}}e^{if^{*}a}|n\rangle\\\
&=e^{-|f|^{2}/2}\sum_{k^{\prime}=0}^{n^{\prime}}\sum_{k=0}^{n}\frac{(if)^{k^{\prime}}(if^{*})^{k}}{k^{\prime}!k!}\langle
n^{\prime}|(a^{\dagger})^{k^{\prime}}a^{k}|n\rangle\\\
&=e^{-|f|^{2}/2}(ie^{i\varphi})^{n^{\prime}-n}(n^{\prime}!n!)^{1/2}\sum_{k=0}^{\min(n^{\prime},n)}\frac{(-1)^{n-k}|f|^{n^{\prime}+n-2k}}{k!(n^{\prime}-k)!(n-k)!},\end{split}$
(72)

where $\varphi$ is the argument of the complex number $f=|f|e^{i\varphi}$.
With this formula we go back to Eq. (70) and get

$\begin{split}\langle
n_{e+}^{\prime}n_{e-}^{\prime},n_{h+}^{\prime}n_{h-}^{\prime}|V_{eh}|n_{e+}n_{e-},n_{h+}n_{h-}\rangle&=-\sqrt{\frac{\beta_{eh}}{2}}\frac{e^{2}}{\epsilon
l_{B}}\left(\prod_{\alpha\lambda}n_{\alpha\lambda}^{\prime}!n_{\alpha\lambda}!\right)^{1/2}(-1)^{n_{e+}^{\prime}+n_{e-}+n_{h+}^{\prime}+n_{h-}}\\\
&\times\sum_{k_{\alpha\lambda}=0}^{\min(n_{\alpha\lambda}^{\prime},n_{\alpha\lambda})}\left(\prod_{\alpha\lambda}\frac{(-1)^{k_{\alpha\lambda}}(\frac{\beta_{eh}}{\beta_{\alpha}})^{\frac{1}{2}(n_{\alpha\lambda}^{\prime}+n_{\alpha\lambda})-k_{\alpha\lambda}}}{k_{\alpha\lambda}!(n_{\alpha\lambda}^{\prime}-k_{\alpha\lambda})!(n_{\alpha\lambda}-k_{\alpha\lambda})!}\right)\cdot
I,\end{split}$ (73)

where the summation indices $\alpha=e,h$ and $\lambda=\pm$ (the summation sign
really means four summations), and we defined
$\beta_{\alpha}=e^{\eta_{\alpha}}=\sqrt{1+(2\omega_{\alpha}/\omega_{\alpha
c})^{2}}$ and $\beta_{eh}=(\beta_{e}^{-1}+\beta_{h}^{-1})^{-1}$, and the
$(n^{\prime},n,k)$-dependent integral

$I=\int_{0}^{\infty}d\tilde{q}^{2}\,e^{-\tilde{q}^{2}}(\tilde{q}^{2})^{p-\frac{1}{2}}e^{-\tilde{q}\tilde{d}}\cdot\int_{0}^{2\pi}\frac{d\phi}{2\pi}\,(ie^{i\phi})^{l-l^{\prime}}e^{i\tilde{q}\tilde{s}\cos\phi},$
(74)

with the dimensionless variables

$\tilde{q}=\frac{l_{B}}{\sqrt{2\beta_{eh}}}|q|,\quad\tilde{d}=\frac{\sqrt{2\beta_{eh}}}{l_{B}}d,\quad\tilde{s}=\frac{\sqrt{2\beta_{eh}}}{l_{B}}s,$
(75)

the index

$p=\sum_{\alpha\lambda}(\frac{n_{\alpha\lambda}^{\prime}+n_{\alpha\lambda}}{2}-k_{\alpha\lambda})$
(76)

and the total angular momentum quantum numbers

$l=n_{e+}-n_{e-}-n_{h+}+n_{h-},\quad
l^{\prime}=n_{e+}^{\prime}-n_{e-}^{\prime}-n_{h+}^{\prime}+n_{h-}^{\prime}.$
(77)

The integral (74) can be simplified in special cases. When $s=0$, the angular
integral is simply $\delta_{l^{\prime}l}$. When $d$ is also zero, the radial
integral yields a gamma function and the result is

$I=\delta_{l^{\prime}l}\cdot\Gamma(p+\frac{1}{2}).$ (78)

When $s=0$ and $d\neq 0$, the radial integral yields the product of a gamma
function and a (Tricomi’s) confluent hypergeometric function

$I=\delta_{l^{\prime}l}\cdot
4^{-p}\,\Gamma(2p+1)\,U(p+\frac{1}{2},\frac{1}{2},\frac{\tilde{d}^{2}}{4}).$
(79)

However, in the general case where $s\neq 0$, only the radial part can be
integrated in closed form and the angular integral remains:

$I=4^{-p}\,\Gamma(2p+1)\int_{0}^{2\pi}\frac{d\phi}{2\pi}\,(ie^{i\phi})^{l-l^{\prime}}\,U(p+\frac{1}{2},\frac{1}{2},\frac{1}{4}(\tilde{d}-i\tilde{s}\cos\phi)^{2}).$
(80)

In practical calculations the integral has to be done numerically. The matrix
elements of $V_{ee}$ and $V_{hh}$ can be calculated in the same way, although
the result looks simpler since $d=s=0$ for the same type of particles. The
result is (without antisymmetrization)

$\begin{split}\langle
n_{1+}^{\prime}n_{1-}^{\prime},n_{2+}^{\prime}n_{2-}^{\prime}|V_{\alpha\alpha}|n_{1+}n_{1-},n_{2+}n_{2-}\rangle&=\frac{\sqrt{\beta_{\alpha}}}{2}\frac{e^{2}}{\epsilon
l_{B}}\left(\prod_{i\lambda}n_{i\lambda}^{\prime}!n_{i\lambda}!\right)^{1/2}(-1)^{n_{1+}^{\prime}+n_{1-}+n_{2+}^{\prime}+n_{2-}}\delta_{l^{\prime}l}\\\
&\times\sum_{k_{i\lambda}=0}^{\min(n_{i\lambda}^{\prime},n_{i\lambda})}\left(\prod_{i\lambda}\frac{(-1)^{k_{i\lambda}}2^{-p}}{k_{i\lambda}!(n_{i\lambda}^{\prime}-k_{i\lambda})!(n_{i\lambda}-k_{i\lambda})!}\right)\cdot\Gamma(p+\frac{1}{2}),\end{split}$
(81)

where $\alpha=e$ or $h$, the summation indices $i=1,2$ is the particle label
(the $e/h$ label is omitted), $\lambda=\pm$, and we defined

$p=\sum_{i\lambda}(\frac{n_{i\lambda}^{\prime}+n_{i\lambda}}{2}-k_{i\lambda}),$
(82) $l=n_{1+}-n_{1-}+n_{2+}-n_{2-},\quad
l^{\prime}=n_{1+}^{\prime}-n_{1-}^{\prime}+n_{2+}^{\prime}-n_{2-}^{\prime}.$
(83)

##  Appendix B The current operator

In the two-band model the current operator is (naively) proportional to
$\sum_{\bm{k}}(e_{\bm{k}}^{\dagger}h_{-\bm{k}}^{\dagger}+h_{-\bm{k}}e_{\bm{k}})$,
which in real space takes the form

$\hat{j}\sim\int
d\bm{r}[e^{\dagger}(\bm{r})h^{\dagger}(\bm{r})+h(\bm{r})e(\bm{r})],$ (84)

where $e$ and $h$ are annihilation operators of electron and hole states. Note
however that in writing Eqs. (2) and (3) with
$\bm{A}_{\alpha}(\bm{r}_{\alpha})=\frac{1}{2}\bm{B}\times\bm{r}_{\alpha}$ for
both $\alpha=e$ and $h$, and with $\bm{r}_{e}$ and $\bm{r}_{h}$ measured from
different origins (see Fig. 2), we have implicitly performed a gauge
transformation

$\bm{A}_{e}(\bm{r}_{e})\to\bm{A}_{e}(\bm{r}_{e})-\frac{Bs}{2}\hat{y},\quad\psi(\\{\bm{r}_{ei}\\},\\{\bm{r}_{hj}\\})\to\psi(\\{\bm{r}_{ei}\\},\\{\bm{r}_{hj}\\})\exp\left(i\frac{eBs}{2\hbar}\sum_{i}y_{ei}\right),$
(85)

compared to the gauge where $\bm{A}(\bm{r})=\frac{1}{2}\bm{B}\times\bm{r}$ for
both electrons and holes, with the origin set at the confining potential
minimum for holes. After the gauge transformation the field operator
$e(\bm{r})$ gains an $\bm{r}$-dependent phase factor, so in this gauge the
current operator becomes

$\hat{j}\sim\int
d\bm{r}\left[e^{\dagger}(\bm{r})h^{\dagger}(\bm{r})\exp\left(i\frac{eBs}{2\hbar}y\right)+h(\bm{r})e(\bm{r})\exp\left(-i\frac{eBs}{2\hbar}y\right)\right].$
(86)

##  Appendix C Basis cutoff

$s$/nm | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---|---|---|---
$N_{x0}$ | 150000 | 3000 | 3000 | 2000 | 2000 | 500 | 200 | 100 | 100
$N_{t0}$ | 500000 | 10000 | 8000 | 8000 | 8000 | 4000 | 4000 | 4000 | 4000
$N_{x1}$ | 20000 | 2000 | 1000 | 500 | 200 | 100 | 100 | 100 | 100
$N_{t1}$ | 50000 | 10000 | 5000 | 5000 | 4000 | 4000 | 4000 | 4000 | 4000
$N_{x2}$ | 10000 | 2000 | 1000 | 500 | 200 | 100 | 100 | 100 | 100
$N_{t2}$ | 20000 | 4000 | 4000 | 4000 | 4000 | 4000 | 4000 | 4000 | 4000
Table 1: Number of basis states for the ED calculations for excitons
($N_{x0},N_{x1},N_{x2}$ for $d=0,1,$2\text{\,}\mathrm{n}\mathrm{m}$$) and
trions ($N_{t0},N_{t1},N_{t2}$) at different in-plane electron-hole
separations $s$.

The energies and optical absorption strengths are obtained by numerically
diagonalizing the Hamiltonian (1) as a matrix in the harmonic oscillator
basis, which has to be cut off under some energy scale. The number of basis
states used in the ED calculations are shown in Table 1. With the exception of
the $s=0$ case, the number of basis states is chosen such that the ground-
state energy is accurate within $1\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ (in
fact the error is much smaller for large $s$), and the optical absorption
strength (for $d=0$) is accurate within 1%. In the magnetic harmonic
oscillator basis, convergence is faster for larger $s$. For $s=0$, the ED
calculations suffer from poor convergence with the number of basis states, as
already noted in Ref. [64] where the exciton problem is solved instead in the
center-of-mass and relative coordinates. However, since we are mostly
interested in the $s\neq 0$ case, and for its easier generalization to more
complicated systems like trions and biexcitons, we solve the $s=0$ case also
in the harmonic oscillator basis, but with a huge number of basis states. The
calculation is manageable when we take advantage of the rotational symmetry of
the system and only diagonalize the block with total angular momentum zero
(which contains the ground state and is optically active).

## References

  * Bistritzer and MacDonald [2011] R. Bistritzerand A. H. MacDonald, Moiré bands in twisted double-layer graphene, Proceedings of the National Academy of Sciences 108, 12233 (2011).
  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, _et al._ , Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature 556, 80 (2018a).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43 (2018b).
  * Lu  _et al._ [2019] X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, _et al._ , Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene, Nature 574, 653 (2019).
  * Yankowitz  _et al._ [2019] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Tuning superconductivity in twisted bilayer graphene, Science 363, 1059 (2019).
  * Po  _et al._ [2018] H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Origin of mott insulating behavior and superconductivity in twisted bilayer graphene, Physical Review X 8, 031089 (2018).
  * Koshino  _et al._ [2018] M. Koshino, N. F. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, Maximally localized Wannier orbitals and the extended Hubbard model for twisted bilayer graphene, Physical Review X 8, 031087 (2018).
  * Kerelsky  _et al._ [2019] A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian, M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, _et al._ , Maximized electron interactions at the magic angle in twisted bilayer graphene, Nature 572, 95 (2019).
  * Serlin  _et al._ [2020] M. Serlin, C. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. Young, Intrinsic quantized anomalous hall effect in a moiré heterostructure, Science 367, 900 (2020).
  * Kang and Vafek [2018] J. Kangand O. Vafek, Symmetry, maximally localized Wannier states, and a low-energy model for twisted bilayer graphene narrow bands, Physical Review X 8, 031088 (2018).
  * Isobe  _et al._ [2018] H. Isobe, N. F. Yuan, and L. Fu, Unconventional superconductivity and density waves in twisted bilayer graphene, Physical Review X 8, 041041 (2018).
  * Wu  _et al._ [2018a] F. Wu, A. MacDonald, and I. Martin, Theory of phonon-mediated superconductivity in twisted bilayer graphene, Physical review letters 121, 257001 (2018a).
  * Lian  _et al._ [2019] B. Lian, Z. Wang, and B. A. Bernevig, Twisted bilayer graphene: a phonon-driven superconductor, Physical review letters 122, 257002 (2019).
  * Xie and MacDonald [2020] M. Xieand A. H. MacDonald, Nature of the correlated insulator states in twisted bilayer graphene, Physical review letters 124, 097601 (2020).
  * Cao  _et al._ [2020] Y. Cao, D. Rodan-Legrain, O. Rubies-Bigorda, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Tunable correlated states and spin-polarized phases in twisted bilayer–bilayer graphene, Nature 583, 215 (2020).
  * Andrei and MacDonald [2020] E. Y. Andreiand A. H. MacDonald, Graphene bilayers with a twist, Nature Materials 19, 1265 (2020).
  * Balents  _et al._ [2020] L. Balents, C. R. Dean, D. K. Efetov, and A. F. Young, Superconductivity and strong correlations in moiré flat bands, Nature Physics 16, 725 (2020).
  * Zondiner  _et al._ [2020] U. Zondiner, A. Rozen, D. Rodan-Legrain, Y. Cao, R. Queiroz, T. Taniguchi, K. Watanabe, Y. Oreg, F. von Oppen, A. Stern, _et al._ , Cascade of phase transitions and dirac revivals in magic-angle graphene, Nature 582, 203 (2020).
  * Wong  _et al._ [2020] D. Wong, K. P. Nuckolls, M. Oh, B. Lian, Y. Xie, S. Jeon, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Cascade of electronic transitions in magic-angle twisted bilayer graphene, Nature 582, 198 (2020).
  * Bultinck  _et al._ [2020] N. Bultinck, E. Khalaf, S. Liu, S. Chatterjee, A. Vishwanath, and M. P. Zaletel, Ground state and hidden symmetry of magic-angle graphene at even integer filling, Physical Review X 10, 031034 (2020).
  * Park  _et al._ [2021] J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Tunable strongly coupled superconductivity in magic-angle twisted trilayer graphene, Nature 590, 249 (2021).
  * Khalaf  _et al._ [2021] E. Khalaf, S. Chatterjee, N. Bultinck, M. P. Zaletel, and A. Vishwanath, Charged skyrmions and topological origin of superconductivity in magic-angle graphene, Science advances 7, eabf5299 (2021).
  * Wu  _et al._ [2018b] F. Wu, T. Lovorn, E. Tutuc, and A. H. MacDonald, Hubbard model physics in transition metal dichalcogenide moiré bands, Physical review letters 121, 026402 (2018b).
  * Regan  _et al._ [2020] E. C. Regan, D. Wang, C. Jin, M. I. B. Utama, B. Gao, X. Wei, S. Zhao, W. Zhao, Z. Zhang, K. Yumigeta, _et al._ , Mott and generalized Wigner crystal states in WSe2/WS2 moiré superlattices, Nature 579, 359 (2020).
  * Tang  _et al._ [2020] Y. Tang, L. Li, T. Li, Y. Xu, S. Liu, K. Barmak, K. Watanabe, T. Taniguchi, A. H. MacDonald, J. Shan, _et al._ , Simulation of hubbard model physics in WSe2/WS2 moiré superlattices, Nature 579, 353 (2020).
  * Xu  _et al._ [2020] Y. Xu, S. Liu, D. A. Rhodes, K. Watanabe, T. Taniguchi, J. Hone, V. Elser, K. F. Mak, and J. Shan, Correlated insulating states at fractional fillings of moiré superlattices, Nature 587, 214 (2020).
  * Wang  _et al._ [2020] L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes, C. Tan, M. Claassen, D. M. Kennes, Y. Bai, B. Kim, _et al._ , Correlated electronic phases in twisted bilayer transition metal dichalcogenides, Nature materials 19, 861 (2020).
  * Kennes  _et al._ [2021] D. M. Kennes, M. Claassen, L. Xian, A. Georges, A. J. Millis, J. Hone, C. R. Dean, D. Basov, A. N. Pasupathy, and A. Rubio, Moiré heterostructures as a condensed-matter quantum simulator, Nature Physics 17, 155 (2021).
  * Jin  _et al._ [2021] C. Jin, Z. Tao, T. Li, Y. Xu, Y. Tang, J. Zhu, S. Liu, K. Watanabe, T. Taniguchi, J. C. Hone, _et al._ , Stripe phases in $\rm{WSe_{2}/WS_{2}}$ moiré superlattices, Nature Materials , 1 (2021).
  * Huang  _et al._ [2021] X. Huang, T. Wang, S. Miao, C. Wang, Z. Li, Z. Lian, T. Taniguchi, K. Watanabe, S. Okamoto, D. Xiao, _et al._ , Correlated insulating states at fractional fillings of the $\rm{WS_{2}/WSe_{2}}$ moiré lattice, Nature Physics 17, 715 (2021).
  * Chu  _et al._ [2020] Z. Chu, E. C. Regan, X. Ma, D. Wang, Z. Xu, M. I. B. Utama, K. Yumigeta, M. Blei, K. Watanabe, T. Taniguchi, _et al._ , Nanoscale conductivity imaging of correlated electronic states in $\rm{Wse_{2}/WS_{2}}$ moiré superlattices, Physical review letters 125, 186803 (2020).
  * Pan  _et al._ [2020] H. Pan, F. Wu, and S. D. Sarma, Quantum phase diagram of a Moiré-Hubbard model, Physical Review B 102, 201104 (2020).
  * Morales-Durán  _et al._ [2021] N. Morales-Durán, A. H. MacDonald, and P. Potasz, Metal-insulator transition in transition metal dichalcogenide heterobilayer moiré superlattices, Physical Review B 103, L241110 (2021).
  * Zhang  _et al._ [2021] Y. Zhang, T. Liu, and L. Fu, Electronic structures, charge transfer, and charge order in twisted transition metal dichalcogenide bilayers, Physical Review B 103, 155142 (2021).
  * Bi and Fu [2021] Z. Biand L. Fu, Excitonic density wave and spin-valley superfluid in bilayer transition metal dichalcogenide, Nature communications 12, 1 (2021).
  * Li  _et al._ [2021a] H. Li, S. Li, E. C. Regan, D. Wang, W. Zhao, S. Kahn, K. Yumigeta, M. Blei, T. Taniguchi, K. Watanabe, _et al._ , Imaging generalized wigner crystal states in a WSe2/WS2 moiré superlattice, arXiv preprint arXiv:2106.10599 (2021a).
  * Angeli and MacDonald [2021] M. Angeliand A. H. MacDonald, $\Gamma$ valley transition metal dichalcogenide moiré bands, Proceedings of the National Academy of Sciences 118 (2021).
  * Wu  _et al._ [2017] F. Wu, T. Lovorn, and A. H. MacDonald, Topological exciton bands in moiré heterojunctions, Physical review letters 118, 147401 (2017).
  * Wu  _et al._ [2018c] F. Wu, T. Lovorn, and A. MacDonald, Theory of optical absorption by interlayer excitons in transition metal dichalcogenide heterobilayers, Physical Review B 97, 035306 (2018c).
  * Yu  _et al._ [2017] 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).
  * Tran  _et al._ [2019] K. Tran, G. Moody, F. Wu, X. Lu, J. Choi, K. Kim, A. Rai, D. A. Sanchez, J. Quan, A. Singh, _et al._ , Evidence for moiré excitons in van der Waals heterostructures, Nature 567, 71 (2019).
  * Jin  _et al._ [2019] C. Jin, E. C. Regan, A. Yan, M. I. B. Utama, D. Wang, S. Zhao, Y. Qin, S. Yang, Z. Zheng, S. Shi, _et al._ , Observation of moiré excitons in WSe2/WS2 heterostructure superlattices, Nature 567, 76 (2019).
  * Seyler  _et al._ [2019] K. L. Seyler, P. Rivera, H. Yu, N. P. Wilson, E. L. Ray, D. G. Mandrus, J. Yan, W. Yao, and X. Xu, Signatures of moiré-trapped valley excitons in MoSe2/WSe2 heterobilayers, Nature 567, 66 (2019).
  * Alexeev  _et al._ [2019] 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, _et al._ , Resonantly hybridized excitons in moiré superlattices in van der Waals heterostructures, Nature 567, 81 (2019).
  * Zhang  _et al._ [2018] N. Zhang, A. Surrente, M. Baranowski, D. K. Maude, P. Gant, A. Castellanos-Gomez, and P. Plochocka, Moiré intralayer excitons in a $\rm{MoSe_{2}/MoS_{2}}$ heterostructure, Nano letters 18, 7651 (2018).
  * Shimazaki  _et al._ [2020] 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).
  * Zhang  _et al._ [2020] L. Zhang, Z. Zhang, F. Wu, D. Wang, R. Gogna, S. Hou, K. Watanabe, T. Taniguchi, K. Kulkarni, T. Kuo, _et al._ , Twist-angle dependence of moiré excitons in WS2/MoSe2 heterobilayers, Nature communications 11, 1 (2020).
  * Baek  _et al._ [2020] H. Baek, M. Brotons-Gisbert, Z. Koong, A. Campbell, M. Rambach, K. Watanabe, T. Taniguchi, and B. D. Gerardot, Highly energy-tunable quantum light from moiré-trapped excitons, Science advances 6, eaba8526 (2020).
  * Ruiz-Tijerina  _et al._ [2020] D. A. Ruiz-Tijerina, I. Soltero, and F. Mireles, Theory of moiré localized excitons in transition metal dichalcogenide heterobilayers, Physical Review B 102, 195403 (2020).
  * Liu  _et al._ [2021] E. Liu, E. Barré, J. van Baren, M. Wilson, T. Taniguchi, K. Watanabe, Y.-T. Cui, N. M. Gabor, T. F. Heinz, Y.-C. Chang, _et al._ , Signatures of moiré trions in WSe2/MoSe2 heterobilayers, Nature 594, 46 (2021).
  * Brotons-Gisbert  _et al._ [2021] M. Brotons-Gisbert, H. Baek, A. Campbell, K. Watanabe, T. Taniguchi, and B. D. Gerardot, Moiré-trapped interlayer trions in a charge-tunable WSe2/MoSe2 heterobilayer, arXiv preprint arXiv:2101.07747 (2021).
  * Baek  _et al._ [2021] H. Baek, M. Brotons-Gisbert, A. Campbell, K. Watanabe, T. Taniguchi, and B. D. Gerardot, Optical read-out of Coulomb staircases in a moiré superlattice via trapped interlayer trions, arXiv preprint arXiv:2102.01358 (2021).
  * Brem  _et al._ [2020] S. Brem, C. Linderälv, P. Erhart, and E. Malic, Tunable phases of moiré excitons in van der Waals heterostructures, Nano letters 20, 8534 (2020).
  * Shabani  _et al._ [2021] 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, Nature Physics 17, 720 (2021).
  * Li  _et al._ [2021b] H. Li, S. Li, M. H. Naik, J. Xie, X. Li, J. Wang, E. Regan, D. Wang, W. Zhao, S. Zhao, _et al._ , Imaging moiré flat bands in three-dimensional reconstructed wse 2/ws 2 superlattices, Nature materials , 1 (2021b).
  * Li  _et al._ [2021c] E. Li, J.-X. Hu, X. Feng, Z. Zhou, L. An, K. T. Law, N. Wang, and N. Lin, Lattice reconstruction induced multiple ultra-flat bands in twisted bilayer WSe2, arXiv preprint arXiv:2103.06479 (2021c).
  * Naik and Jain [2018] M. H. Naikand M. Jain, Ultraflatbands and shear solitons in moiré patterns of twisted bilayer transition metal dichalcogenides, Physical review letters 121, 266401 (2018).
  * Guo  _et al._ [2020] H. Guo, X. Zhang, and G. Lu, Shedding light on moiré excitons: A first-principles perspective, Science advances 6, eabc5638 (2020).
  * Naik  _et al._ [2020] M. H. Naik, S. Kundu, I. Maity, and M. Jain, Origin and evolution of ultraflat bands in twisted bilayer transition metal dichalcogenides: Realization of triangular quantum dots, Physical Review B 102, 075413 (2020).
  * Zhan  _et al._ [2020] Z. Zhan, Y. Zhang, P. Lv, H. Zhong, G. Yu, F. Guinea, J. Á. Silva-Guillén, and S. Yuan, Tunability of multiple ultraflat bands and effect of spin-orbit coupling in twisted bilayer transition metal dichalcogenides, Physical Review B 102, 241106 (2020).
  * Maity  _et al._ [2021] I. Maity, P. K. Maiti, H. Krishnamurthy, and M. Jain, Reconstruction of moiré lattices in twisted transition metal dichalcogenide bilayers, Physical Review B 103, L121102 (2021).
  * Kundu  _et al._ [2021] S. Kundu, I. Maity, R. Bajaj, H. Krishnamurthy, and M. Jain, Atomic reconstruction and flat bands in strain engineered transition metal dichalcogenide bilayer moiré systems, arXiv preprint arXiv:2108.01112 (2021).
  * Wojs and Hawrylak [1995] A. Wojsand P. Hawrylak, Negatively charged magnetoexcitons in quantum dots, Physical Review B 51, 10880 (1995).
  * Halonen  _et al._ [1992] V. Halonen, T. Chakraborty, and P. Pietiläinen, Excitons in a parabolic quantum dot in magnetic fields, Physical Review B 45, 5980 (1992).
  * Que [1992] W. Que, Excitons in quantum dots with parabolic confinement, Physical Review B 45, 11036 (1992).
  * Maksym and Chakraborty [1990] P. Maksymand T. Chakraborty, Quantum dots in a magnetic field: Role of electron-electron interactions, Physical review letters 65, 108 (1990).
  * Yu  _et al._ [2014] H. Yu, G.-B. Liu, P. Gong, X. Xu, and W. Yao, Dirac cones and Dirac saddle points of bright excitons in monolayer transition metal dichalcogenides, Nature communications 5, 1 (2014).
  * Glazov  _et al._ [2014] M. Glazov, T. Amand, X. Marie, D. Lagarde, L. Bouet, and B. Urbaszek, Exciton fine structure and spin decoherence in monolayers of transition metal dichalcogenides, Physical Review B 89, 201302 (2014).
  * Yu and Wu [2014] T. Yuand M. Wu, Valley depolarization due to intervalley and intravalley electron-hole exchange interactions in monolayer MoS2, Physical Review B 89, 205303 (2014).
  * Wu  _et al._ [2015] F. Wu, F. Qu, and A. H. Macdonald, Exciton band structure of monolayer MoS2, Physical Review B 91, 075310 (2015).
  * Kormányos  _et al._ [2015] A. Kormányos, G. Burkard, M. Gmitra, J. Fabian, V. Zólyomi, N. D. Drummond, and V. Fal’ko, $k\cdot p$ theory for two-dimensional transition metal dichalcogenide semiconductors, 2D Materials 2, 022001 (2015).
  * Gu  _et al._ [2021] J. Gu, L. Ma, S. Liu, K. Watanabe, T. Taniguchi, J. C. Hone, J. Shan, and K. F. Mak, Dipolar excitonic insulator in a moiré lattice, arXiv preprint arXiv:2108.06588 (2021).