# Effect of intravalley and intervalley electron-hole exchange on the
nonlinear optical response of monolayer MoSe2

N.H. Kwong  Wyant College of Optical Sciences, The University of Arizona,
Tucson, AZ 85721  J.R. Schaibley  Department of Physics, The University of
Arizona, Tucson, AZ 85721  R. Binder  Wyant College of Optical Sciences, The
University of Arizona, Tucson, AZ 85721  Department of Physics, The University
of Arizona, Tucson, AZ 85721

###### Abstract

The coherent third-order nonlinear response of monolayer transition-metal
dichalcogenide (TMD) semiconductors, such as MoSe2 is dominated by the
nonlinear exciton response, as well as biexciton and trion resonances. The
fact that these resonances may be spectrally close together makes
identification of the signatures, for example in differential transmission
(DT), challenging. Instead of focusing on explaining a given set of
experimental data, a systematic study aimed at elucidating the roles of
intravalley and intervalley electron-hole (e-h) exchange on the DT spectra is
presented. Previous works have shown that the e-h exchange introduces a linear
leading-order term in the exciton dispersion. Based on a generalized Lippmann-
Schwinger equation, we show that the presence of this linear dispersion term
can reduce the biexciton binding energy to zero, contrary to the conventional
situation of quadratic dispersion where an arbitrarily weak (well-behaved)
attractive interaction always supports bound state(s). The effects of spin-
scattering and the spin-orbit interaction caused by e-h exchange is also
clarified, and the DT lineshape at the exciton and trion resonance are studied
as a function of e-h exchange strength. In particular, as the exciton
lineshape is determined by the interplay of linear exciton susceptibility and
the bound-state two-exciton resonance in the T-matrix, the lineshape at the
trion is similarly determined by the interplay of the linear trion
susceptibility and the bound-state exciton-trion resonance in the T-matrix.

##  I Introduction

The physics of monolayer transition-metal dichalcogenide (TMD) semiconductors,
such as MoSe2, has has been widely investigated. Review articles include Refs.
wang-etal.12, ; schaibley-etal.16, ; choi-etal.16, ; manzeli-etal.17, ;
berkelbach-reichman.18, ; mueller-malic.2018, . Examples of studies of these
materials include their growth, preparation and structural properties lieth.77
; duerloo-etal.14 ; li-etal.16 ; hu-etal.18 , as well as thermoelectric
wickramaratne-etal.14 , piezoelectric duerloo-etal.12 , electronic
mattheiss.73 ; mak-etal.10 ; ramasubtamaniam-etal.11 ; xiao-etal.12 ;
cheiwchanchamnangij-lambrecht.12 ; liu-etal.13 ; li-niu.13 ; kosmider-etal.13
; mai-etal.13 ; macneill-etal.15 ; dery.16 , electron transport and transfer
li-etal.17 ; ceballos-etal.17nanolett ; ceballos-etal.17prm ,
mechanical/phononic ataca-etal.11 ; molina-wirtz.11 ; sanchez-wirtz.11 ;
kaasbjerg-etal.12 ; kaasbjerg-etal.13 , and optical, exciton and trion
properties eda-etal.11 ; jones-etal.13 ; qiu-etal.13 ; berkelbach-etal.13 ;
liu-etal.14apl ; liu-etal.14 ; chernikov-etal.14 ; ye-etal.14 ; yu-etal.14 ;
steinhoff-etal.14 ; tran-etal.14 ; bellus-etal.15 ; wu-etal.15prb ; you-
etal.15 ; stroucken-koch.15 ; moody-etal.15 ; kylanpaa-komsa.15 ; jones-
etal.15 ; hsu-etal.15 ; yu-etal.15 ; rivera-etal.16 ; hao-etal.16NaturePhysics
; singh-etal.16 ; hao-etal.16nanolett ; hao-etal.17 ; chow-etal.17 ; scuri-
etal.18 ; meckbach-etal.2018prb ; stier-etal.18 ; villari-etal.2018 ; vantuan-
etal.2018 ; mahon-etal.2019 , and also exploration for future device
applications wang-etal.12 ; jariwala-etal.14 ; li-etal.15 ; seyler-etal.15 ;
ye-etal.15 ; wu-etal.15Nature ; li-etal.17 . Similar to conventional III-V
semiconductors, TMDs have direct electronic bandgaps and host excitons (albeit
with larger exciton binding energies than their III-V counterparts). A crucial
aspect that is different between III-V and TMD semiconductors is the size of
the electron-hole (e-h) exchange, which is on the order of meV in TMDs while
on $\mu$eV in III-V semiconductors. Another aspect that is, at least in
practice, different between typical III-V and TMD semiconductors is the
presence of trion resonances below the exciton in TMDs. This is a consequence
of intentional or unintentional doping of TMDs. The spectroscopic challenge
that arises from the trion resonances is that they are generally close to the
expected position of the biexciton resonance, a few tens of meV below the
exciton. This makes identification of the resonances, which would be easy if
they would be substantially spectrally separated, more difficult. It also
makes the identification of lasing processes, that have been observed in TMDs
ye-etal.15 ; wu-etal.15Nature ; li-etal.17 , more difficult, as there is no
strongly spectrally isolated signature that identifies the lasing as
excitonic, trion-assisted or biexciton-assisted lasing. Previous work hao-
etal.16nanolett ; hao-etal.17 ; steinhoff-etal.2018  has been successful at
using microscopic theories, including e-h exchange steinhoff-etal.2018 , to
identify spectral signatures in experimental differential transmission or
absorption spectra to originate from trions or biexcitons. Also, the effect of
e-h exchange interaction on the exciton dispersion is by now well understood
yu-etal.14 ; wu-etal.15prb ; qiu-etal.15 ; steinhoff-etal.2018 ; deilmann-
thygesen.19 ; schneider-etal.2019 ; hong-etal.20 , as are intervalley exciton
scattering dynamics jiang-etal.2021 . But a systematic study explaining the
underlying principles of how e-h exchange influences the spectral positions
and lineshape in differential transmission (DT) is still lacking.

In the following, we use a T-matrix model and develop a generalized Lippmann-
Schwinger equation engelbrecht-randeria.92 ; takayama-etal.02 , and combine it
with the theory of differential transmission in coherent third-order nonlinear
response regime takayama-etal.04  (so-called $\chi^{(3)}$), to elucidate the
effects of both intravalley and interalley e-h exchange on the nonlinear
exciton response and the resonances corresponding to the various two-particle
complexes, where by ‘particle’ we mean exciton and trion, i.e. on the
2-particle complexes of biexcitons, exciton-trion and trion-trion bound
states. To keep the theory simple and transparent, we use a separable
potential for the 2-particle interaction potential.

There are two central aspects to this analysis. First, the e-h exchange can
modify the dispersion relation of the particles (i.e. the center-of-mass
dispersion of the excitons and trions). It has been shown theoretically wu-
etal.15prb ; qiu-etal.15 ; steinhoff-etal.2018 ; deilmann-thygesen.19  that
the e-h exchange interaction leads to a term linear in the exciton momentum as
the leading order in the dispersion relation of certain branch(es) of the
lowest-lying exciton states. The presence of this linear dispersion behavior
has received some experimental support schneider-etal.2019 ; hong-etal.20 . We
study how this change in the analytic property of the dispersion at low
momenta affects the ability for the particles to form bound states. We first
consider a simple model of two particles with linear dispersion interacting
via an attractive separable potential in a single-channel setting in 2D. We
show analytically that there exists a threshold that the strength of the
interaction needs to exceed for it to support a bound state. This is in
contrast to the conventional case of particles with quadratic energy-momentum
dispersion in 2D where an attractive interaction satisfying some broad
conditions always supports bound state(s) no matter how weak the interaction
is (e.g. see Refs. simon.76, ; chadan-etal.03, ; yang-de-llano.89, ;
apenko.98, ). Typically, replacing the quadratic dispersion with a linear one
weakens the binding of the bound state. We then show that the effect of
weakening or precluding a bound biexciton also holds in a model that has the
spin-valley configuration structure of excitons (and trions) in a TMD
(specifically ${\rm MoSe}_{2}$) where one exciton branch has a linear leading-
order dispersion.

Second, the intervalley e-h exchange results in the coupling of the two
different spin states of each particle. In other words, the spin does not
remain a good quantum number. This coupling also involves spin-orbit coupling,
since, say, a scattering of two particles from the initial state of both
particle being spin-down states to the final state of both being spin-up needs
to be compensated by a corresponding decrease of the orbital angular momentum
(OAM) of the the two-particle complex. Of course, in the lineshape of the DT
all these effects are entangled, but we provide an analysis where we
systematically switch on the various physical processes, thus clarifying their
respective role on the DT spectra.

The paper is organized as follows. In Sec. II we first analyze a scalar model
of the T-matrix and obtain an an analytical solution for particles with linear
dispersion. We also formulate the general T-matrix theory applicable to our
system and derive and expression for the DT in terms of the T-matrix. In Sec.
III we present numerical results for the DT and the relevant T-matric
components. We summarize our findings in Sec. IV.

##  II Theoretical formalism

In the first subsection, we show by a simple single-channel model that two
particles in two dimensions with linear kinetic energy dispersion and
interacting through an attractive force do not necessarily have a bound state.
A threshold may exist which the strength of the attraction needs to exceed in
order for bound state(s) to be supported. This is in contrast to the common
case of quadratic dispersion where an attractive interaction, no matter how
weak, always supports bound state(s) in 2D. In formulating the problem, we use
the Lippmann-Schwinger equation for the two-body scattering T-matrix and
assume a separable interaction. In the second subsection, we generalize the
T-matrix formalism to a model with multiple particle species and spin states,
which is applied to the system of excitons and trions in monolayer Transition
metal dichalcogenides (TMD). Our model configuration is appropriate for ${\rm
MoSe}_{2}$. Electron-hole (e-h) exchanges, both intravalley and intervalley,
are included in the single-particle part of the Hamiltonian. In the last
subsection, the third order ($\chi^{(3)}$) differential transmission of light
through a TMD layer in a pump-probe setting is related to the T-matrix
developed in the second subsection. Numerical results showing how the bound
state energies, the differential transmission etc. are affected by using
dispersions with linear leading-order terms in the e-h exchanges are
investigated in Sec.III.

###  II.1 Single-channel bound states in 2D with linear free-particle
dispersion

We consider the scattering of two particles in the frame where the total
momentum is zero. The scattering T-matrix
$T(\mathbf{k},\mathbf{k}^{\prime},\Omega)$ is given by the Lippmann-Schwinger
equation (e.g. taylor.72 )

$T(\mathbf{k},\mathbf{k}^{\prime},\Omega)=V(\mathbf{k},\mathbf{k}^{\prime})+\sum_{\mathbf{q}}V(\mathbf{k},\mathbf{q})G_{0}^{R}(\mathbf{q},\Omega)T(\mathbf{q},\mathbf{k}^{\prime},\Omega)$
(1)

where $\mathbf{k}$ and $\mathbf{k}^{\prime}$ are the final and initial
relative momenta respectively, $\hbar\Omega$ is the total energy of the
particles, $V(\mathbf{k},\mathbf{k}^{\prime})$ is the interaction, and
$G_{0}^{R}(\mathbf{q},\Omega)$ is the retarded, interaction-free, two-particle
propagator. The interaction is taken to be separable

$V(\mathbf{k},\mathbf{k}^{\prime})=\lambda
u(\mathbf{k})u(\mathbf{k}^{\prime})$ (2)

where $u(\mathbf{k})$ is a square-integrable function, and the interaction is
repulsive (attractive) for $\lambda=1(-1)$. With this interaction, we solve
Eq. (1) analytically in the usual way. Eq. (1) is written as

$T(\mathbf{k},\mathbf{k}^{\prime},\Omega)=\lambda
u(\mathbf{k})\left[u(\mathbf{k}^{\prime})+J(\mathbf{k}^{\prime},\Omega)\right]$
(3)

where

$J(\mathbf{k}^{\prime},\Omega)=\sum_{\mathbf{q}}u(\mathbf{q})G_{0}^{R}(\mathbf{q},\Omega)T(\mathbf{q},\mathbf{k}^{\prime},\Omega)$
(4)

Multiplying both sides of Eq. (3) by
$u(\mathbf{k})G_{0}^{R}(\mathbf{k},\Omega)$ and summing over $\mathbf{k}$
gives

$J(\mathbf{k}^{\prime},\Omega)=\lambda
I(\Omega)\left[u(\mathbf{k}^{\prime})+J(\mathbf{k}^{\prime},\Omega)\right]$
(5)

where

$I(\Omega)=\sum_{\mathbf{k}}u^{2}(\mathbf{k})G_{0}^{R}(\mathbf{k},\Omega)$ (6)

From Eq. (5), $J(\mathbf{k}^{\prime},\Omega)$ is obtained as

$J(\mathbf{k}^{\prime},\Omega)=\frac{\lambda
I(\Omega)u(\mathbf{k}^{\prime})}{1-\lambda I(\Omega)}$ (7)

Substituting Eq. (7) into Eq. (3), we get after a little algebra

$T(\mathbf{k},\mathbf{k}^{\prime},\Omega)=\frac{\lambda
u(\mathbf{k})u(\mathbf{k}^{\prime})}{1-\lambda I(\Omega)}$ (8)

The spectral properties of the two-particle Hamiltonian are related to the
poles and branch cuts of the T-matrix as a function $\Omega$ which is treated
as a complex variable. The energies of the bound states are given by the pole
positions of the T-matrix, which, in our separable model, are the zeros of the
denominator in Eq. (8):

$1-\lambda
I(\Omega)=0\quad\Rightarrow\quad\lambda\sum_{\mathbf{k}}u^{2}(\mathbf{k})G_{0}^{R}(\mathbf{k},\Omega)=1$
(9)

It is known that if the free-particle energy-momentum dispersion is parabolic,
an attractive interaction always supports bound states in two dimensions while
it supports bound states in three dimensions when its strength exceeds a
threshold (see, for example, Refs. simon.76, ; chadan-etal.03, ; yang-de-
llano.89, ; apenko.98, ). We will show that if the free-particle dispersion is
linear, a strength threshold also exists in two dimensions for an attractive
interaction to support bound states.

The retarded Green’s function is

$G_{0}^{R}(\mathbf{k},\Omega)=\frac{1}{\hbar\Omega-\epsilon(\mathbf{k})+i\gamma}$
(10)

where $\epsilon(\mathbf{k})$ is the total free-particle energy of the particle
pair, and $\gamma$ is a loss width. For algebraic convenience, we consider a
simple model for $u(\mathbf{k})$ in which it is a constant over a range of
values of $k=|\mathbf{k}|$ and vanishes elsewhere:
$u(\mathbf{k})=u_{0}\theta(k)\theta(k_{m})$. The pole condition Eq. (9) then
becomes (with the conversion $\sum_{\mathbf{k}}\rightarrow\int
d^{2}k/(2\pi)^{2}$)

$\lambda
I(\Omega)=\lambda\frac{u^{2}_{0}}{2\pi}\int^{k_{m}}_{0}dk\frac{k}{\hbar\Omega-\epsilon(\mathbf{k})+i\gamma}=1$
(11)

We consider the effects of both parabolic and linear dispersions on the bound
state energy.

Parabolic dispersion. The free-particle energy has the form
$\epsilon(\mathbf{k})=ak^{2}$. The particles have momenta $\pm\mathbf{k}/2$.
The constant $a$ is given by $a=\hbar^{2}/(8m_{r})$ where $m_{r}$ is the
reduced mass. Carrying out the integral in Eq. (11) with this
$\epsilon(\mathbf{k})$ reduces the equation to

$-\lambda\frac{u^{2}_{0}}{4\pi a}\ln\left[\frac{\hbar\Omega-
ak^{2}_{m}+i\gamma}{\hbar\Omega+i\gamma}\right]=1$ (12)

We consider the limit $\gamma\downarrow 0$. The natural log function in Eq.
(12) (and hence $I(\Omega)$) has a branch cut along a line segment which lies
infinitesimally below, by $-i\gamma$, the positive real $\Omega$ axis between
$\Omega=0$ and $\Omega=ak^{2}_{m}/\hbar$. We seek solutions to Eq. (12) along
the axis ${\rm Re}\Omega<0$. The log function is continuous along this axis,
and so we set $\gamma=0$. The solution is

$\hbar\Omega=\frac{ak^{2}_{m}}{1-e^{-\frac{4\pi a}{\lambda
u^{2}_{0}}}}\quad,\quad\Omega<0$ (13)

Eq. (13) shows that, as expected, Eq. (12) does not have a solution on the
negative real $\Omega$ axis for $\lambda>0$, and it has one solution for
$\lambda<0$ however weak $u_{0}$ is.

Linear dispersion. In this case, the free-particle energy has the form
$\epsilon(\mathbf{k})=bk$. Using this form, $\lambda I(\Omega)$ becomes

$\displaystyle\lambda I(\Omega)$ $\displaystyle=-\frac{\lambda u^{2}_{0}}{2\pi
b}\int^{k_{m}}_{0}dk\frac{k}{k+z}\quad,\quad
z\equiv-\frac{1}{b}(\hbar\Omega+i\gamma)$ (14) $\displaystyle=-\frac{\lambda
u^{2}_{0}}{2\pi b}\left[k_{m}-z\ln\left(\frac{k_{m}+z}{z}\right)\right]$ (15)

Eq. (15) shows that $I(\Omega)$ has a branch cut which is the line segment
$0\leq\hbar\Omega\leq bk_{m}$ shifted downwards by $-i\gamma$. The limit
$\gamma\downarrow 0$ is again taken. To seek solutions to Eq. (9) on the
negative ${\rm Re}\Omega$ axis, we set $\gamma=0$. On this axis, Eq. (9)
becomes

$\lambda
I(\Omega)=-\lambda\Gamma\left[1-\frac{z}{k_{m}}\ln\left(\frac{k_{m}+z}{z}\right)\right]=1,$
(16)

with $z=-\frac{\hbar\Omega}{b}>0$, where we have defined a strength parameter
$\Gamma\equiv\frac{u^{2}_{0}k_{m}}{2\pi b}$. We have not been able to find an
explicit expression for $z$ from Eq. (16). It is more convenient to use Eq.
(14) to deduce the conditions for existence of solutions for Eq. (16). First,
the form of the integrand in Eq. (14) shows that for $z$ real and positive
(i.e. for $\gamma=0$ and $\Omega$ being real and negative), the integral is
positive. Therefore, for a repulsive interaction $\lambda>0$, $\lambda
I(\Omega)$ is negative, which implies that Eq. (16) does not have a solution.
It is also clear that the integral decreases as $z$ increases. For an
attractive interaction, $\lambda=-1$, $\lambda I(\Omega)$ has a maximum value
of $\Gamma$ at $\Omega=0$ and decreases monotonically to zero as $\Omega$ goes
along the negative real axis to $-\infty$. Therefore the strength threshold
for Eq. (16) to have a solution is $\Gamma=1$. A weak attraction below this
threshold ($\Gamma<1$) does not support a bound state.

###  II.2 Biexcitons in photo-excited 2D transition metal dichalcogenides

In the previous subsection, we have gained some understanding of the effect of
a linear dispersion on the two-body bound state in the single-channel model.
We extend the formal treatment here to the multi-channel case of excitons and
trions in the weakly nonlinear optics of 2D transition metal dichalcogenides
(TMD).

In (direct-gap) semiconductors and TMD, an important driver of $\chi^{(3)}$
nonlinearity is the scattering, including resonant bound-state formation,
among excitons, and possibly also trions. We treat the excitons and trions as
point particles and formulate the scattering using a two-body T-matrix, which
includes in this case the spin/valley degrees of freedom. We consider
specifically $\text{MoSe}_{2}$, which has a relatively simple spin/valley band
structure near the band gap. Only the lowest-lying exciton and trion subbands
are included in our model. The single-particle basis in this subspace is
specified by the parameter set $(\beta,\sigma,\mathbf{p})$, where
$\beta=x\textrm{(for exciton)},t\textrm{(for trion)}$ labels the particle
species, $\sigma=+/-$ is the spin, and $\mathbf{p}$ is the 2D momentum in the
monolayer plane. We limit our considerations to electron-doped samples so that
the trion is made up of an electron-hole pair in one valley, $K$ for example,
and an electron in an inequivalent valley $K^{\prime}$. $\sigma$ for the trion
refers to the spin state of the valence electron corresponding to the hole in
the trion: $\sigma=+(-)$ for valence electron spin being up (down). There is a
one-to-one correspondence between the pair of inequivalent valleys and the
exciton/trion spin states. So no additional parameter labeling the valley in
the single-particle basis is needed.

The exciton-trion Hamiltonian in our model is

$\displaystyle\hat{H}=\sum_{\begin{subarray}{c}\beta=x,t;\mathbf{p},\\\
\sigma,\sigma^{\prime}=\pm\end{subarray}}\langle\sigma|h^{\beta}($
$\displaystyle\mathbf{p})|\sigma^{\prime}\rangle
a^{{\dagger}}_{\beta\mathbf{p}\sigma}a_{\beta\mathbf{p}\sigma^{\prime}}+\sum_{\begin{subarray}{c}\beta_{1},\beta_{2}=x,t;\mathbf{Q},\mathbf{k},\mathbf{k}^{\prime},\\\
\sigma_{1},\sigma_{2},\sigma^{\prime}_{1},\sigma^{\prime}_{2}=\pm\end{subarray}}\langle\sigma_{1}\sigma_{2}|V^{\beta_{1}\beta_{2}}(\mathbf{k},\mathbf{k}^{\prime})|\sigma^{\prime}_{1}\sigma^{\prime}_{2}\rangle$
$\displaystyle\cdot
a^{{\dagger}}_{\beta_{1}(\mathbf{Q}+\mathbf{k}/2)\sigma_{1}}a^{{\dagger}}_{\beta_{2}(\mathbf{Q}-\mathbf{k}/2)\sigma_{2}}a_{\beta_{2}(\mathbf{Q}-\mathbf{k}^{\prime}/2)\sigma^{\prime}_{2}}a_{\beta_{1}(\mathbf{Q}+\mathbf{k}^{\prime}/2)\sigma^{\prime}_{1}}$
(17)

Written as a matrix in the spin sub-basis $(+,-)^{T}$, the one-body
Hamiltonian has the form

$\hat{h}^{\beta}(\mathbf{p})=\begin{pmatrix}\epsilon_{\beta}(p)&J^{inter}_{\beta}(p)e^{-2i\theta_{p}}\\\
J^{inter\ast}_{\beta}(p)e^{2i\theta_{p}}&\epsilon_{\beta}(p)\end{pmatrix}\quad,\quad\mathbf{p}=(p,\theta_{p})$
(18)

Each diagonal term consists of a part derived from the electron-hole band
structure and an intravalley e-h exchange. At the low-momentum limit, the
leading order (in momentum) of the band-structure part is quadratic while that
of the intravalley exchange is linear. The off-diagonal terms flip the spin
and effect intervalley e-h exchange. The spin-flipping comes with a
compensating change in orbital angular momentum (OAM) so that the total
angular momentum of the particle is conserved. For example, in the transition
from $\sigma=-$ to $\sigma=+$, the factor $e^{-2i\theta_{p}}$ reduces the OAM
by 2. The leading order of the function $J^{inter}_{\beta}(p)$ is also linear.
Knowledge of the functional forms of $\epsilon_{\beta}(p)$ and
$J^{inter}_{\beta}(p)$ is not needed for the development in this subsection,
and so we will leave them unspecified. For the trion, since an intervalley
exchange would put the two electrons and the hole in the same valley, and this
configuration does not support a bound state at a comparable energy to the
regular trion, we set the intervalley exchange to zero, $J^{inter}_{t}(p)=0$.
We expect the intravalley e-h exchange occurring inside a trion would affect
the dispersion relation between the trion energy and its center-of-mass
momentum. The e-h pair, however, is a subsystem of the three-body bound state,
and the leading-order behavior of the trion dispersion caused by the exchange
is not clear. Hampered by this uncertainty, we will not consider the effects
of e-h exchange on the bound states containing the trions and will neglect the
intravalley exchange of the trion.

In the interaction term, $\mathbf{k}$ and $\mathbf{k}^{\prime}$ denote the two
particles’ relative momenta in the outgoing and incoming states respectively,
and $\mathbf{Q}$ denotes the total momentum of the two particles. The
interaction is assumed to be independent of $\mathbf{Q}$ and to conserve
individual spins:

$\displaystyle\langle\sigma_{1}\sigma_{2}$ $\displaystyle|$ $\displaystyle
V^{\beta_{1}\beta_{2}}(\mathbf{k},\mathbf{k}^{\prime})|\sigma^{\prime}_{1}\sigma^{\prime}_{2}\rangle=$
(19)
$\displaystyle\delta_{\sigma_{1}\sigma^{\prime}_{1}}\delta_{\sigma_{2}\sigma^{\prime}_{2}}\langle\sigma_{1}\sigma_{2}|V^{\beta_{1}\beta_{2}}(\mathbf{k},\mathbf{k}^{\prime})|\sigma_{1}\sigma_{2}\rangle$

We also write the interaction as a $4\times 4$ matrix in the two-particle spin
sub-basis $(++,+-,-+,--)$, denoting the matrix by
$\hat{V}^{\beta_{1}\beta_{2}}(\mathbf{k},\mathbf{k}^{\prime})$.

The scattering of two excitons/trions is formulated in terms of the T-matrix,
a generic element of which is written as
$\langle\sigma_{1}\sigma_{2}|T^{\beta_{1}\beta_{2}}(\mathbf{k},\mathbf{k}^{\prime},\Omega)|\sigma^{\prime}_{1}\sigma^{\prime}_{2}\rangle$.
The subscripts 1 and 2 refer to the two particles in the scattering. The
symbols for momenta, spins, and particle species are the same as those in the
interaction. We work in the frame of zero total momentum for the two particles
(normal incident beams). $\Omega$ is the total two-particle energy in this
frame. We also use
$\hat{T}^{\beta_{1}\beta_{2}}(\mathbf{k},\mathbf{k}^{\prime},\Omega)$ to
denote the T-matrix written as a $4\times 4$ matrix in the two-particle spin
sub-basis. The Lippmann-Schwinger equation for the T-matrix is

$\displaystyle\hat{T}^{\beta_{1}\beta_{2}}$
$\displaystyle(\mathbf{k},\mathbf{k}^{\prime},\Omega)=\hat{V}^{\beta_{1}\beta_{2}}(\mathbf{k},\mathbf{k}^{\prime})$
$\displaystyle+$
$\displaystyle\sum_{\mathbf{q}}\hat{V}^{\beta_{1}\beta_{2}}(\mathbf{k},\mathbf{q})\hat{G}_{0}^{\beta_{1}\beta_{2}R}(\mathbf{q},\Omega)\hat{T}^{\beta_{1}\beta_{2}}(\mathbf{q},\mathbf{k}^{\prime},\Omega)$

$\hat{G}_{0}^{\beta_{1}\beta_{2}R}(\mathbf{q},\Omega)$ is the two-particle
free retarded Green’s function written as a matrix in the two-particle spin
sub-basis. It is given by

$\hat{G}_{0}^{\beta_{1}\beta_{2}R}(\mathbf{q},\Omega)=\left[\hbar\Omega-\hat{h}^{\beta_{1}}(\mathbf{q}/2)-\hat{h}^{\beta_{2}}(-\mathbf{q}/2)+i\gamma\right]^{-1}$
(21)

In Appendix A, the explicit expression of
$\hat{G}_{0}^{\beta_{1}\beta_{2}R}(\mathbf{q},\Omega)$ is derived from Eq.
(21). The interaction
$\hat{V}^{\beta_{1}\beta_{2}}(\mathbf{k},\mathbf{k}^{\prime})$ is assumed to
conserve the relative orbital angular momentum (OAM) of the two particles. It
can therefore be expanded as a sum over a single OAM integer index:

$\hat{V}^{\beta_{1}\beta_{2}}(\mathbf{k},\mathbf{k}^{\prime})=\sum_{\mu=-\infty}^{\infty}\hat{V}^{\beta_{1}\beta_{2}}_{\mu}(k,k^{\prime})e^{i\mu(\theta_{k}-\theta_{k^{\prime}})}$
(22)

where $\mathbf{k}=(k,\theta_{k})$. The retarded Green’s function and the
T-matrix are likewise expanded,

$\displaystyle\hat{G}^{\beta_{1}\beta_{2}R}_{0}(\mathbf{q},\Omega)$
$\displaystyle=\sum_{\mu=-\infty}^{\infty}\hat{G}^{\beta_{1}\beta_{2}R}_{0\mu}(q,\Omega)e^{i\mu\theta_{q}}\quad,\quad\mathbf{q}=(q,\theta_{q})$
(23)
$\displaystyle\hat{T}^{\beta_{1}\beta_{2}}(\mathbf{k},\mathbf{k}^{\prime},\Omega)$
$\displaystyle=\sum_{\mu,\mu^{\prime}=-\infty}^{\infty}\hat{T}^{\beta_{1}\beta_{2}}_{\mu\mu^{\prime}}(k,k^{\prime},\Omega)e^{i(\mu\theta_{k}-\mu^{\prime}\theta_{k^{\prime}})}$
(24)

Substituting the expansions Eqs. (22)-(24) into Eq. (II.2), we obtain

$\displaystyle\sum_{\mu,\mu^{\prime}=-\infty}^{\infty}\hat{T}^{\beta_{1}\beta_{2}}_{\mu\mu^{\prime}}(k,k^{\prime},\Omega)$
$\displaystyle
e^{i(\mu\theta_{k}-\mu^{\prime}\theta_{k^{\prime}})}=\sum_{\mu=-\infty}^{\infty}\delta_{\mu\mu^{\prime}}\hat{V}^{\beta_{1}\beta_{2}}_{\mu}(k,k^{\prime})e^{i(\mu\theta_{k}-\mu^{\prime}\theta_{k^{\prime}})}$
(25)
$\displaystyle+\sum^{\infty}_{\mu,\mu^{\prime},\mu^{\prime\prime},\mu_{q}=-\infty}$
$\displaystyle\frac{1}{(2\pi)^{2}}\int^{\infty}_{0}dq\int^{2\pi}_{0}d\theta_{q}\hat{V}^{\beta_{1}\beta_{2}}_{\mu}(k,q)e^{i\mu(\theta_{k}-\theta_{q})}$
$\displaystyle\cdot\hat{G}^{\beta_{1}\beta_{2}R}_{0\mu_{q}}(q,\Omega)e^{i\mu_{q}\theta_{q}}\hat{T}^{\beta_{1}\beta_{2}}_{\mu^{\prime\prime}\mu^{\prime}}(q,k^{\prime},\Omega)e^{i(\mu^{\prime\prime}\theta_{q}-\mu^{\prime}\theta_{k^{\prime}})}$

The integral over $\theta_{q}$ gives
$\int^{2\pi}_{0}d\theta_{q}e^{i(-\mu+\mu_{q}+\mu^{\prime\prime})\theta_{q}}=2\pi\delta_{\mu_{q},\mu-\mu^{\prime\prime}}$.
Since Eq. (25) is valid for any values of $\theta_{k}$ and
$\theta_{k^{\prime}}$, it implies the following equation

$\displaystyle\hat{T}^{\beta_{1}\beta_{2}}_{\mu\mu^{\prime}}(k,k^{\prime},\Omega)=\delta_{\mu\mu^{\prime}}\hat{V}^{\beta_{1}\beta_{2}}_{\mu}(k,k^{\prime})+\sum^{\infty}_{\mu^{\prime\prime}=-\infty}\frac{1}{2\pi}\int^{\infty}_{0}dq\,$
$\displaystyle\hat{V}^{\beta_{1}\beta_{2}}_{\mu}(k,q)\hat{G}^{\beta_{1}\beta_{2}R}_{0\mu-\mu^{\prime\prime}}(q,\Omega)$
$\displaystyle\cdot\hat{T}^{\beta_{1}\beta_{2}}_{\mu^{\prime\prime}\mu^{\prime}}(q,k^{\prime},\Omega)$
(26)

The total angular momentum $m_{tot}$ of the two-particle state is the sum of
the spins and the relative OAM, $m_{tot}=\sigma_{1}+\sigma_{2}+\mu$. Since
both $V^{\beta_{1}\beta_{2}}$ and $G^{\beta_{1}\beta_{2}R}_{0}$ conserve
$m_{tot}$, the T-matrix is block-diagonal in $m_{tot}$, which motivates the
switching from the $\mu$ sub-basis to the $m_{tot}$ sub-basis. We define the
interaction and the T-matrix elements in the $m_{tot}$ sub-basis in terms of
those in the $\mu$ sub-basis as (suppressing the subscript $tot$ on $m_{tot}$
for notational simplicity in the remainder of this section)

$\displaystyle\langle\sigma_{1}\sigma_{2}|\mathbb{V}^{\beta_{1}\beta_{2}}_{m}(k,k^{\prime})|\sigma^{\prime}_{1}\sigma^{\prime}_{2}\rangle$
$\displaystyle=\delta_{\sigma_{1}\sigma^{\prime}_{1}}\delta_{\sigma_{2}\sigma^{\prime}_{2}}\langle\sigma_{1}\sigma_{2}|V^{\beta_{1}\beta_{2}}_{m-\sigma_{1}-\sigma_{2}}(k,k^{\prime})|\sigma_{1}\sigma_{2}\rangle$
(27)
$\displaystyle\langle\sigma_{1}\sigma_{2}|\mathbb{T}^{\beta_{1}\beta_{2}}_{m}(k,k^{\prime},\Omega)|\sigma^{\prime}_{1}\sigma^{\prime}_{2}\rangle$
$\displaystyle=\langle\sigma_{1}\sigma_{2}|T^{\beta_{1}\beta_{2}}_{m-\sigma_{1}-\sigma_{2},m-\sigma^{\prime}_{1}-\sigma^{\prime}_{2}}(k,k^{\prime},\Omega)|\sigma^{\prime}_{1}\sigma^{\prime}_{2}\rangle$
(28)

We also introduce the notation

$\langle\sigma_{1}\sigma_{2}|\mathbb{G}^{\beta_{1}\beta_{2}R}_{0}(q,\Omega)|\sigma^{\prime}_{1}\sigma^{\prime}_{2}\rangle=\langle\sigma_{1}\sigma_{2}|G^{\beta_{1}\beta_{2}R}_{0,\sigma^{\prime}_{1}+\sigma^{\prime}_{2}-\sigma_{1}-\sigma_{2}}(q,\Omega)|\sigma^{\prime}_{1}\sigma^{\prime}_{2}\rangle$
(29)

This matrix does not depend on $m$. The $4\times 4$ matrices (in the
$|\sigma_{1}\sigma_{2}\rangle$ sub-basis) of which the above quantities are
elements are denoted by
$\hat{\mathbb{V}}^{\beta_{1}\beta_{2}}_{m}(k,k^{\prime})$,
$\hat{\mathbb{T}}^{\beta_{1}\beta_{2}}_{m}(k,k^{\prime},\Omega)$, and
$\hat{\mathbb{G}}^{\beta_{1}\beta_{2}R}_{0}(q,\Omega)$. In terms of these
matrices, the Lippmann-Schwinger equation (26) is rewritten in a compact form
as

$\displaystyle\hat{\mathbb{T}}^{\beta_{1}\beta_{2}}_{m}(k,k^{\prime},\Omega)=\hat{\mathbb{V}}^{\beta_{1}\beta_{2}}_{m}(k,k^{\prime})+\frac{1}{2\pi}\int^{\infty}_{0}dq\,$
$\displaystyle\hat{\mathbb{V}}^{\beta_{1}\beta_{2}}_{m}(k,q)\hat{\mathbb{G}}^{\beta_{1}\beta_{2}R}_{0}(q,\Omega)$
$\displaystyle\cdot\hat{\mathbb{T}}^{\beta_{1}\beta_{2}}_{m}(q,k^{\prime},\Omega)$
(30)

For the momentum dependence of the interaction, we adopt a separable form in
each OAM and spin channel

$\langle\sigma_{1}\sigma_{2}|V^{\beta_{1}\beta_{2}}_{\mu}(k,k^{\prime})|\sigma_{1}\sigma_{2}\rangle=\lambda^{\beta_{1}\beta_{2}}_{\mu\sigma_{1}\sigma_{2}}u^{\beta_{1}\beta_{2}}_{\mu\sigma_{1}\sigma_{2}}(k)u^{\beta_{1}\beta_{2}}_{\mu\sigma_{1}\sigma_{2}}(k^{\prime})$
(31)

with $\lambda^{\beta_{1}\beta_{2}}_{\mu\sigma_{1}\sigma_{2}}=1$ or $-1$. For
formal manipulations, $\lambda^{\beta_{1}\beta_{2}}_{\mu\sigma_{1}\sigma_{2}}$
and $u^{\beta_{1}\beta_{2}}_{\mu\sigma_{1}\sigma_{2}}(k)$ are also written in
$m$-sub-basis matrix form

$\displaystyle\langle\sigma_{1}\sigma_{2}|\hat{\Lambda}^{\beta_{1}\beta_{2}}_{m}|\sigma^{\prime}_{1}\sigma^{\prime}_{2}\rangle$
$\displaystyle=\delta_{\sigma_{1}\sigma^{\prime}_{1}}\delta_{\sigma_{2}\sigma^{\prime}_{2}}\lambda^{\beta_{1}\beta_{2}}_{m-\sigma_{1}-\sigma_{2},\sigma_{1}\sigma_{2}}$
(32)
$\displaystyle\langle\sigma_{1}\sigma_{2}|\hat{\mathbb{U}}^{\beta_{1}\beta_{2}}_{m}(k)|\sigma^{\prime}_{1}\sigma^{\prime}_{2}\rangle$
$\displaystyle=\delta_{\sigma_{1}\sigma^{\prime}_{1}}\delta_{\sigma_{2}\sigma^{\prime}_{2}}u^{\beta_{1}\beta_{2}}_{m-\sigma_{1}-\sigma_{2},\sigma_{1}\sigma_{2}}(k)$
(33)

in terms of which Eq. (31) can be written as

$\hat{\mathbb{V}}^{\beta_{1}\beta_{2}}_{m}(k,k^{\prime})=\hat{\Lambda}^{\beta_{1}\beta_{2}}_{m}\hat{\mathbb{U}}^{\beta_{1}\beta_{2}}_{m}(k)\hat{\mathbb{U}}^{\beta_{1}\beta_{2}}_{m}(k^{\prime})$
(34)

With this interaction, the T-matrix equation, Eq. (30) can be formally solved
in a similar way as its single-channel counterpart Eq. (1), proper care being
taken of the non-commutativity of (most of) the matrices involved. The
solution is

$\hat{\mathbb{T}}^{\beta_{1}\beta_{2}}_{m}(k,k^{\prime},\Omega)=\hat{\Lambda}^{\beta_{1}\beta_{2}}_{m}\hat{\mathbb{U}}^{\beta_{1}\beta_{2}}_{m}(k)\left[1-\hat{\mathbb{I}}^{\beta_{1}\beta_{2}}_{m}(\Omega)\right]^{-1}\hat{\mathbb{U}}^{\beta_{1}\beta_{2}}_{m}(k^{\prime})$
(35)

where

$\hat{\mathbb{I}}^{\beta_{1}\beta_{2}}_{m}(\Omega)=\frac{1}{2\pi}\int^{\infty}_{0}dkk\hat{\mathbb{U}}^{\beta_{1}\beta_{2}}_{m}(k)\hat{\mathbb{G}}^{\beta_{1}\beta_{2}R}_{0}(k,\Omega)\hat{\Lambda}^{\beta_{1}\beta_{2}}_{m}\hat{\mathbb{U}}^{\beta_{1}\beta_{2}}_{m}(k)$
(36)

The trion being a fermion, the two-trion scattering and bound-state
wavefunctions observe total antisymmetry under exchange. When the spin state
is symmetric (antisymmetric), the spatial orbital wavefunction is a
superposition of odd (even) OAM states. In our formalism, this antisymmetry
condition is imposed by setting the even-OAM components of the trion-trion
interaction in the parallel-spin channel to zero:

$\langle\sigma\sigma|V^{tt}_{\mu}(k,k^{\prime})|\sigma\sigma\rangle=0\quad\mathrm{when}\,\mu\,\mathrm{is\,even}$
(37)

###  II.3 Differential tranmission

In this subsection, we relate the $\text{MoSe}_{2}$ exciton/trion T-matrix to
the $\chi^{(3)}$ pump-probe differential transmission (DT). We consider
normally incident pump and probe on the monolayer. (One of the beams may be
slightly oblique to introduce a directional separation. We will ignore the
effect of the small deviation from the normal.) Up to $\chi^{(3)}$, the
transmission of the probe is given in the frequency domain by

$\displaystyle T_{i}(\omega)$
$\displaystyle=\frac{|E^{(t)}_{pr,i}(\omega)|^{2}}{|E^{(inc)}_{pr,i}(\omega)|^{2}}$
$\displaystyle=1-\frac{4\pi\omega}{nc}\frac{\mathrm{Im}\left[E^{(inc)\ast}_{pr,i}(\omega)\left(P^{(1)}_{pr,i}(\omega)+P^{(3)}_{pr,i}(\omega)\right)\right]}{|E^{(inc)}_{pr,i}(\omega)|^{2}}$
(38)

The subscript $i=+,-$ labels the circular polarized components of the various
quantities, and the subscript $pr$ means the probe. $E^{(inc)}_{pr,i}$ and
$E^{(t)}_{pr,i}$ are the incident and transmitted probed fields respectively.
$P^{(k)}_{pr,i}$ is the $\chi^{(k)}$ induced polarization. $n$ is the
refractive index of the medium on the two sides of the monolayer. Through the
induced polarization, the linear ($\chi^{(1)}$) transmission and $\chi^{(3)}$
differential transmission are related to the susceptibility by

$\displaystyle T^{(1)}_{i}(\omega)$
$\displaystyle=1-\frac{4\pi\omega}{nc}\sum_{j=+,-}\mathrm{Im}\left[\chi^{(1)}_{ij}(\omega)\frac{E^{(t)\ast}_{pr,i}(\omega)E^{(t)}_{pr,j}(\omega)}{|E^{(inc)}_{pr,i}(\omega)|^{2}}\right]$
(39) $\displaystyle DT^{(3)}_{i}(\omega)$
$\displaystyle=-\frac{4\pi\omega}{nc}\sum_{j,k,\ell=+,-}\mathrm{Im}\left[\frac{1}{(2\pi)^{2}E^{(inc)}_{pr,i}(\omega)}\int
d\omega_{1}d\omega_{2}d\omega_{3}\right.$ (40)
$\displaystyle\quad\quad\quad\cdot\delta(\omega-\omega_{1}-\omega_{2}+\omega_{3})\chi^{(3)}_{ijk\ell}(\omega_{1},\omega_{2},\omega_{3})$
$\displaystyle\quad\quad\quad\left.\cdot\left(E^{(t)}_{pr,j}(\omega_{1})E^{(t)}_{p,k}(\omega_{2})E^{(t)\ast}_{p,\ell}(\omega_{3})+E^{(t)}_{p,j}(\omega_{1})E^{(t)}_{pr,k}(\omega_{2})E^{(t)\ast}_{p,\ell}(\omega_{3})\right)\right]$

where the subscript $p$ means the pump. We specialize to monochromatic fields
with pump frequency $\omega_{p}$ and probe frequency $\omega_{pr}$:

$E^{(inc)}_{p,i}(\omega)=2\pi\delta(\omega-\omega_{p})E_{p,i}\quad,\quad
E^{(inc)}_{pr,i}(\omega)=2\pi\delta(\omega-\omega_{pr})E_{pr,i}$ (41)

Assuming weak response, we approximate the transmitted fields by the incident
fields on the right hand sides of Eqs. (39) and (40). With these
simplifications, the transmission equations become

$\displaystyle T^{(1)}_{i}(\omega_{pr})$
$\displaystyle=1-\frac{4\pi\omega}{nc}\sum_{j=+,-}\mathrm{Im}\left[\chi^{(1)}_{ij}(\omega_{pr})\frac{E^{\ast}_{pr,i}E_{pr,j}}{|E_{pr,i}|^{2}}\right]$
(42) $\displaystyle DT^{(3)}_{i}(\omega)$
$\displaystyle=-\frac{4\pi\omega}{nc}\sum_{j,k,\ell=+,-}\mathrm{Im}\left[\frac{E^{\ast}_{pr,i}E_{pr,j}}{|E_{pr,i}|^{2}}\left(\chi^{(3)}_{ijk\ell}(\omega_{pr},\omega_{p},\omega_{p})+\chi^{(3)}_{ikj\ell}(\omega_{p},\omega_{pr},\omega_{p})\right)\right.$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\left.\cdot
E_{pr,k}E^{\ast}_{pr,\ell}\right]$ (43)

In our model, the linear susceptibility is the sum of contributions from the
exciton and the negative trion:

$\chi^{(1)}_{ij}(\omega_{pr})=\delta_{ij}\sum_{\beta=x,t}\chi^{\beta(1)}_{i}(\omega_{pr})$
(44)

where

$\chi^{x(1)}_{i}(\omega_{pr})=-\frac{|\mu_{i}|^{2}}{\hbar\omega_{pr}-\epsilon_{x}(0)+i\gamma_{x}}\quad,\quad\chi^{t(1)}_{i}(\omega_{pr})=-\frac{|f_{i}|^{2}n^{\prime}_{e,i}}{\hbar\omega_{pr}-\epsilon_{t}(0)+i\gamma_{t}}$
(45)

$\mu_{i}$ is the exciton electric interband dipole moment, $f_{i}$ is the
trion excitation amplitude, and $n^{\prime}_{e,i}$ is the doped electron
density in a valley inequivalent to that where the other two constituent
particles, an e-h pair, of the trion reside. The two-particle scattering
contribution to the $\chi^{(3)}$ susceptibility is given by

$\displaystyle\chi^{(3)}_{ijk\ell}(\omega_{pr},\omega_{p},\omega_{p})=-\sum_{\beta_{1},\beta_{2}=x,t}$
$\displaystyle\left(\frac{\chi^{\beta_{1}(1)}_{i}(\omega_{pr})}{\zeta_{\beta_{1},i}}\right)\left(\frac{\chi^{\beta_{1}(1)}_{j}(\omega_{pr})}{\zeta^{\ast}_{\beta_{1},j}}\right)\left(\frac{\chi^{\beta_{2}(1)}_{k}(\omega_{p})}{\zeta^{\ast}_{\beta_{2},k}}\right)$
(46)
$\displaystyle\cdot\left(\frac{\chi^{\beta_{2}(1)\ast}_{\ell}(\omega_{p})}{\zeta_{\beta_{2},\ell}}\right)\langle
i\ell|T^{\beta_{1}\beta_{2}}_{00}(0,0,\omega_{pr}+\omega_{p}|jk\rangle$

where

$\zeta_{x,i}=\mu_{i}\quad,\quad\zeta_{t,i}=f_{i}\sqrt{n^{\prime}_{e,i}}$ (47)

##  III Numerical results and discussion

In this section, we present numerical results of the theory developed in Sec.
II.

Unless otherwise noted, we use the following parameters. The effective
electron mass follows from a parabolic approximation of a hyperbolic band
model, and is given by $m_{e}=\hbar^{2}E_{g}/(2a_{L}^{2}t_{h}^{2})$, where we
use $a_{L}=3.313$ Å for the lattice constant, $t_{h}=0.94$eV for the hopping
matrix element, and $E_{g}=2.126$eV for the bandgap xiao-etal.12 ; yu-etal.14
. This yields $m_{e}=0.83m_{0}$, where $m_{0}$ is the electron mass in vacuum,
which is consistent with Ref. larentis-etal.2018, . We assume the hole and
electron mass to be equal, so that the exciton mass is $m_{x}=2m_{e}$, and the
trion mass is $m_{trion}=3m_{e}$. We use the Coulomb potential in the MoSe2
layer developed in Ref. vantuan-etal.2018, , Eq. (5), which can be viewed as
an improved version of the Rytova-Keldysh potential. Specifically, we use the
dielectric function given as Eq. (6)-(7) in Ref. vantuan-etal.2018, , with the
monolayer thickness $d=0.6$nm and $\ell_{+}=\ell_{-}=5d$. We call this
dielectric function $\epsilon_{D}(k)$. Hence, using Gaussian units, the
Coulomb interaction is $V^{c}(k)=\frac{2\pi e^{2}}{k\epsilon_{D}(k)}$. We use
the separable model Eq. (31) for the two-particle ($xx$, $xt$, $tt$)
interactions. Each component is assumed to be of constant strength over a
common range of relative momentum:

$\displaystyle u^{\beta_{1}\beta_{2}}_{\mu\sigma_{1}\sigma_{2}}(k)$
$\displaystyle=\tilde{u}^{\beta_{1}\beta_{2}}_{\mu\sigma_{1}\sigma_{2}}\quad,\quad
0\leq k\leq k_{max}$ (48) $\displaystyle=0\quad\quad\quad,\quad{\rm
otherwise}$

where the $\tilde{u}^{\beta_{1}\beta_{2}}_{\mu\sigma_{1}\sigma_{2}}$’s are
strength constants, which we set to be dependent on the particle species but
not on the spin states. We chose the strength of the exciton-exciton
interaction to be
$\left(\tilde{u}^{xx}_{\mu\sigma_{1}\sigma_{2}}\right)^{2}=u_{0}^{2}=6$ eV Å2
and the upper limit of the wave vector range, $k_{max}=6/a_{B}$, where the 2D
Bohr radius is $a_{B}=1$ nm. The exciton-trion and trion-trion interactions
are set at
$\left(\tilde{u}^{xt}_{\mu\sigma_{1}\sigma_{2}}\right)^{2}=0.6u_{0}^{2}$,
$\left(\tilde{u}^{tt}_{\mu\sigma_{1}\sigma_{2}}\right)^{2}=0.2u_{0}^{2}$,
respectively. In each channel, the parameter
$\lambda^{\beta_{1}\beta_{2}}_{\mu\sigma_{1}\sigma_{2}}=\pm 1$ determines
whether the interaction is repulsive or attractive. The trion is set to be 26
meV below the exciton resonance.

The interband dipole matrix element entering the first-order exciton
susceptibility, $er_{cv}$ (where $e$ is the absolute value of the electron
charge in vacuum), is determined according to Ref. yu-etal.14,  as
$r_{cv}=\sqrt{2}a_{L}t_{h}/E_{g}$. For simplicity, we use the hydrogen-like
form for the exciton wave function at zero electron-hole separation,
$\phi_{1s}(r=0)=\sqrt{8/\pi}\ a_{B}$. The oscillator strength of the exciton
is proportional to the square of the exciton dipole moment,
$\mu_{i}=er_{cv}\phi_{1s}(r=0)$, where $i=+,-$ indicates the exciton spin, and
the linear susceptibility is independent of spin. For clarity, we take the
oscillator strength of the trion to be relatively large, 20% of that of the
exciton, i.e. $|f_{i}|^{2}n_{ei}^{\prime}=0.2|\mu_{i}|^{2}$. In actual
simulations of typical samples the trion oscillator strength would typically
be chosen to be much smaller, for example 1%. The dephasing of the exciton is
taken to be $\gamma_{0}\equiv\gamma_{exciton}=2$meV, and we take the trion
dephasing to be $\gamma_{-}\equiv\gamma_{trion}=0.5\gamma_{exciton}$, which is
motivated by the fact that trions are more likely to be localized and
localized exciations have usually smaller dephasing rates than delocalized
ones. Finally, unless otherwise noted, the following results assume orbital
angular momentum (OAM) filtering, by which we mean that both the incoming and
outgoing waves in the T-matrix have only zero OAM.

Figure 1:  (Color online.) (a) Exciton (upper curve) and trion (lower curve)
dispersion without e-h exchange. (b) Exciton dispersion with intravalley
exchange. The two exciton-spin states are degenerate and the two dispersions
cannot be distinguished. (c) Exciton dispersion with intervalley exchange,
shown are the upper and lower dispersion branches. (d) Exciton dispersion
(upper and lower branches) with intravalley and intervalley exchange.

The single-particle Hamiltonian is given in Eq. (18). For the exciton, the
intervalley e-h exchange is given by yu-etal.14

$J^{inter}_{x}(k)=3V^{c}(k)\frac{1}{2}r_{cv}^{2}k^{2}|\phi_{1s}(r=0)|^{2}$
(49)

which is real-valued in our case. Here, $k$ is the magnitude of the exciton
wave vector k, the factor of 3 comes from the fact that there are three
equivalent valleys, and the factor 1/2 comes from the factor of $1/\sqrt{2}$
in the matrix element of the lattice-periodic parts of the Bloch wave
functions, $\langle
u_{c\textbf{k}}|i\vec{\nabla}|u_{c\textbf{k}}\rangle=r_{cv}\frac{1}{\sqrt{2}}(\hat{x}+i\hat{y})$.
As mentioned above, the phase factors that depend on the direction of the wave
vector in Eq. (18) show the spin-orbit coupling: flipping the exciton spin
from plus (minus) to minus (plus) is associated with an increase (decrease) of
the orbital angular momentum of 2 (in units of $\hbar$).

In contrast to the intervalley exchange, the intravalley exchange does not
couple inequivalent Dirac valleys and hence does not couple opposite exciton
spins. Therefore, the intravalley exchange modifies only the diagonal terms,
in the Hamiltonian Eq. (18) which is given as

$\epsilon_{x}(k)=\hbar^{2}k^{2}/(2m_{x})+J^{intra}_{x}(k)$ (50)

where the intravalley exchange is taken to be the same as the intervalley
exchange $J^{intra}_{x}(k)=J^{inter}_{x}(k)$.

Figure 2:  (Color online.) (a) Differential transmission (DT) in the cross-
circular configuration without intervalley e-h exchange. The intravalley
exchange is increased equidistantly with $\eta_{intra}$=0, 0.2, 0.4, 0.6, 0.8,
1 from top to bottom. Except for the bottom curve ($\eta$=1), the spectra are
shifted vertically for clarity. (b) Same data as in (a) but zoomed in to the
biexciton dip, indicated by vertical arrows. The blue shift of the biexciton
with increasing intravelley exchange can clearly be seen. The spectra in (b)
are not vertically shifted.  Figure 3:  (Color online.) (a) The 2-exciton
T-matrix component $T_{-+-+}$ corresponding to the spectra in Fig. 2. The
exchange interaction is varied in the same way as in Fig. 2. The biexciton
peak can clearly be seen. (b) Same as (a) but for $T_{-++-}$. Here, the
T-matrix component is zero independent of the e-h exchange strength. The
spectra are again shifted vertically for clarity.

Figure 1 summarizes the effects of the e-h exchange on the dispersion
relations, which are in agreement with Refs. yu-etal.14, ; wu-etal.15prb, ;
steinhoff-etal.2018, ; schneider-etal.2019, . Note that, because of the
k-dependent dielectric constant $\epsilon_{D}(k)$, the Coulomb potential is
not simply proportional to $1/k$, as would be the case if $\epsilon_{D}(k)$
were k-independent, and therefore the exchange interaction is not simply
proportional to $k$. This becomes apparent when we look at the exciton
dispersion, Fig. 1. Figure 1a shows the bare exciton dispersion without
exchange effect, $\varepsilon(k)=\hbar^{2}k^{2}/(2m_{x})$. Figure 1b shows the
exciton dispersion for the case where we only have intravalley e-h exchange.
In this case, the two spin states are degenerate, and given by
$\varepsilon(k)=\hbar^{2}k^{2}/(2m_{x})+J^{intra}_{x}(k)$. We see that, at
small wave vectors, the dispersion is approximately linear, corresponding to
linear-dispersion model discussed above (Sec. II.1). However, due to the
k-dependence of the dielectric function, the dispersion starts to deviate
significantly from the linear form at approximately $k\approx 10^{6}$cm-1.
Figure 1c shows the exciton dispersion for the case where we only have
intervalley e-h exchange. In this case, the degeneracy is lifted, and we have
two branches similar to Ref. yu-etal.14, , with the lower branch showing a
minimum at non-zero wave vector,
$\varepsilon(k)_{\pm}=\hbar^{2}k^{2}/(2m_{x})\pm J^{inter}_{x}(k)$. Finally,
Fig. 1d shows the exciton dispersion for the case where we have both, inter
and intravalley e-h exchange. In this case, the lower branch is again the un-
modified exciton dispersion,
$\varepsilon(k)_{-}=\hbar^{2}k^{2}/(2m_{x})+J^{intra}_{x}(k)-J^{inter}_{x}(k)=\hbar^{2}k^{2}/(2m_{x})$,
while the upper branch contains the exchange interaction twice,
$\varepsilon(k)_{+}=\hbar^{2}k^{2}/(2m_{x})+J^{intra}_{x}(k)+J^{inter}_{x}(k)=\hbar^{2}k^{2}/(2m_{x})+2J(k)$
where $J^{intra}_{x}(k)=J^{inter}_{x}(k)\equiv J(k)$ In the following
numerical analysis, we assume that we can adjust the strength of the inter and
intravalley exchange independently,
$J^{inter}_{x}(k)\rightarrow\eta_{inter}J^{inter}_{x}(k)$ and
$J^{intra}_{x}(k)\rightarrow\eta_{intra}J^{intra}_{x}(k)$, and we vary $\eta$
between 0 and 1.

Figure 4:  (Color online.) Differential transmission (DT) in the cross-
circular configuration without intravalley e-h exchange. The intervalley
exchange is increased equidistantly with $\eta_{inter}$=0, 0.2, 0.4, 0.6, 0.8,
1 from top to bottom. Except for the bottom curve ($\eta$=1), the spectra are
shifted vertically for clarity. For clarity, the biexciton resonance is
indicated by the dash-dotted line. Figure 5:  (Color online.) (a) The
2-exciton T-matrix component $T_{-+-+}$ corresponding to the spectra in Fig.
4. The exchange interaction is varied in the same way as in Fig. 4. The
biexciton peak can clearly be seen as it splits and shifts with increasing e-h
exchange coupling. (b) Same as (a) but for $T_{-++-}$. The spectra (a) and (b)
are shifted vertically for clarity.  Figure 6:  (Color online.) Same data as
in Fig. 5, but only showing the two spectra at full exchange strength (purple
curves) and for comparion the curve for zero exchange strength. The split
biexciton corresponds to the two disperion curves shown in Fig. 1c. The sum of
the two dispersions, entering the T-matrix, is equivalent to the sum of two
exciton dispersion without e-h exchange and yield the upper of the split
biexciton peaks, which coincides with the biexction without exchange.

We begin the analysis for the simplest case, in which we have no trions, and
we have only intravalley exchange. This reduces the system to the scalar case
discussed in Sec. II.1. Figure 2 shows the differential transmission (DT) in
the cross-circular configuration, with the pump being $\sigma_{+}$ polarized
and spectrally centered at the exciton, and the probe $\sigma_{-}$ and its
frequency is scanned in the figure. We see that without exchange, the DT shows
a feature similar to a blue shift: positive DT just below the exciton,
negative DT just above it; however in this case the negative and positive
parts are not equal in height, thus deviating from an exact blue shift signal,
and a biexciton dip (induced absorption) at approximately 20 meV below the
exciton (i.e. the biexciton binding energy is 20meV). This is in qualitative
agreement with the case of GaAs, which has been shown in Ref. sieh-etal.99,
(see Fig. 4d in Ref. sieh-etal.99,  and note that the induced biexciton
absorption peak (corresponding to a dip in DT) in GaAs is not well separated
from the exciton feature).

When we increase the e-h exchange from zero to its full strength, we see the
biexciton binding energy decreasing and the biexciton dip eventually merges
with the features in the vicinity of the exciton, making the exciton DT line
shift more similar to a pump-induced red shift (this is not an exact red
shift, because the DT signal at the exciton is, from left to right, negative-
positive-negative, instead of negative-positive in the exact redshift). In
order to gain more insight into this behavior, we show the underlying
2-exciton T-matrix in Fig. 3. For the discussion of results shown in Figs. 2
to 9, we simplify the notation and write
$\langle\sigma_{1}\sigma_{2}|T^{xx}_{00}|\sigma^{\prime}_{1}\sigma^{\prime}_{2}\rangle$
as $T_{\sigma_{1}\sigma_{2}\sigma^{\prime}_{1}\sigma^{\prime}_{2}}$. We note
that in general (Eqs. (43) and (46)), DT in the cross-circular configuration
receives contributions from $T_{-+-+}$ (i.e. the T-matrix where there is no
net change in the spin of either exciton) and $T_{-++-}$ (the T-matrix where
the two excitons flip their spins, one from -1 to 1 and the other from +1 to
-1, so that the total spin of the two excitons remains unchanged). In
$T_{-+-+}$, as expected, the biexciton peak moves from -20 meV in the case
without e-h exchange to a frequency slightly above the two-exciton continuum
edge in the case of full exchange, Fig. 3a. This corroborates our finding in
Sec. II.1 that a linear exciton dispersion decreases the bound state binding
energy and can prohibit biexcitons from forming. Since in the case of only
intravalley e-h exchange we have no mixing of the exciton spin states, the
T-matrix in the ’-++-” channel is zero, regardless of the strength of the e-h
intravalley exchange.

Next, we show results where we have only intervalley e-h exchange but no
intravalley exchange. Figure 4 shows the evolution of the cross-circular DT in
increasing intervalley exchange. Without exchange, we have again the blue-
shift (or close to blue shift) feature at the exciton, and the biexciton dip
at -20meV. At full intervalley exchange, the exciton feature has developed a
slightly negative feature just below the exciton, and the negative feature
just above the exciton has slightly diminished. The biexciton dip has moved to
lower frequencies (i.e. the biexciton binding energy has increased). Again, to
obtain more insight into this behavior, we look at the T-matrix. Figure 5a
shows, as before, a prominent biexciton peak in the case without exchange.
Since the intervalley exchange couples the exciton spins and leads to a
splitting in the exciton dispersion, Fig. 1c, the biexciton peak in $T_{-+-+}$
splits into two when the intervalley exchange is non-zero. However, there is
no such splitting in the DT spectrum. In order to understand this, we plot in
Fig. 5b the T-matrix in the ’-++-” channel. When intervalley exchange is
switched on, $T_{-++-}$ becomes nonzero, and moreover the spin-mixing is
accompanied by spin-orbit coupling. In this channel we find that the two
biexciton resonances correspond to different signs in the T-matrix. The usual
case, where Im(T) is negative, corresponds to a biexciton resonance that leads
to induced absorption (or a dip in DT). In the present case, the upper of the
two split biexciton states has the opposite sign, which leads to a
cancellation against the upper peak in $T_{-+-+}$ and the absence of the
second biexciton peak in DT. Moreover, we see from Fig. 6 that the upper of
the two biexciton peaks is exactly at the biexciton resonance that we have in
the absence of the intervalley exchange. This can be understood from the fact
that if the two excitons that participate in the scattering are from different
branches, one from the upper and one from the lower, then the sum of the two
dispersion relations yields exactly the dispersions of two excitons without
modification due to exchange, Fig. 1c.

In Fig. 7 we finally switch on both exchange contributions. We start from the
case of only intervalley exchange and increase the intravalley exchange from
zero to full strength. In this case, the biexciton dip moves from -30 meV
towards the exciton and there merges with the exciton scattering feature,
resulting in a negative-positive-(weakly)negative DT lineshape. The
corresponding T-matrix results, Fig. 8, shows the biexciton, as it moves
toward the two-exciton continuum, keeps its splitting such that the upper of
the two biexciton peaks remains cancelled in the DT regardless of the strength
of the intravalley exchange.

Figure 7:  (Color online.) Differential transmission (DT) in the cross-
circular configuration ranging from only intervalley exchange to full
intravalley and intervalley exchange, $\eta_{inter}$=1 and $\eta_{intra}$=0,
0.2, 0.4, 0.6, 0.8, 1 from top to bottom. Except for the bottom curve
($\eta$=1), the spectra are shifted vertically for clarity. (b) Same data as
in (a) but zoomed in to the biexciton dip, indicated by vertical arrows. The
blue shift of the biexciton with increasing intravelley exchange can clearly
be seen. The spectra in (b) are not vertically shifted.  Figure 8:  (Color
online.) (a) The 2-exciton T-matrix component $T_{-+-+}$ corresponding to the
spectra in Fig. 7. The exchange interaction is varied in the same way as in
Fig. 7. The biexciton peak can clearly be seen as it splits and shifts with
increasing e-h exchange coupling. (b) Same as (a) but for $T_{-++-}$. The
spectra (a) and (b) are shifted vertically for clarity.

The fact that the e-h intervalley exchange interaction leads to spin coupling
suggests that a biexciton resonance could be visible in the DT signal in the
co-circular configuration. Without e-h intervalley exchange interaction, the
DT in the co-circular configuration shows a exciton blue shift, which is a
consequence of the repulsive exciton-exciton interaction in the spin-triplet
channel, and theoretically is dominated by the Hartree-Fock (first order in
the exciton-exciton interaction) contribution to the T-matrix $T_{++++}$. When
the e-h intervalley exchange interaction is included in the calculation,
$T_{++++}$ is coupled to $T_{-+-+}$, the latter possibly having a pronounced
biexciton resonance. Figure 9 shows the corresponding T-matrix components and
DT spectra for our parameters are full exchange interaction strength. We see
that $T_{++++}$ indeed exhibits a biexciton resonance at the same energy as
the biexciton resonance in $T_{-+-+}$. This is in contrast to the case of
GaAs, where the co-circular $T_{++++}$ has only the onset of a two-exciton
continuum, see Ref. takayama-etal.02, . However, the biexciton resonance in
the co-circular channel does not flip the sign of the exciton shift; Fig. 9c
still shows a blue shift in the co-circular configuration. The reason for this
can be seen from Fig. 9b, which shows that the real part of $T_{++++}$ is much
larger than the imaginary part, and it is the real part that is responsible
for the exciton blue shift.

Figure 9:  (Color online.) (a) Real (magenta curve) and imaginary (green
curve) part of the 2-exciton T-matrix component $T_{-+-+}$ that dominantly
determines the DT in the cross-circular configuration. The imaginary part is
the same as in Fig. 8. The biexciton resonance is approximately at -3.5 meV.
(b) Same as (a) but for $T_{++++}$. The imaginary part shows a biexciton
resonance which is a consequence of the spin coupling due to e-h exchange and
is not present without e-h exchange, and the 2-exciton continuum that would be
present even without e-h exchange. Note that the T-matrix in (a) is much
larger than in (b). (c) Differential transmission in the cross-circular
configuration, dash-dotted red line, which is the same as in Fig. 7, and the
co-circular configuration (solid blue line). All results in this figure are at
full e-h exchange strength.

We now evaluate the full theory developed in Sec. II.2. We add the trion and
show results where we include the T-matrices describing exciton-exciton,
exciton-trion, and trion-trion coupling. In Fig. 10, we show the cross-
circular DT under simultaneous increase of intervalley and intravalley
exchanges. In other words, for each curve the two exchange energies are
multiplied by the same factor of $\eta=\eta_{inter}=\eta_{intra}$. At zero
exchange, we have again the (almost ideal) blue shift signal at the exciton
(i.e. positive-negative), a relatively shallow biexciton dip at -20 meV, and
an almost ideal bleaching signal at the trion (mostly positive DT). At full
exchange strength, we find, as above, the exciton to have the negative-
positive-(weakly)negative signal, and the trion develops an almost ideal red
shift. To analyse this further, we show in Fig. 11 the relevant T-matrix
components. The evolution of the exciton-exciton T-matrix is similar to the
case without trions discussed above. In the absence of exchange, the exciton-
trion T-matrix shows a peak, corresponding to a bound state of an exciton and
a trion, at approximately -30 meV. This peaks shifts to -25meV and undergoes a
strong lineshape deformation as the exchange is increased to full strength.
The superposition of the trion-exciton bound state at -25 meV and the trion
resonance at -26 meV then leads to the development of the blue shift signal at
the trion frequency seen in Fig. 10. For clarity, we show in Fig. 12 a
comparison of the relevant T-matrix components for the case of full exchange.
The exciton-exciton interaction determines only the signal in the vicinity of
the exciton, there is no biexciton peak. The trion-exciton bound state is
represented by a strongly asymmetric peak at -25 meV, and the trion-trion
bound state gives an peak at -63meV. In DT, there is no feature visible at the
trion-trion bound state frequency. However, there is in principle a dip in DT
at that frequency. We have verified numerically that an artificial enhancement
of the trion-trion oscillator strength yields a clearly visible dip in DT at
-63 meV, but for our normal parameters this dip is too small to be visible.

Figure 10:  (Color online.) Differential transmission (DT) in the cross-
circular configuration ranging from zero intravalley and intervalley exchange
to full intravalley and intervalley exchange, $\eta_{inter}=\eta_{intra}$=0,
0.2, 0.4, 0.6, 0.8, 1 from top to bottom. Except for the bottom curve
($\eta$=1), the spectra are shifted vertically for clarity. In contrast to the
previous DT figures, here the trion is taken into account. The DT response of
the trion can clearly be seen at approximately 26 meV below the exciton.
Effects of e-h exchange in the trion is omitted.  Figure 11:  (Color online.)
(a) The 2-exciton T-matrix component $T_{-+-+}$ corresponding to the spectra
in Fig. 10. The exchange interaction is varied in the same way as in Fig. 10.
The biexciton peak can clearly be seen as it splits and shifts with increasing
e-h exchange coupling. (b) Same as (a) but for $T_{-++-}$. (c) Same as in (a)
but for the exciton-trion T-matrix component. The spectra (a), (b) and (c) are
shifted vertically for clarity.

The results including trions discussed so far are without exchange effects in
the trions. As discussed in Sec. II.2, we omit intervalley e-h exchange
because in a trion that involves the electron from the valley where the
exciton resides, and an electron of opposite spin in a non-equivalent valley
(2-electron singlet state), e-h exchange would transfer the exciton to the
non-equivalent valley, leading to a state with two electrons of the same spin
(2-electron spin triplet state) and possible disintegration of the 3-particle
state into an exciton and a free electron. We have also omitted intravalley
e-h exchange of the exciton that is part of the trion. Regarding the
intravelley e-h exchange effect for the trion, we note that Ref. bayer-
etal.2002  pointed out that the e-h exchange interaction vanishes in charged
excitons in III-V quantum dots. It might therefore be possible that the
intravalley e-h exchange in trions in TMDs is also small. Nevertheless, we
have also evaluated a purely hypothetical model of intravalley e-h exchange of
the exciton that is part of the trion. In this model, the trion dispersion
becomes $\varepsilon(k)_{trion}=\hbar^{2}k^{2}/(3m_{e})+J^{intra}_{t}(k)$.
This affects the exciton-trion and trion-trion channel in the T-matrix, As
expected, we find that within this hypothetical model (not shown) that the
binding energy of the trion-trion complex is now reduced, and that of the
exciton-trion complex is slightly reduced.

The trion nonlinearity discussed so far simulates only direct interactions
between pump and probe induced trions. In principle, the pump-induced excitons
can lead to pump-induced free charge carriers if at least some of the exciton
dissociate. In that case, the pump-induced free charge carriers make forming
trions by the probe more probable. To model this scenario, we present in
Appendix B a phenomenological model for the 3rd-order response from pump-
induced free-carrier generation. Figure 13 shows that the main effect on the
DT is a reduction (or pulling down) of the DT feature at the trion, consistent
with the notion that the pump-induced carriers facilitate trion absorption of
the probe.

Figure 12:  (Color online.) Comparison of several T-matrix components for the
case of Fig. 11. Only results for full exchange strength are shown. In
addition to the exciton-exciton and exciton-trion component, the trion-trion
component is shown. A peak corresponding to a bound trion-trion complex can
clearly be seen at approximately -52 meV. Note that the trion-trion T-matrix
is small on this scale and has been multiplied by a factor of 5 for better
visibility.  Figure 13:  (Color online.) The effect of pump-induced free-
carrier generation on the cross-circular DT. The solid blue curve is the same
as in Fig. 10 at full exchange strength, which is without free-carrier
generation, and the dashed red curve is with free-carrier generation.

We now comment on the OAM characteristics of the DT. In the presence of
intervalley e-h exchange, spin-orbit coupling involved in the exciton
scattering creates OAM states that are not necessarily included in the
incident light field. In the results presented so far, we assume the incident
field has zero OAM and the transmitted field is filtered to also include only
zero OAM. In a typical experiment, the incident field has indeed only zero
OAM, but the transmitted field is not necessarily filtered. As discussed in
Sec. II.2, the T-matrix coupling is block diagonal in the total angular
momentum ($m_{tot}$) basis. The total angular momentum is related to the OAM
and the spins $\sigma_{1}$ and $\sigma_{2}$ of the 2 particles involved in the
scattering process via $m_{tot}=\mu+\sigma_{1}+\sigma_{2}$.

In the case of the incident field having only zero OAM, $\mu^{\prime}=0$,
there are only three total angular momentum states that contribute to the DT,
namely $m_{tot}=0,2,-2$. In the $m_{tot}=0$ subspace, the states written as
$\\{\sigma_{1},\sigma_{2},\mu\\}$ are $\\{+,+,-2\\}$, $\\{+,-,0\\}$,
$\\{-,+,0\\}$, $\\{-,-,2\\}$, Similarly, in the $m_{tot}=2$ subspace, the
states are $\\{+,+,0\\}$, $\\{+,-,2\\}$, $\\{-,+,2\\}$, $\\{-,-,4\\}$, and in
the $m_{tot}=-2$ subspace, the states are $\\{+,+,-4\\}$, $\\{+,-,-2\\}$,
$\\{-,+,-2\\}$, $\\{-,-,0\\}$. In the co-circular configuration, the created
exciton pair is in the $m_{tot}=2$ subspace. When intervalley exchange flips
one spin, the pair switches to the $\\{+,-,2\\}$ or $\\{-,+,2\\}$ state. Thus
the biexciton contributing to $T_{++++}$ in this case has OAM $\mu=2$. Within
each subspace, the T-matrix is a 4x4 matrix in the basis of the states just
given. Hence, each T-matrix component can be labeled by two OAM numbers, in
the notation of Sec. II.2 $T=T_{\mu,\mu^{\prime}}$ (suppressing all other
dependencies of the T-matrix). In the numerical evaluation, we restrict
ourselves to $\hat{T}^{\beta_{1}\beta_{2}}_{0,0}$ because we consider only
incident and detected beams in the direction normal to the monolayer plane.
The matrix elements $\hat{T}^{\beta_{1}\beta_{2}}_{\mu\mu^{\prime}}$ with non-
zero OAM $\mu$ and/or $\mu^{\prime}$ contribute to quantities of beams
prepared and/or detected in oblique directions. $T=T_{0,0}$ , $\mu$
$\mu^{\prime}=0$.

In Fig. 14, we show DT in the four (co-circular, cross-circular, co-linear,
cross-linear) configurations with the pump tuned at the exciton frequency. We
note that in Fig. 14 the usual relation between the DT in the co-circular and
cross-linear configuration, $DT_{co-circular}=2DT_{cross-linear}$, is not
fulfilled. This relation holds for the coherent 3rd-order response in
configurations where one resonance, for example the exciton, dominates the
linear response and can readily be verified from the relations given in Ref.
takayama-etal.04, . It can be shown analytically that the relation does not
hold in the case where more than one resonances, for example the exciton and
the trion, contribute to the response.

Figure 14:  (Color online.) Differential transmission (DT) in the co- and
cross-circular and co- and cross-linear configuration at full intravalley and
intervalley exchange.

##  IV Conclusion

In conclusion, we have shown how intravalley and intervalley e-h exchange
affect the nonlinear transmission spectra via the modification of the
dispersion relations of the particles involved in two-particle scattering, in
particular excitons. Using a model of a separable 2-particle potential, we
have shown analytically that, in 2 dimensions, a linear dispersion reduces the
2-particle bound state energy, and, depending on the interaction strength, can
eliminate the bound state. We have derived and numerically solved a
generalized Lippmann-Schwinger equation describing 2-particle scattering,
where the two particles can be two excitons, an exciton and a trion, or two
trions. The Lippmann-Schwinger equation accounts for intravalley e-h exchange
via the dispersion relations, and for intervalley e-h exchange resulting in
spin scattering and spin-orbit interaction. We found the T-matrix to be block
diagonal with each block labeled by the total angular momentum of the two
particles. We have analyzed the DT spectra and underlying T-matrix components
for the case of pump at exciton as a function of intravalley and intervalley
exchange coupling strength, and have seen that if both are at full strength,
the biexciton dip in DT (corresponding pump-induced absorption) merges with
the two-exciton continuum-scattering contribution, changing the DT lineshape
at the exciton from a positive-negative signal (indicating an exciton blue
shift) without exchange to one that is negative-positive-(weakly)negative
(similar to an exciton red shift) at full exchange strength.

Quite generally, the DT signal in $\chi^{(3)}$ is large when the light fields
are tuned to the linear (1st-order) resonances. That is why the DT signal is
largest when the probe is at the exciton. The specific lineshape of the DT
signal results from a complex interplay of the 1st-order susceptibilities with
the 2-particle T-matrix.

The trion also provides a 1st-order resonance, hence the DT signal becomes
large when the probe field is at the trion. We found that the lineshape at the
trion changes from mostly positive DT at zero e-h exchange in the exciton to a
redshift signal (negative-positive) at full exchange. The fact that the
lineshape at the trion is affected by the e-h exchange in the exciton is a
consequence of the bound trion-exciton resonance, which for our parameter
values shifts from below the linear trion resonance to above it as the e-h
exchange in the exciton is switched on. Note that the present calculations do
not include pump-induced doping, i.e. a mechanism whereby the pump field
provides additional free electrons that facilitate trion formation.

In the future, it would be desirable to extend the present analysis to the
case of incoherent nonlinear optical response, i.e. a regime where incoherent
scattering, in particular intervalley scattering and spin relaxation,
contributes to the optical response. It would also be desirable to extent the
present analysis beyond the separable interaction potential, similar to the
full numerical solution of the exciton-exciton T-matrix performed in Ref.
takayama-etal.02, .

###### Acknowledgements.

RB gratefully acknowledges financial support from the US National Science
Foundation (NSF) under grant number DMR 1839570, and the use of High
Performance Computing (HPC) resources supported by the University of Arizona.
JS acknowledges support from AFOSR Grant: FA9550-20-1-0217, NSF Grants:
DMR-2039329, ECCS-2054572, and ARO Grant: W911NF2010215 .

##  Appendix A The two-particle free, retarded Green’s function

The two-particle free, retarded Green’s function
$\hat{G}_{0}^{\beta_{1}\beta_{2}R}(\mathbf{q},\Omega)$ is defined by Eq. (21).
We derive here explicit expressions of this Green’s function through an
expansion in the eigenbasis of the single-particle Hamiltonian
$\hat{h}^{\beta}(\mathbf{p})$ Eq. (18). We denote the two eigenvalues of
$\hat{h}^{\beta}(\mathbf{p})$ by $\xi_{\beta,j}(p)$,
$j=u\,{\rm(upper)},\ell\,{\rm(lower)}$ and their corresponding eigenvectors by
$|\beta,u\rangle$ and $|\beta,\ell\rangle$. (We have omitted the momentum
labels $\mathbf{p}$ in the eigenvectors to reduce notational clutter.) They
are given by

$\displaystyle\xi_{\beta,u/\ell}(p)$
$\displaystyle=\epsilon_{\beta}(p)\pm|J^{inter}_{\beta}(p)|$ (51)
$\displaystyle\begin{pmatrix}\langle+|\beta,u\rangle\\\
\langle-|\beta,u\rangle\end{pmatrix}$
$\displaystyle=\frac{1}{\sqrt{2}}\begin{pmatrix}e^{-i\theta_{p}}\\\
e^{i\theta_{p}}\end{pmatrix}\quad,\quad\begin{pmatrix}\langle+|\beta,\ell\rangle\\\
\langle-|\beta,\ell\rangle\end{pmatrix}=\frac{1}{\sqrt{2}}\begin{pmatrix}e^{-i\theta_{p}}\\\
-e^{i\theta_{p}}\end{pmatrix}$ (52)

Expanded in the product eigenbasis of the two particles, a Green’s function
matrix element in the spin sub-basis becomes

$\langle\sigma_{1}\sigma_{2}|{G}_{0}^{\beta_{1}\beta_{2}R}(\mathbf{q},\Omega)|\sigma^{\prime}_{1}\sigma^{\prime}_{2}\rangle=\sum_{m,n=u,\ell}\frac{\langle\sigma_{1}\sigma_{2}|\beta_{1},m;\beta_{2},n\rangle\langle\beta_{1},m;\beta_{2},n|\sigma^{\prime}_{1}\sigma^{\prime}_{2}\rangle}{\hbar\Omega-\xi_{\beta_{1},m}({q}/2)-\xi_{\beta_{2},n}({q}/2)+i\gamma}$
(53)

Using the eigenvectors Eq. (52), we can write the Green’s function more
explicitly. When both particles are excitons, the Green’s function matrix in
the spin basis $(++,+-,-+,--)$ is as follows:

$\hat{G}_{0}^{xxR}(\mathbf{q},\Omega)=\frac{1}{4}\begin{pmatrix}s_{1}&s_{3}e^{-2i\theta_{q}}&s_{3}e^{-2i\theta_{q}}&s_{2}e^{-4i\theta_{q}}\\\
s_{3}e^{2i\theta_{q}}&s_{1}&s_{2}&s_{3}e^{-2i\theta_{q}}\\\
s_{3}e^{2i\theta_{q}}&s_{2}&s_{1}&s_{3}e^{-2i\theta_{q}}\\\
s_{2}e^{4i\theta_{q}}&s_{3}e^{2i\theta_{q}}&s_{3}e^{2i\theta_{q}}&s_{1}\end{pmatrix}$
(54)

where

$\displaystyle s_{1}$
$\displaystyle=g^{xx}_{uu}+g^{xx}_{\ell\ell}+2g^{xx}_{u\ell}\quad,\quad
s_{2}=g^{xx}_{uu}+g^{xx}_{\ell\ell}-2g^{xx}_{u\ell},$ $\displaystyle s_{3}$
$\displaystyle=g^{xx}_{uu}-g^{xx}_{\ell\ell}$ (55)

and

$\displaystyle g^{\beta_{1}\beta_{2}}_{mn}$
$\displaystyle=\frac{1}{\hbar\Omega-\xi_{\beta_{1},m}({q}/2)-\xi_{\beta_{2},n}({q}/2)+i\gamma}$
(56)

Because spin-flipping is absent in the trion, many elements in the exciton-
trion and trion-trion Green’s function matrices equal zero. The exciton-trion
Green’s function is

$\hat{G}_{0}^{xtR}(\mathbf{q},\Omega)=\frac{1}{4}\begin{pmatrix}r_{1}&0&r_{3}e^{-2i\theta_{q}}&0\\\
0&r_{1}&0&r_{3}e^{-2i\theta_{q}}\\\ r_{3}e^{2i\theta_{q}}&0&r_{1}&0\\\
0&r_{3}e^{2i\theta_{q}}&0&r_{1}\end{pmatrix}$ (57)

where

$\displaystyle r_{1}$
$\displaystyle=g^{xt}_{uu}+g^{xt}_{\ell\ell}+2g^{xt}_{u\ell}\quad,\quad
r_{3}=g^{xt}_{uu}-g^{xt}_{\ell\ell}$ (58)

and $g^{xt}_{\beta_{1}\beta_{2}}$ is defined in Eq. (55). The trion-trion
Green’s function matrix is the identity matrix multiplied by the factor
$g^{tt}_{uu}(=g^{tt}_{\ell\ell})$.

##  Appendix B Phenomenological model for pump-induced free-carrier generation

In this appendix we present a simple model that can account for pump-induced
free charge carriers, which in turn modify the trion susceptibility. This
phenomenological model does not contain microscopic details of how the pump
light creates free carriers. Some possibilities include two-photon absorption
of the pump, or thermal ionization of the excitons created by the pump, with
the latter being less likely if the experiment is performed at low
temperatures. Assuming that the free-carrier density is changed by the pump by
the amount $\Delta n_{e,i}^{\prime}$, i.e. $n_{e,i}^{\prime}\rightarrow
n_{e,i}^{\prime}+\Delta n_{e,i}^{\prime}$, the trion susceptibility, wich is
linear in the carrier density,will change by the amount

$\Delta\chi_{ij}^{t}(\omega_{pr})=-\frac{|f_{i}|^{2}\Delta
n_{e,ij}^{\prime}}{\hbar\omega_{pr}-\epsilon_{t}(0)+i\gamma_{t}}=\chi_{ij}^{t(1)}(\omega_{pr})\frac{\Delta
n_{e,ij}^{\prime}}{n_{e,i}^{\prime}}$ (59)

which is contributes to the third-order nonlinear response if we restrict
ourselves to the lowest order in the pump amplitude, where the change of the
density will be of second order. We account for possible resonance enhancement
of the pump-induced free carrier density with the following model for the
pump-induced density,

$\frac{\Delta
n_{e,ij}^{\prime}(\omega_{p},\omega_{p})}{n_{e,i}^{\prime}}=\sum_{kl}g_{ijkl}^{e}\sum_{\beta=x,t}\left(\frac{\chi_{k}^{\beta(1)}(\omega_{p})}{\zeta_{\beta,k}^{\ast}}\right)\left(\frac{\chi_{\ell}^{\beta(1)\ast}(\omega_{p})}{\zeta_{\beta,\ell}}\right)E_{p,k}E_{p,\ell}^{\ast}$
(60)

and the corresponding pump-induced free carrier contribution (superscript $f$)
to the third-order susceptibility,

$\chi_{ijkl}^{f(3)}(\omega_{pr},\omega_{p},\omega_{p})=\chi_{ij}^{t(1)}(\omega_{pr})g_{ijkl}^{e}\sum_{\beta=x,t}\left(\frac{\chi_{k}^{\beta(1)}(\omega_{p})}{\zeta_{\beta,k}^{\ast}}\right)\left(\frac{\chi_{\ell}^{\beta(1)\ast}(\omega_{p})}{\zeta_{\beta,\ell}}\right)$
(61)

which is added to the original $\chi^{(3)}$ terms in Eq. 46.

Since the probe does not mix the spins, and under the assumption that the pump
does not create intervalley spin coherence and that the spin $k$ in the valley
containing the pump-induced carriers must be different from the one associated
with the polarization $i$, we have
$g_{ijkl}^{e}=g_{0}^{e}\delta_{ij}\delta_{k\ell}[1-\delta_{ik}]$. In our
numerical evaluations, we take $g_{0}^{e}=6\times 10^{-13}$cm2.

For simple estimates, the following model for the exciton and trion oscillator
strengths are also useful. The effective exciton dipole matrix element is

$\mu_{i}=-q_{e}\frac{i\hbar\wp_{i}}{m_{0}\epsilon_{x}(0)}\phi_{x}(r=0)$ (62)

where $q_{e}$ is the absolute value of the electron charge in vacuum, $m_{0}$
the electron mass in vacuum, $\wp_{i}$ the interband momentum matrix element,
which in our case is independent of the spin subscript $i$ (see Ref. gu-
etal.13,  for a detailed analysis of the relation between momentum and dipole
matrix elements), and $\phi_{x}(r=0)$ the exciton wave function evaluated as
zero relative e-h coordinate $\mathbf{r}=\mathbf{r}_{e}-\mathbf{r}_{h}$.

According to Eq. (14) of Ref. glazov.2020, , the trion oscillator strength
differs from that of the exciton by a factor $4\pi
a^{2}_{tr}n_{e,i}^{\prime}$, where $a^{2}_{tr}$ is the trion Bohr radius.
Hence, we can write the trion susceptibility as

$\chi_{ij}^{t}(\omega_{pr})=-\frac{|\mu_{i}|^{2}4\pi
a^{2}_{tr}n_{e,ij}^{\prime}}{\hbar\omega_{pr}-\epsilon_{t}(0)+i\gamma_{t}}$
(63)

## References

  * (1) Jian Wang, Jeng-Yuan Yang, Irfan M. Fazal, Nisar Ahmed, Yan Yan, Hao Huang, Yongxiong Ren, Yang Yue, Samuel Dolinar, Moshe Tur, and Alan E. Willner, Nature Photonics 6, 488 (2012).
  * (2) J. R. Schaibley, H. Yu, G. Clark, P. Rivera, J. S. Ross, K. L. Seyler, W. Yao, and X. Xu, Nature Reviews 1, 1 (2016).
  * (3) W. Choi, N. Choudhary, G. H. Han, J. Park, D. Akinwande, and Y. H. Lee, Mater. Today 20, 116 (2017).
  * (4) S. Manzeli, D. Ovchinnikov, D. Pasquier, O. V. Yazyev, and A. Kis, Nat. Rev. 2, 17033 (2017).
  * (5) T. C. Berkelbach and D. R. Reichman, Ann. Rev. Condens. Matter Phys. 9, 379 (2018).
  * (6) Th. Mueller and E. Malic, NPJ 2D Materials and Applications 2, 29 (2018).
  * (7) Preparation and Crystal Growth of Materials with Layered Structures, edited by R.M.A. Lieth (Springer, New York, 1977).
  * (8) K. A. Duerloo, Y. Li, and E. J. Reed, Nat. Commun. 5, 4214 (2014).
  * (9) Y. Li, K. A. N. Duerloo, K. Wauson, and E. J. Reed, Nat. Commun. 7, 10671 (2016).
  * (10) Z. Hu, Z. Wu, C. Han, J. He, Z. Ni, and W. Chen, Chem. Soc. Rev. 47, 3100 (2018).
  * (11) D. Wickramaratne, F. Zahid, and R. K. Lake, J. Chem. Phys. 140, 124710 (2014).
  * (12) K. A. Duerloo, M. T. Ong, and E. J. Reed, J. Phys. Chem. Lett. 3, 2871 (2012).
  * (13) L. F. Mattheiss, Phys. Rev. B 8, 3719 (1973).
  * (14) K. F. Mak, C. Lee, J. Hone, J. Shan, and T. Heinz, Phys. Rev. Lett. 105, 136805 (2010).
  * (15) A. Ramasubramaniam, D. Naveh, and E. Towe, Phys. Rev. B 84, 205325 (2011).
  * (16) D. Xiao, G. Liu, W. Feng, X. Xu, and W. Yao, Phys. Rev. Lett. 108, 196802 (2012).
  * (17) T. Cheiwchanchamnangij and W. R. L. Lambrecht, Phys. Rev. B 85, 205302 (2012).
  * (18) G. Liu, W. Shan, Y. Yao, W. Yao, and D. Xiao, Phys. Rev. B 88, 085433 (2013).
  * (19) X. Li, F. Zhang, and Q. Niu, Phys. Rev. Lett. 110, 066803 (2013).
  * (20) K. Kosmider, J.W. Gonzales, and J. Fernandez-Rossier, Phys. Rev. B 88, 245436 (2013).
  * (21) C. Mai, A. Barrette, Y. Yu, Y. G. Semenov, K. W. Kim, L. Cao, and K. Gundogdu, Nano Lett. 14, 202 (2013).
  * (22) David MacNeill, Colin Heikes, Kin Fai Mak, Zachary Anderson, Andor Kormanyos, Viktor Zolyomi, Jiwoong Park, and Daniel C. Ralph, Phys. Rev. Lett. 114, 037401 (2015).
  * (23) H. Dery, Phys. Rev. B 94, 075421 (2016).
  * (24) Y. Li, J. Zhang, D. Huang, H. Sun, F. Fan, J. Feng, Z. Wang, and C. Z. Ning, Nature Nanotechnology 12, 987 (2017).
  * (25) F. Ceballos, M. G. Ju, S. D. Lane, X. C. Zeng, and H. Zhao, Nano Lett. 17, 1623 (2017).
  * (26) F. Ceballos, P. Zereshki, and H. Zhao, Phys. Rev. Materials 1, 044001 (2017).
  * (27) C. Ataca, M. Topsakal, E. Aktürk, and S. Ciraci, The Journal of Physical Chemistry C 115, 16354 (2011).
  * (28) A. Molina-Sanchez and L. Wirtz, Phys. Rev. B 84, 155413 (2011).
  * (29) A. Molina-Sanchez and L. Wirtz, Phys. Rev. B 84, 155413 (2011).
  * (30) K. Kaasbjerg, K. Thygesen, and K. Jacobsen, Phys. Rev. B 85, 165440 (2012).
  * (31) K. Kaasbjerg, K. Thygesen, and A. Jauho, Phys. Rev. B 87, 235312 (2013).
  * (32) G. Eda, H. Yamaguchi, D. Voiry, T. Fujita, M. Chen, and M. Chhowalla, Nano Lett. 11, 5111 (2011).
  * (33) A. M. Jones, H. Yu, N. J. Ghimire, S. Wu, F. Aivazian, J. S. Ross, B. Zhao, J. Yan, D. G. Mandrus, D. Xiao, W. Yao, and X. Xu, Nature Nanotechnology 8, 634 (2013).
  * (34) D. Y. Qiu, F. H. da Jornada, and S. G. Louie, Phys. Rev. Lett. 111, 216805 (2013).
  * (35) T. C. Berkelbach, M. S. Hybertsen, and D. R. Reichman, Phys. Rev. B 88, 045318 (2013).
  * (36) H. L. Liu, C. C. Shen, S. H. Su, C. L. Hsu, M. Y. Li, and L. J. Li, Appl. Phys. Lett. 105, 201905 (2014).
  * (37) X. Liu and T. Galfsky and Zh. Sun and F. Xia and E. Lin and Y. Lee and S. Kena-Cohen and V. Menon, Nature Photonics 9, 30 (2014).
  * (38) A. Chernikov, T. C. Berkelbach, H. M. Hill, A. Rigosi, Y. Li, O. B. Aslan, D. R. Reichman, M. S. Hybertsen, and T. F. Heinz, Phys. Rev. Lett. 113, 076802 (2014).
  * (39) Z. Ye, T. Cao, K O’Brien, H. Zhu, X. Yin, Y. Wang, S. Louie, and X. Zhang, Nature 513, 214 (2014).
  * (40) H. Yu, G. Liu, P. Gong, X. Xu, and W. Yao, Nature Communications 5, 3876 (2014).
  * (41) A. Steinhoff, M. Rosner, F. Jahnke, T. O. Wehling, and C. Gies, Nano Lett. 14, 3743 (2014).
  * (42) V. Tran, R. Soklaski, Y. Liang, and L. Yang, Phys. Rev. B 89, 235319 (2014).
  * (43) M. Z. Bellus, F. Ceballos, H. Y. Chiu, and H. Zhao, ACS Nano 9, 6459 (2015).
  * (44) F. Wu, F. Qu, and A. H. MacDonald, Phys. Rev. B 91, 075310 (2015).
  * (45) Y. You, X. Zhang, T. C. Berkelbach, M. S. Hybertsen, D. R. Reichman, and T. F. Heinz, Nature Physics 11, 477 (2015).
  * (46) T. Stroucken and S. W. Koch, Journal of Physics: Condensed Matter 27, 345003 (2015).
  * (47) G. Moody, C.K. Dass, K. Hao, C. Chen, L. Li, A. Singh, K. Tran, G. Clark, X. Xu, G. Berghäuser, E. Malic, A. Knorr, and X. Li, Nature Communications 6, (2015).
  * (48) I. Kylanpaa and H. Komsa, Physcial Review B 92, 205418 (2015).
  * (49) A. M. Jones, H. Yu, J. R. Schaibley, J. Yan, D. G. Mandrus, T. Taniguchi, K. Watanabe, H. Dery, W. Yao, and X. Xu, Nature Physics 12, 323 (2015).
  * (50) F. Hsu, W. Xie, Y. Lee, S. Lin, and C. Lai, Phys. Rev. B 91, 195312 (2014).
  * (51) H. Yu, X. Cui, X. Xu, and W. Yao, National Science Review 2, 57 (2015).
  * (52) P. Rivera, K. L. Seyler, H. Yu, J. R. Schaibley, J. Yan, D. G. Mandrus, W. Yao, and X. Xu, Science 351, 688 (2016).
  * (53) K. Hao, G. Moody, F. Wu, C. K. Dass, L Xu., C. H. Chen, L. Sun, M. Y. Li, L. J. Li, A. H. MacDonald, and X. Li, Nature Physics 12, 677 (2016).
  * (54) A. Singh, G. Moody, K. Tran, M. E. Scott, V. Overbeck, G. Berghäuser, J. Schaibley, E. J. Seifert, D. Pleskot, N. M. B Gabor, J. Yan, D. G. Mandrus, M. Richter, E. Malic, X. Xu, and X. Li, Phys. Rev. B 93, 041401 (2016).
  * (55) K. Hao, L. Xu, P. Nagler, A. Singh, K. Tran, C. Dass, C. Schüller, T. Korn, X. Li, and G. Moody, Nano Letters 16, 5109 (2016).
  * (56) K. Hao, L. Xu, J. F. Specht, P. Nagler, K. Tran, A. Singh, C. Dass, C. Schüller, T. Korn, M. Richter, A. Knorr, X. Li, and G. Moody, Nat. Comm. 8, 15553 (2017).
  * (57) C. M. Chow, H. Yu, A. M. Jones, J. R. Schaibley, M. Koehler, D. G. Mandrus, R. Merlin, W. Yao, and X. Xu, NPJ 2D Materials and Applications 1, 33 (2017).
  * (58) G. Scuri, Y. Zhou, A. A. High, D. S. Wild, C. Shu, K. De Greve, L. A. Jauregui, T. Taniguchi, K. Watanabe, P. Kim, M. D. Lukin, and H. Park, Phys. Rev. Lett. 120, 037402 (2018).
  * (59) L. Meckbach, T. Stroucken, and S. W. Koch, Phys. Rev. B 97, 035425 (2018).
  * (60) A. V. Stier, N. P. Wilson, K. A. Velizhanin, J. Kono, X. Xu, and S. A. Crooker, Phys. Rev. Lett. 120, 057405 (2018).
  * (61) Leone Di Mauro Villari, Ian Galbraith, and Fabio Biancalana, Phys. Rev. B 98, 205402 (2018).
  * (62) Dinh Van Tuan, Min Yang, and Hanan Dery, Phys. Rev. B 98, 125308 (2018).
  * (63) P.T. Mahon, R.A. Muniz, and J.E. Sipe, Phys. Rev. B 99, 235140 (2019).
  * (64) D. Jariwala, V. K. Sanwan, L. J. Lauhon, T. J. Marks, and M. C. Hersam, ACS Nano 8, 1102 (2014).
  * (65) H. Li, Y. Shi, M. H. Chiu, and L. J. Li, Nano Energy 18, 293 (2015).
  * (66) K. L. Seyler, J. R. Schaibley, P. Gong, P. Rivera, A. M. Jones, S. Wu, J. Yan, D. G. Mandrus, W. Yao, and X. Xu, Nature Nanotechnology 10, 407 (2015).
  * (67) Y. Ye, Z. J. Wong, X. Lu, X. Ni, H. Zhu, X. Chen, Y. Wang, and X. Zhang, Nature Photonics 9, 733 (2015).
  * (68) S. Wu, S. Buckley, J. R. Schaibley, L. Fengand, J. Yan, D. G. Mandrus, F. Hatami, W. Yao, J. Vucković, A. Majumbar, and X. Xu, Nature 520, 69 (2015).
  * (69) Alexander Steinhoff, Matthias Florian, Akshay Singh, Kha Tran, Mirco Kolarczik, Sophia Helmrich, Alexander W. Achtstein, Ulrike Woggon, Nina Owschimikow, Frank Jahnke, and Xiaoqin Li, Nature Phys. 14, 1199 (2018).
  * (70) D.Y. Qiu, T. Cao, and S.G. Louie, Phys. Rev. Lett. 115, 176801 (2015).
  * (71) T. Deilmann and K. S. Thygesen, 2D Materials 6, 035003 (2019).
  * (72) L. M. Schneider, S. S. Esdaille, D. A. Rhodes, K. Barmak, J. C. Hone, and A. Rahimi-Iman, Opt. Exp. 27, 37131 (2019).
  * (73) J. Hong, R. Senga, T. Pichler, and K. Suenaga, Phys. Rev. Lett. 124, 087401 (2020).
  * (74) Xiang Jiang, Qijing Zheng, Zhenggang Lan, Wissam A. Saidi, Xinguo Ren, and Jin Zhao, Science Advances 7, (2021).
  * (75) Jan R. Engelbrecht and Mohit Randeria, Phys. Rev. B 45, 12419 (1992).
  * (76) R. Takayama, N. H. Kwong, I. Rumyantsev, M. Kuwata-Gonokami, and R. Binder, Eur. Phys. J. B 25, 445 (2002).
  * (77) R. Takayama, N. H. Kwong, I. Rumyantsev, M. Kuwata-Gonokami, and R. Binder, J. Opt. Soc. Am. B 21, 2164 (2004).
  * (78) B. Simon, Ann. Phys. 97, 279 (1976).
  * (79) K. Chadan, N. N. Khuri, A. Martin, and T. T. Wu, J. Math. Phys. 44, 406 (2003).
  * (80) K. Yang and M. de Llano, Am. J. Phys. 57, 85 (1989).
  * (81) S. M. Apenko, J. Phys. A 31, 1553 (1998).
  * (82) J. R. Taylor, Scattering Theory: The Quantum Theory on Nonrelativistic Collisions (John Wiley and Sons, Inc., New York, 1972).
  * (83) Stefano Larentis, Hema C. P. Movva, Babak Fallahazad, Kyounghwan Kim, Armand Behroozi, Takashi Taniguchi, Kenji Watanabe, Sanjay K. Banerjee, and Emanuel Tutuc, Phys. Rev. B 97, 201407 (2018).
  * (84) C. Sieh, T. Meier, F. Jahnke, A. Knorr, S. W. Koch, P. Brick, M. Hübner, C. Ell, J. Prineas, G. Khitrova, and H.M. Gibbs, Phys. Rev. Lett. 82, 3112 (1999).
  * (85) M. Bayer, G. Ortner, O. Stern, A. Kuther, A. A. Gorbunov, A. Forchel, P. Hawrylak, S. Fafard, K. Hinzer, T. L. Reinecke, S. N. Walck, J. P. Reithmaier, F. Klopf, and F. Schäfer, Phys. Rev. B 65, 195315 (2002).
  * (86) B. Gu, N.H. Kwong, and R. Binder, Phys. Rev. B 87, 125301 (2013).
  * (87) M.M. Glazov, J. Chem. Phys. 153, 034703 (2020).