# Trion clustering structure and binding energy in 2D semiconductor materials:
Faddeev equations approach

K. Mohseni  Departamento de Física, Instituto Tecnológico de Aeronáutica,
DCTA,
12228-900, São José dos Campos, Brazil  M. R. Hadizadeh
<EMAIL_ADDRESS>College of Engineering, Science, Technology and
Agriculture, Central State University, Wilberforce, OH 45384, USA  Department
of Physics and Astronomy, Ohio University, Athens, OH 45701, USA  T. Frederico
Departamento de Física, Instituto Tecnológico de Aeronáutica, DCTA,
12228-900, São José dos Campos, Brazil  D. R. da Costa  Departamento de
Física, Universidade Federal do Ceará, Campus do Pici, 60455-900 Fortaleza,
Ceará, Brazil  A. J. Chaves  Departamento de Física, Instituto Tecnológico de
Aeronáutica, DCTA,
12228-900, São José dos Campos, Brazil

###### Abstract

In this work, we develop the basic formalism to study trions in semiconductor
layered materials using the Faddeev equations in momentum space for three
different particles lying in two dimensions. We solve the trion Faddeev
coupled integral equations for both short-range one-term separable Yamaguchi
potential and Rytova-Keldysh (RK) interaction applied to the MoS2 layer. We
devise two distinct regularization methods to overcome the challenge posed by
the repulsive electron-electron RK potential in the numerical solution of the
Faddeev equations in momentum space. The first method regulates the repulsive
interaction in the infrared region, while the second regulates it in the
ultraviolet region. By extrapolating the trion energy to the situation without
screening, the two methods gave consistent results for the MoS2 layer with a
trion binding energy of $-49.5(1)$ meV for the exciton energy of $-753.3$ meV.
We analyzed the trion structure for the RK and Yamaguchi potentials in detail,
showing their overall similarities and the dominant cluster structure, where
the strongly bound exciton is weakly bound to an electron. We found that this
property is manifested in the dominance of two of the Faddeev components over
the one where the hole is a spectator of the interacting electron pair.

##  I Introduction

Few-body problems appear in physics at different scales, ranging from
subatomic to celestial bodies. For semiconductors, the electron (e) and hole
(h) can form bound states due to the electrostatic attraction. Speculated
since the 30’s [1], the exciton appears in the absorption spectrum of
crystals, such as splitting of the lines in molecular crystals [2] or in band
edge absorption features [3]. Albeit the exciton is weakly bound due to the
intrinsic screening of traditional semiconductors, the electron-hole
interaction is fundamental to the understanding of the optical properties in
semiconductors and insulators [4], as it follows from the works of G.
Dresselhaus [5] and Elliot [6].

More complex few-body systems composed of holes and electrons were proposed by
Lampert in 1958 [7], such as the trion (eeh or ehh) and the biexciton.
However, the weak binding energy of the trion, which results from the strong
screening of the Coulomb interaction in ordinary materials, hindered its study
until the advent of quantum wells. Although observed in 1977 in the asymmetric
tail of exciton luminescence [8], a trion peak was only observed in 1993 [9]
in quantum wells, effectively a two-dimensional (2D) system, whose energies
were predicted to be an order of magnitude greater than in the three-
dimensional (3D) case [10] due to the quantum confinement effect.

One interesting aspect of trion physics already noted by Lampert is the
different limits as the hole and electron mass ratio changes, we can have the
analog of H${}_{2}^{+}$, H-, and e-e-e+. One should note that trions in
semiconductors differ from traditional three-body systems such as the triton
or the 4He3 atomic trimer, as the constituents of the trion have two
attractive interactions and one repulsive. We will explore this distinctive
feature in this paper.

As we already mentioned, it was only with the dimension reduction that trions
were first detected by observing the asymmetric tail of exciton luminescence
[8]. Three-particle bound states also appear in cold atom physics, where
through trapping, the dimension can be reduced continuously from 3D to 2D
[11]. The consequence is the disappearance of the Efimov effect and, together,
the log-periodicity of the wave function, which turns into a power law [12].

With the synthesis of 2D semiconductors [13], it was found that excitons can
have huge binding energies [14] as also trions [15], when compared to
traditional materials. This happens due to the reduced screening, as the
electric field lines lie outside the 2D semiconductor[16]. In those systems,
the strength of the interaction can be externally controlled by suitable
dielectric engineering [17]. Charge carriers in transition metal
dichalcogenides (TMDs) interact mainly via the screened Coulomb interaction
that in the classical regime is given by the Rytova-Keldysh potential,
obtained as the solution of the Poisson equation for an infinitesimal thin
dielectric slab [18].

Exciting prospects appear for few-body systems in novel 2D materials. There is
a plethora of different materials that hosts excitons, trions, and biexcitons,
such as TMDs, hexagon boron nitride, and graphene. In addition to that,
excitons and trions can strongly couple with light, forming exciton-polaritons
[19, 20] and trion-polaritons [21], respectively. As the Fermi energy
increases, e.g., with electrostatic doping, there is a transition of the trion
to an exciton-Fermi polaron [22]. The proximity effect [23], which originates
from short-range interactions, can also be used to tune the properties of
excitons and trions through the suitable choice of van der Waals
heterostructures, for example, the valley manipulation of excitons in TMDs due
to the coupling with CrI3 [24], whose magnetization can be controlled by an
external magnetic field, that breaks the time-reversal symmetry and the valley
degeneracy.

Several experiments have already probed trions in 2D materials. Observations
of large trion binding energies in MoS2 reported experimental values between
20–43 meV for samples deposited on SiO2 substrates [25, 26, 27, 28] and 80 meV
on suspended samples [29], while Ref. [30] measured for different substrates
and found an extrapolation curve for the suspended case of 44 meV. Besides the
dependence on the dielectric environment [30], the trion binding energy
depends on the doping [15], and also on the temperature [31]. In this work
(Sec. V), we will discuss the trion in an undoped suspended MoS2 layer at zero
temperature, thus we do not expect an exact agreement with experimental
measurements that are performed in a finite temperature and with residual
doping.

There are already several theoretical calculations on trion binding energies
[32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44]. The authors in Ref. [37]
found a good agreement between the multiband and effective mass models, thus
justifying our choice of using the effective mass approach in this work. In
Ref. [43], it is reported a calculated binding energy for the trion in MoS2 of
33.6 meV with a variationally optimized orbital approach ($m_{e}/m_{0}=$0.47,
$m_{h}/m_{0}=$0.54 and $r=\,44.68$ Å), in Ref. [34] the value of 33.7 meV was
obtained by the stochastic variational method, and Ref. [44] reported the
value of 32.1 meV by using an imaginary time evolution method for numerically
solving the trion Schrödinger-like Hamiltonian. Based on the ab-initio many-
body theory, a converged negatively charged intralayer trion binding energy
was found to be of 58 meV [38] with an exciton binding energy of $-0.76$ eV.
The use of Faddeev equations in configuration space to calculate the charge
positive and negative trion energies in various TMDs was reported in Refs.
[39, 40, 41, 42].

Our goal in this work is to study negatively charged trions within the Faddeev
equations approach in momentum space and explore both the binding and
structural properties of the trion in a MoS2 layer. Within the adopted method,
each Faddeev component is computed, which sums up the total wave function and
carries information about each pair that composes the three-body state. The
numerical convergence due to the repulsion between the electrons is a
challenge, and to overcome this, we use two different approaches to regularize
the electron-electron interaction at both long and short distances, to weaken
the repulsion and turn the numerical calculations more accurate and finally,
we extrapolate both results to compute the trion energy accurately.
Furthermore, we cross-check the accuracy of our calculated trion energy by
computing the expectation value of the Hamiltonian with the wave function. Our
work addresses the following main points: (i) we provide a general discussion
of the wave function properties for a 2D trion obtained within the Faddeev
equations approach; (ii) a theoretical-numerical calculation of the trion
binding energy in freestanding monolayer MoS2 with different regularization
schemes, with the accuracy checked by computing the expectation value of the
Hamiltonian; and (iii) the degree of clusterization of the trion weakly bound
state.

The assumed theoretical framework is presented in Sec. II, where we derive the
Faddeev equations in 2D considering three different particles. In Sec. III, we
present results for the Yamaguchi model, a non-local separable and short-range
potential [45, 46], considering three attractive interactions and also for two
attractive and one repulsive potential. In Sec. IV, we illustrate the cluster
structure of the Yamaguchi model for trions. Sec. V is devoted to presenting
the results for the Rytova-Keldysh potential, where two different
regularization procedures are introduced to compute the trion binding energy.
In Sec. VI, it is illustrated the cluster structure of the trion by showing
results for the total wave function and its Faddeev components, which we
compare with the structure of the wave functions obtained by the two potential
models. In Sec. VII, we summarize the main findings of our study. This work is
accompanied by six appendices where we detail our framework and numerical
methods.

##  II  Faddeev equations for 3B bound states in two-dimension

We consider the effective mass Hamiltonian for three-different particles

$H=\sum_{i=1}^{3}\left(\frac{k_{i}^{2}}{2m_{i}}+V_{i}\right),$ (1)

with $V_{i}\equiv V_{i}(r_{j}-r_{k})$, $i\neq j\neq k$ and $m_{i}$ being the
mass of the $i$-th particle. In the case of trions, this corresponds to the
Wannier-Mott model. The Schrödinger equation for the bound state of three
different particles interacting with pairwise interactions $V_{i}\equiv
V_{jk}$ is given by

$\displaystyle\Psi=\sum_{i=1}^{3}G_{0}V_{i}\Psi=\sum_{i=1}^{3}\psi_{i},$ (2)

where $\psi_{i}=G_{0}V_{i}\Psi$ are the Faddeev components,
$G_{0}=(E-H_{0})^{-1}$ is the free propagator with three-body (3B) binding
energy $E$ and free Hamiltonian $H_{0}$. Three Faddeev components $\psi_{i}$
satisfy the following coupled equations

$\psi_{i}=G_{0}\,t_{i}\,(\psi_{j}+\psi_{k}\,),$ (3)

where $\\{i,j,k\\}$ is a cyclic permutation of $\\{1,2,3\\}$. The two-body
(2B) transition operators $t_{i}$ are defined by the Lippmann-Schwinger
equation

$t_{i}=V_{i}+V_{i}G_{0}t_{i}.$ (4)

Figure 1: Jacobi momenta convention used through the paper. $\\{i,j,k\\}$
correspond to the indices associated with the three particles, and $p_{i}$ and
$q_{i}$ are their Jacobi momenta given in Eq. (5).

To solve the coupled Faddeev equations (3) in momentum space, we consider the
3B basis states $|{\bf p}_{i}{\bf q}_{i}\rangle$, composed of two Jacobi
momenta (see Fig. 1), which are defined in terms of the single particle
momenta as

$\displaystyle{\bf p}_{i}$ $\displaystyle=$ $\displaystyle\frac{m_{k}{\bf
k}_{j}-m_{j}{\bf k}_{k}}{m_{j}+m_{k}},$ (5) $\displaystyle{\bf q}_{i}$
$\displaystyle=$ $\displaystyle\frac{m_{j}+m_{k}}{m_{i}+m_{j}+m_{k}}\left({\bf
k}_{i}-\frac{m_{i}}{m_{j}+m_{k}}({\bf k}_{j}+{\bf k}_{k})\right),$ (6)

where ${\bf p}_{i}$ is the relative momentum of the pair $jk$, and ${\bf
q}_{i}$ is the relative momentum of the third particle $i$ with respect to the
pair $jk$. The completeness relation of 3B basis states in 2D is defined as

$\displaystyle\int d^{2}p_{i}\,\int d^{2}q_{i}\,|{\bf p}_{i}{\bf
q}_{i}\rangle\langle{\bf p}_{i}{\bf q}_{i}|=1.$ (7)

The projection of coupled Faddeev equations (3) on 3B basis states $|{\bf
p}_{i}{\bf q}_{i}\rangle$ leads to three coupled 2D integral equations

$\displaystyle\psi_{i}(p_{i},q_{i},\phi_{i})=\dfrac{1}{E_{3B}-\dfrac{p_{i}^{2}}{2\mu_{jk}}-\dfrac{q_{i}^{2}}{2\mu_{i,jk}}}\int_{0}^{\infty}dp^{\prime}_{i}\,p^{\prime}_{i}\int_{0}^{2\pi}d\phi^{\prime}_{i}$
(8) $\displaystyle\times\
t_{i}(p_{i},p^{\prime}_{i},\phi^{\prime}_{i};\epsilon_{i})\Big{[}\psi_{j}(\mathcal{P}_{ji},\mathcal{Q}_{ji},\phi_{ji})+\psi_{k}(\mathcal{P}_{ki},\mathcal{Q}_{ki},\phi_{ki})\Big{]},\quad$
(9)

where $\mu_{jk}=\dfrac{m_{j}m_{k}}{m_{j}+m_{k}}$ and
$\mu_{i,jk}=\dfrac{m_{i}(m_{j}+m_{k})}{m_{i}+m_{j}+m_{k}}$ are 2B and 3B
reduced masses. The shifted momenta and angle quantities ${\mathcal{P}}_{ij}$,
$\mathcal{Q}_{ij}$, and $\phi_{ij}$ are defined in Eq. (45). The details of
the derivation are given in Appendix A.

The non-partial-wave 2B $t-$matrices
$t_{i}(p_{i},p^{\prime}_{i},\phi^{\prime}_{i};\epsilon_{i})$, with 2B
subsystem energies $\epsilon_{i}=E_{3B}-\dfrac{q^{2}_{i}}{2\mu_{i,jk}}$, can
be obtained from the summation of partial wave (PW) $t-$matrices
$t_{m}(p_{i},p^{\prime}_{i};\epsilon_{i})$ as

$\displaystyle t_{i}(p_{i},p^{\prime}_{i},\phi^{\prime}_{i};\epsilon_{i})$
$\displaystyle=$
$\displaystyle\frac{1}{2\pi}\,\sum^{\infty}_{m=0}\varepsilon_{m}\cos(m\phi^{\prime}_{i})\,t_{m}(p_{i},p^{\prime}_{i};\epsilon_{i}),\quad$
$\displaystyle\mbox{with}\quad\varepsilon_{m}=\begin{cases}1&m=0\\\ 2&m\neq
0\end{cases},$

where PW projected 2B $t-$matrices in channel $m$, i.e.,
$t_{m}(p_{i},p^{\prime}_{i},\epsilon_{i})$, can be obtained from the solution
of inhomogeneous Lippmann-Schwinger integral equation as

$\displaystyle t_{m}(p_{i},p^{\prime}_{i};\epsilon_{i})$ $\displaystyle=$
$\displaystyle
V_{m}(p_{i},p^{\prime}_{i})+\int_{0}^{\infty}dp^{\prime\prime}_{i}\,p^{\prime\prime}_{i}\,V_{m}(p_{i},p^{\prime\prime}_{i})\,$
(11) $\displaystyle\times$
$\displaystyle\frac{1}{\epsilon_{i}-\dfrac{p^{\prime\prime
2}_{i}}{2\mu_{jk}}}\,t_{m}(p^{\prime\prime}_{i},p^{\prime}_{i};\epsilon_{i}),$

with PW projected interactions obtained from

$\displaystyle
V_{m}(p_{i},p^{\prime}_{i})=\int_{0}^{2\pi}\,d\phi^{\prime}_{i}\,V(p_{i},p^{\prime}_{i},\phi^{\prime}_{i})\,\cos(m\phi^{\prime}_{i}).$
(12)

In the following sections, we present our numerical results for the solution
of the coupled Faddeev equations (8) for trions with two different potentials:
the short-range Yamaguchi potential and the long-range Rytova-Keldysh
potential. The numerical solution details are provided in Appendix E.

##  III Trions: Yamaguchi potential

To test the formulation of the coupled Faddeev integral equations (8) and to
validate our numerical solution, we first utilize the one-term separable
potential with Yamaguchi-type form factors [45, 46]

$V(p,p^{\prime})=-\lambda g(p)g(p^{\prime}),\quad
g(p)=\frac{1}{(\beta^{2}+p^{2})^{m}},$ (13)

where the potential strength $\lambda$ can be obtained from the pole property
of the 2B $t-$matrix at the 2B binding energy. We present our numerical
results for 3B binding energies and wave functions using two different
interaction combinations: (i) three attractive, and (ii) two attractive and
one repulsive Yamaguchi-type potential, considering three particles with
identical masses.

We use three attractive Yamaguchi interactions to evaluate our formalism and
computer codes for solving the general form of three coupled Faddeev integral
equations in 2D. We should mention that this case is paradigmatic in cold-atom
physics [47], and it also appears in the formulation of the three-magnon bound
state problem [48]. In our context, these calculations are valuable as a
preparation for the practical application involving trions.

In order to clarify the assumed notation here, we denote the 3B binding energy
as $E_{3B}$, being defined as the eigenvalue of the 3B Hamiltonian in Eq. (1),
whereas the trion binding energy $E_{\mathrm{t}}$ is defined as the splitting
between the 2B and 3B binding energies

$E_{\mathrm{t}}=E_{3B}-E_{2B}.$ (14)

In Table 1, we present 3B binding energies obtained from the solution of the
three coupled Faddeev integral equations (8) for three identical particles
(mass = 1) interacting with three attractive Yamaguchi interactions. The input
2B $t-$matrices are obtained from the $s-$wave interactions. The calculated 3B
and 2B binding energy ratios with different potential strengths $\lambda$ and
form factor powers $m$ are in excellent agreement with the corresponding
results from Ref. [49].

By solving the coupled Faddeev integral equations and having 3B binding energy
and the Faddeev components, one can calculate the 3B wave function as a
summation of three Faddeev components. In Appendix B, we show the details of
the derivation of the 3B wave function in momentum space.

Table 1: 2B and 3B binding energies $E_{2B}$ and $E_{3B}$ calculated for three attractive Yamaguchi-type potentials with form factor parameter $\beta=1$ and different powers $m$. The potential strength $\lambda$ is fitted to reproduce the desired 2B binding energy $E_{2B}$. The ratio of 3B and 2B binding energies $E_{3B}/E_{2B}$ are compared with corresponding results from Ref. [49]. The calculations are done with $\hbar c=\text{ mass}=1$.  $\lambda$ |  | $E_{2B}$ |  | $E_{3B}/E_{2B}$ |  |  $E_{3B}/E_{2B}$ [49]
---|---|---|---|---|---|---
$m=1$
$0.0602$ |  | $-0.0019$ |  | $9.21$ |  | $9.21$
$0.0863$ |  | $-0.0100$ |  | $6.83$ |  | $6.83$
$0.1838$ |  | $-0.1000$ |  | $4.58$ |  | $4.58$
$m=2$
$0.0801$ |  | $-0.0032$ |  | $7.30$ |  | $7.30$
$0.1400$ |  | $-0.0211$ |  | $5.14$ |  | $5.14$
$m=4$
$0.0481$ |  | $-0.0001$ |  | $11.54$ |  | $11.53$
$0.0731$ |  | $-0.0010$ |  | $7.91$ |  | $7.91$
$0.1861$ |  | $-0.0200$ |  | $4.55$ |  | $4.55$
$m=10$
$0.0561$ |  | $-0.0001$ |  | $10.05$ |  | $10.05$
$0.0923$ |  | $-0.0010$ |  | $6.61$ |  | $6.61$
$0.1562$ |  | $-0.0050$ |  | $4.89$ |  | $4.99$
Table 2: Expectation values (EV) of 3B free Hamiltonian $\langle H_{0}\rangle$, pair interactions $\langle V_{i}\rangle$, total 2B interactions $\langle V\rangle$, 3B Hamiltonian $\langle H\rangle$, and eigenvalue ${E}_{3B}$ calculated for Yamaguchi-type potentials (three attractive (3A) in the second column, two attractive plus one repulsive interaction (2A+R) in the third column), given in Eq. (13) with form factor parameters $\beta=m=1$, and the potential strength $\lambda$ that reproduces dimer binding energy $E_{2B}=-0.1$. The relative percentage difference is Error=$|(\langle H\rangle-{E}_{3B})/{E}_{3B}|\times 100\%$. The calculations are done with $\hbar c=\text{mass}=1$. EV | 3A | 2A + R
---|---|---
$\quad\langle H_{0}\rangle\quad$ | $+0.46887$ | $+0.15756$
$\langle V_{1}\rangle$ | $-0.30904$ | $+0.03260$
$\langle V_{2}\rangle$ | $-0.30904$ | $-0.15526$
$\langle V_{3}\rangle$ | $-0.30904$ | $-0.15523$
$\langle V\rangle$ | $-0.92712$ | $-0.27789$
$\langle H\rangle$ | $-0.45825$ | $-0.12033$
${E}_{3B}$ | $-0.45824$ | $-0.12034$
Error | $+0.00218$ | $+0.00831$

To test the accuracy of the 3B wave function in momentum space, in Table 2, we
compare the expectation values of 3B Hamiltonian with the calculated 3B
binding energy for the factor parameters $\beta=m=1$ and the potential
strength $\lambda$ that reproduces 2B binding energy $E_{2B}=-0.1$. The
separable potential strength is obtained by introducing Eq. (13) in Eq. (11),
and considering that the t-matrix has a pole at $\epsilon_{i}=E_{2B}$, then
the potential strength can be obtained by

$\lambda^{-1}=-2\pi\int_{0}^{\infty}dp_{i}^{\prime\prime}p_{i}^{\prime\prime}\frac{|g(p^{\prime\prime})|^{2}}{E_{\mathrm{2B}}-\frac{{p_{i}^{\prime\prime}}^{2}}{2\mu_{jk}}}.$
(15)

The expectation values of the kinetic energy and potential in the exciton
state are, in this case, given by

$\langle H_{0}\rangle=0.138407\quad\text{and}\quad\langle V\rangle=-0.238407.$
(16)

As we can see in Table 2, the 3B binding energy and the expectation value of
Hamiltonian are in excellent agreement. The details of the calculation of
expectation values of Hamiltonian $\langle H\rangle$ from the expectation
value of 3B free Hamiltonian $\langle H_{0}\rangle$ and 2B interactions
$\langle V_{i}\rangle$ are given in Appendix C.

Some interesting qualitative aspects can be seen in Table 2. When the sign of
the potential $V_{1}$ is changed and becomes repulsive, the state swells due
to the dramatic decrease in the splitting of the 2B and 3B energies, namely
from $|E_{3B}-E_{2B}|=0.3582$ to $0.0203$. Consequently, the kinetic energy is
also reduced to about one-third of the value obtained with only attractive
potentials. Due to the repulsion, the wave function is depleted when the
relative distance between particles 2 and 3 lies in the range of the
potential, and the expectation value $\langle V_{1}\rangle$ turns to be
negative and reduced to one-tenth with respect to the attractive case.
Furthermore, the expectation values of $\langle V_{2}\rangle$ and $\langle
V_{3}\rangle$ are also halved, and due to our choice of mesh points, the
equality $\langle V_{2}\rangle=\langle V_{3}\rangle$ is fulfilled to 0.02%,
which is reflected in error around 0.008% in the computation of $\langle
H\rangle$, which is four times larger than the error in the attractive case.

---
Figure 2: Left panel: 3D plot of the 3B wave function for three attractive Yamaguchi-type potentials (3A). Right panel: relative error. All plots are for the angle $\phi_{1}=0$. The results are obtained with form factor parameters $\beta=m=1$, and the potential strength $\lambda$ that reproduces 2B binding energy $E_{2B}=-0.1$. The calculations are done with $\hbar c=\text{mass}=1$. |
---|---
Figure 3: Left panel: 3D plot of the 3B wave function for one repulsive
$(V_{1})$ and two attractive Yamaguchi-type potentials (2A + R). Right panel:
relative error. All plots are for the angle $\phi_{1}=0$. The results are
obtained with form factor parameters $\beta=m=1$, and the potential strength
$\lambda$ that reproduces 2B binding energy $E_{2B}=-0.1$. The calculations
are done with $\hbar c=\text{mass}=1$.

The cluster structure of the trion is indicated by its kinetic and potential
energies when comparing the results of $\langle H_{0}\rangle$, $\langle
V_{2}\rangle$ and $\langle V_{3}\rangle$ from Table 2, with the expectation
values of the kinetic and potential energies of the exciton given in Eq. (16).
The screening of the hole that composes the exciton in the trion [41] weakens
the electron-exciton interaction favoring the formation of the remarkable
cluster structure. Note that the electrons should be in a singlet spin state
or an antisymmetric combination of different valley states.

##  IV Yamaguchi Trion Clusterization

In Figs. 2 and 3, we show the magnitude of the 3B wave function obtained in
two cases with Yamaguchi interactions as a function of the magnitude of Jacobi
momenta $p_{1}$ and $q_{1}$ for the angle between them $\phi_{1}=0$. The first
case shown in Fig. 2 is the eigenstate of the 3B system with three attractive
potentials (3A), and the second case shown in Fig. 3 is the wave function for
the weakly bound state obtained with one repulsive potential, $V_{1}$, and two
attractive ones (2A+R). We also show the relative percentage error for the
verification of the Schrödinger equation, defined in Appendix D, with 3B wave
function and binding energy.

The plot of the wave function for the 3A case for ${\bf p}_{1}$ and ${\bf
q}_{1}$ aligned presented in the left panel of Fig. 2 shows that the momentum
distribution is somewhat symmetric due to the identical masses and the bosonic
symmetry of the system. When the repulsive potential $V_{1}$ is introduced in
the case 2A+R, seen in the left panel of Fig. 3, the wave function develops a
node line as well becomes more sharply peaked around the origin, the latter
due to the small binding energy (see Table 2). The numerical accuracy of our
calculations is checked through the ratio $|(E\Psi-H\Psi)/(E\Psi)|$ and
presented in the right panels of Figs. 2 and 3 for the 3A and 2A+R cases,
respectively. As expected, the results for the 3A case show quite good
numerical accuracy, while the 2A+R results, mainly outside the node, are also
accurate. As expected, the region of the largest errors for the 2A+R case in
the right panel of Fig. 3 follows the node of the 3B wave function.

|
---|---
Figure 4: Contour plot of the 3B wave function calculated for three Yamaguchi-
type potentials as a function of the magnitude of Jacobi momenta. $V_{1}$
repulsive, $V_{2}$ and $V_{3}$ attractive. The plots are for the angle
$\phi_{1}=\phi_{2}=0$. The results are obtained with form factor parameters
$\beta=m=1$ and the potential strength $\lambda$ that reproduces 2B binding
energy $E_{2B}=-0.1$. The calculations are done with $\hbar c=\text{mass}=1$.

In Fig. 4, we present the contour plots of the 2A+R model in the ($p_{1}\times
q_{1}$)–plane for $\phi_{1}=0$ (left panel) and ($p_{2}\times q_{2}$)–plane
with $\phi_{2}=0$ (right panel). The node line, visible in the left panel of
Fig. 4, similarly to the left panel in Fig. 3, comes from the cancellation
between $\psi_{1}$ and $\psi_{2}+\psi_{3}$ due to the reversed sign of
$\psi_{1}$ with respect to $\psi_{2}$ and $\psi_{3}$ from the repulsive
potential $V_{1}$, then

$\psi_{1}(p_{1},q_{1},0)=-\psi_{2}(p_{2},q_{2},\phi_{2})-\psi_{3}(p_{3},q_{3},\phi_{3})\,,$
(17)

where $\phi_{2,3}=0$, and these relations implicitly define the node line,
understood by rewriting the momenta labeled by 2 and 3, according to

$\displaystyle{\bf p}_{2}=-\frac{1}{2}{\bf p}_{1}-\frac{3}{4}{\bf
q}_{1},\quad{\bf q}_{2}={\bf p}_{1}-\frac{1}{2}{\bf q}_{1},$ (18a)
$\displaystyle{\bf p}_{3}=-\frac{1}{2}{\bf p}_{1}+\frac{3}{4}{\bf
q}_{1},\quad{\bf q}_{3}=-{\bf p}_{1}-\frac{1}{2}{\bf q}_{1}.$ (18b)

The node line is barely seen in the right panel of Fig. 4, with the momenta
expressed in terms of $p_{2}$ and $q_{2}$ with $\phi_{2}=0$.

Another property of the wave function is the well-defined maximum seen in both
Figs. 3 and 4, which can be qualitatively understood by a semi-classical
picture and the prevalence of the cluster structure. This dominant
configuration suggests that the electron and hole (denoted as particle 1) are
very close and “moving together”, which relates the momenta $p_{1}$ and
$q_{1}$ and provides an interpretation of the pattern of the maximum found in
the left panel of Fig. 3 and around the whitish-yellow color of Fig. 4,
namely, along the line $p_{1}\propto q_{1}$ for the 2A+R model. In other
words, the relative velocity between the two electrons is the same as the one
between the far-apart electron and the hole, which forms the strongly bound
exciton, which is very clear in the situation where $\phi_{1}=0$. This also
explains the obtained pattern of the maximum of the wave function, with two
branches observed in the right panel of Fig. 3 in the ($p_{2}\times
q_{2}$)–plane, as we shall discuss in more detail in what follows. We would
like to draw the reader’s attention to the practical significance of such
plots, which provide insights into the regions where the wave function is more
substantial. This information is crucial for distributing mesh points
appropriately to obtain precise solutions to the Faddeev equations.

Table 3: The inner product of the Faddeev components $\langle\psi_{i}|\psi_{j}\rangle$ and their contribution to the normalization of the 3B wave function $|\Psi\rangle$. | $j=1$ | $j=2$ | $j=3$
---|---|---|---
$i=1$ | $0.0560$ | $-0.1047$ | $-0.1047$
$i=2$ | $-0.1047$ | $0.4400$ | $0.2413$
$i=3$ | $-0.1047$ | $0.2413$ | $0.4403$

We would like to emphasize that each Faddeev component of the wave function in
our system carries the asymptotic form of the total wave function in each
pairwise interaction channel [50]. Specifically, in our context, the Faddeev
component $\psi_{1}$ at asymptotically large distances of the hole to the
center of mass of the electron-electron interacting pair decays exponentially,
indicating that the two electrons are in a continuum state. Similarly, at
asymptotically large distances of the spectator electron (particle 2) to the
center of mass of the electron-hole system (particles 3 and 1), the Faddeev
component $\psi_{2}$ decays exponentially, signifying that this pair
necessarily forms the strongly bound exciton state. The same reasoning applies
to $\psi_{3}$, where the electron-hole pair is formed by particles 1 and 2.

A scheme illustrates the clustering of the wave function: $\psi_{2}\sim$
[3(e)1(h)]—2(e) and $\psi_{3}\sim$ [1(h)2(e)]—3(e), which should be the two
dominant configurations, with the electrons in the spin singlet state or an
antisymmetric combination of different valley states. Indeed in Table 3, one
observes that the Faddeev component $\psi_{1}$ is suppressed with respect to
$\psi_{2}$ and $\psi_{3}$ by one order of magnitude considering the inner
products. With that in mind, we should now look to Fig. 4 (right panel) for
the modulus of the total wave function, $|\Psi|$ in the ($p_{2}\times
q_{2}$)–plane. We identify two branches where $|\Psi|$ is larger: one for
small $q_{2}$ and a diagonal one. The lower branch corresponds to the
contribution of $\psi_{2}$ for $q_{2}\approx 0$, which is the relative
momentum of the weakly bound spectator particle 2(e) with respect to the
strongly correlated pair of particles 1 and 3. The spread in the values of
$p_{2}$ is associated with the small size of the strongly bound exciton in the
[3(e)1(h)]—2(e) configuration. The diagonal branch, $p_{2}\propto q_{2}$,
where the momentum probability density is enhanced, corresponds to the
dominance of $\psi_{3}$ associated with the [1(h)2(e)]—3(e) configuration. In
this case, electron 2 moves together with hole 1, as the exciton is strongly
bound, and electron 3 is the spectator.

##  V Trions: Rytova-Keldysh potential

Building on our understanding gained from the 2A+R Yamaguchi potential model
in 2D, we now study the trion binding energy and structure for the MoS2 layer
with the Rytova-Keldysh potential. The Rytova-Keldysh electron-hole (e-h)
[electron-electron (e-e)] interaction in momentum space is given by [51]

$V_{\genfrac{}{}{0.0pt}{}{eh}{ee}}(q)=\pm\frac{1}{4\pi^{2}}\left(\frac{1}{4\pi\epsilon_{0}}\frac{2\pi
e^{2}}{q(1+r_{0}q)}\right),$ (19)

where the momentum transfer is defined by $|{\bf q}|=|{\bf p}-{\bf
p}^{\prime}|$. The parameters of the e-e and e-h potentials for the MoS2 layer
are given in Table 4. The value of the screening length $r_{0}$ is fitted to
give an exciton binding energy of $-753$ meV in agreement with the value
obtained from the measurement of the exciton position in the absorption
spectrum of a suspended MoS2 layer [52] and the corresponding GW bandgap [53].
For our reference, the expectation values of the kinetic and potential
energies in the exciton state are

$\langle H_{0}\rangle=214.64\,\text{meV}\quad\text{and}\quad\langle
V\rangle=-967.96\,\text{meV}\,,$ (20)

which are related to the manifestation of this strongly bound two-particle
system.

Table 4: The used parameters in our calculations for the Rytova-Keldysh electron-hole and electron-electron interactions, defined in Eq. (19), for monolayer MoS2.  $r_{0}$ |  $27.05$ Å
---|---
$\epsilon_{0}/e^{2}$ |  $\dfrac{1}{4\pi\alpha}\cdot\dfrac{1}{\hbar c}$ K${}^{-1}\cdot$ Å-1
$\alpha$ | $137.035999084$
$m_{e}$ |  $0.47\ m_{0}$ [54]
$m_{h}$ |  $0.54\ m_{0}$ [54]
$m_{0}$ |  $0.510998950$ MeV
1 eV |  $1.160451812\cdot 10^{4}$ K
$\hbar c$ |  $1973.269804$ eV$\cdot$ Å

The example studied in Sec. IV has already taught us that the accuracy of our
numerical solution of the Faddeev equations decreases in the case of 2A+R
Yamaguchi potential with respect to the 3A attractive case (cf. Table 2). This
expected behavior of our numerical solutions is due to the small trion binding
energy and the node in the wave function. On top of that, considering that the
Rytova-Keldysh potential is of a longer range when contrasted to the Yamaguchi
model, the numerical solution becomes more challenging due to the competition
between attraction and repulsion with the same strength. To make this issue
numerically amenable, the repulsive Rytova-Keldysh potential between the
electrons is screened by two different regulators [55], namely

$V(q)\to(1-e^{-l_{0}q})V_{ee}(q)\quad\text{ or}\quad
e^{-l_{0}\,q}V_{ee}(q)\,,$ (21)

where in the first case, the Rytova-Keldysh potential is damped at small
momentum or large distance, while the second one is at large momentum or small
distance.

Figure 5: The screening function $V(q)/V_{ee}(q)$ as a function of the momentum transfer $q$ for (blue curve) $V(q)\to e^{-l_{0}\,q}V_{ee}(q)$ and (red curve) $V(q)\to(1-e^{-l_{0}q})V_{ee}(q)$ with screening parameter $l_{0}=10$ Å. Table 5: Expectation values (EV) in meV of the 3B free Hamiltonian $\langle H_{0}\rangle$, pair interactions $\langle V_{i}\rangle$, total 2B interactions $\langle V\rangle$, 3B Hamiltonian $\langle H\rangle$, and binding energy $E_{3B}$ in meV calculated for Rytova-Keldysh potentials (two attractive plus one repulsive interaction (2A+R)) given in Eq. (19) with screening parameter $l_{0}=100$ Å in the scheme $V_{1}(q)\to(1-e^{-l_{0}q})V_{ee}(q)$. The relative percentage difference is Error=$|(\langle H\rangle-E_{3B})/E_{3B}|\times 100\%$. EV | 2A + R
---|---
$\quad\langle H_{0}\rangle\quad$ | $+247.66$
$\langle V_{1}\rangle$ | $+443.39$
$\langle V_{2}\rangle$ | $-825.90$
$\langle V_{3}\rangle$ | $-825.76$
$\langle V\rangle$ | $-1208.27$
$\langle H\rangle$ | $-960.61$
$E_{3B}$ | $-960.58$
Error | $+0.00312$

Fig. 5 illustrates quantitatively both screenings [Eq. (21)] with $l_{0}=10$
Å. In our actual calculations, the results for the trion binding energy will
be obtained by performing the extrapolation to $l_{0}=0$. Before that, the
results of the expectation values for the Rytova-Keldysh potential screened at
low momenta for $l_{0}=100$ Å are depicted in Table 5. Due to the contribution
of the spectator electron external to the exciton, the expectation value of
the kinetic energy is somewhat larger than the one found for the exciton given
in Eq. (19) with a value of 214.64 meV compared to the trion one of 247.66
meV. The expectation values of the attractive potentials $V_{2}$ and $V_{3}$
are somewhat less in magnitude than the one for the exciton of -967.96 meV. In
the trion magnitude of the potential energy of the repulsive potential is
about one-half of the attractive one. This last feature can be understood as
the electrons should be more separated than the relative distance within the
strongly bound electron-hole pair. While the trion and exciton splitting is
207.29 meV, it shows a weakly bound trion with respect to the exciton. Table 5
also indicates a good accuracy found in our solution by comparing the results
from the expectation value of the Hamiltonian $\langle H\rangle$ and the
energy $E_{3B}$ obtained by solving the coupled Faddeev integral equations,
which shows a deviation of only 0.003% between these two values.

Table 6: Trion ground state binding energies ($E_{3B}$) for different screening parameter $l_{0}$ obtained from two screening schemes shown in Fig. 5 and given in Eq. (21). $V(q)\to e^{-l_{0}\,q}V_{ee}(q)$ | $V(q)\to(1-e^{-l_{0}q})V_{ee}(q)$
---|---
$l_{0}$ (Å) |  $E_{3B}$ (meV) |  $l_{0}$ (Å) |  $E_{3B}$ (meV)
$25$ | $-1195.1$ | $1$ | $-1444.8$
$20$ | $-1150.6$ | $5$ | $-1309.6$
$17$ | $-1117.3$ | $10$ | $-1212.0$
$15$ | $-1091.5$ | $15$ | $-1147.7$
$13$ | $-1062.2$ | $20$ | $-1101.2$
$11$ | $-1028.9$ | $30$ | $-1037.3$
$10$ | $-1010.6$ | $50$ | $-965.6$
$9$ | $-991.0$ | $70$ | $-924.8$
$8$ | $-970.3$ | $90$ | $-898.2$
$7$ | $-948.0$ | $100$ | $-888.0$

In Appendix F, we present a convergence study of the trion energy, as
summarized in Table 8, which requires the extrapolation in the number of
quadrature points. Noteworthy that the results presented in Table 5 are not
converged in terms of the number o quadrature points but are good enough to
compute the expectation value of the Hamiltonian, which should be interpreted
as a lower bound. The extrapolated results from Table 8 are collected in Table
6 for the two forms of the screening implemented for the Rytova-Keldysh
electron-electron repulsive potential. Fig. 6 shows this extrapolation as a
function of (left panel) $l_{0}^{-1}$, for the short distance screening trion
energies of the repulsive potential and (right panel) $l_{0}$, for the large
distance screening. As shown in Fig. 6, the results exhibit a perfect linear
behavior which allows an accurate extrapolation to the trion binding energy.
The linear extrapolation on binding energies obtained from the first screening
(left panel) on the domain $[70,100]$ Å-1 leads to a trion binding energy of
$-49.6$ meV, while a linear extrapolation on the second screening (right
panel) on the domain $[7,10]$ Å leads to a trion binding energy of $-49.4$
meV. These results lie in the range of previous experiments reported in Refs.
[30, 29].

Figure 6: Extraction of trion binding energy $E_{t}$ with a linear extrapolation on energies obtained from two screenings (see Table 6) at the physical points $l^{-1}_{0}=0$ Å-1 (left panel) and $l_{0}=0$ Å(right panel). Table 7: The inner product of the Faddeev components $\langle\psi_{i}|\psi_{j}\rangle$ and their contributions in the normalization of the 3B wave function $|\Psi\rangle$ obtained for the screening parameter $l_{0}=100$ Å in the screening scheme $V_{1}(q)\to(1-e^{-l_{0}q})V_{ee}(q)$. | $j=1$ | $j=2$ | $j=3$
---|---|---|---
$i=1$ | $0.1772$ | $-0.2808$ | $-0.2806$
$i=2$ | $-0.2808$ | $0.5197$ | $0.4534$
$i=3$ | $-0.2806$ | $0.4534$ | $0.5190$

Figure 7: The evolution of the Faddeev components (top panels)
$\psi_{1}(p_{1},q_{1},\phi_{1}=0)$ and (bottom panels)
$\psi_{2}(p_{1},q_{1},\phi_{1}=0)$ obtained for the screening parameter
$l_{0}=1,30,50,70,100$ Å in the screening scheme
$V_{1}(q)\to(1-e^{-l_{0}q})V_{ee}(q)$.

It is worth noting that in excitonic physics, the electron-hole interaction
comprises both attractive screened interaction and repulsive exchange
interaction. However, for strongly bound excitons, the exchange interaction
generally has a small contribution, as demonstrated in Ref. [56]. For the
trion, apart from the electron-hole interaction, there is also an exchange
term for the electron-electron interaction. However, considering that the
electron in the trion is weakly bound to the exciton and associated with small
momenta as compared to the reciprocal vector, the exchange terms are
presumably much less significant in determining the trion binding energy as
compared to the contribution to the exciton energy, which is already small in
this case. Despite this, the Faddeev approach to solving the Hamiltonian
eigenvalue problem in momentum space is suitable for dealing with non-local
exchange terms in a similar way for the exciton Hamiltonian [57], which is
left for a future study.

##  VI Rytova-Keldysh Trion Clusterization

Trion structure is studied for the screened electron-electron potential
$V_{1}(q)\to(1-e^{-l_{0}q})V_{ee}(q)$ with $l_{0}=100$ Å. We chose this
particular model since the electron potential is screened at large distances,
which acts together with the natural screening of the exciton interaction with
the spectator electron. In this sense, the two effects act coherently, making
the trion to be overbinding with an energy of $-207.3$ meV compared to the
extrapolated one of $-49.6$ meV. We should keep in mind that features
associated with small trion binding energy, with respect to exciton, will be
further highlighted towards the converged trion with Rytova-Keldysh potential.
Our analysis is based on the screened electron-electron potential, which for
the moment, is a limitation of our numerical method applied to the repulsive
Rytova-Keldysh potential. Despite that, we study the structure of the trion
within the screened model to shed light on its structure and compare it with
the 2A+R Yamaguchi model.

We should emphasize that we consider a negatively charged trion with one hole
and two electrons, where the two electrons will have the same effective mass.
In general, for TMDs, the electrons have the same mass if they belong to the
same band/minimum point. This happens for 1) intravalley electrons with the
same spin or 2) intervalley and opposite spin electrons, however, they will
have the same mass if we neglect the spin-orbit coupling for the conduction
band, which is the case in our work.

Figure 8: The evolution of the total wave function
$\Psi(p_{1},q_{1},\phi_{1}=0)$ (top panels) and $\Psi(p_{2},q_{2},\phi_{2}=0)$
(bottom panels) obtained for the screening parameter $l_{0}=1,30,50,70,100$ Å
in the screening scheme $V_{1}(q)\to(1-e^{-l_{0}q})V_{ee}(q)$.

The overlaps between the Faddeev components of the wave function are given in
Table 7. As expected, the relative normalization of the component
$\langle\psi_{1}|\psi_{1}\rangle$, where the hole is the spectator particle of
the interacting electron-electron pair, is almost three times smaller than
$\langle\psi_{2}|\psi_{2}\rangle=\langle\psi_{3}|\psi_{3}\rangle$. Similarly,
we have also observed it in the 2A+R Yamaguchi model (cf. Table 3), which
shows a quite small overlap $\langle\psi_{1}|\psi_{1}\rangle$ with respect to
the total normalization of the wave function. The configuration where the hole
is a spectator of the interacting electron-electron pair is suppressed,
favoring the clusterization of the wave function where the electron and hole
are close, forming essentially the exciton and a distant spectator electron.
We also observe the opposite sign of $\langle\psi_{1}|\psi_{2,3}\rangle$ with
respect to $\langle\psi_{2}|\psi_{3}\rangle$, as a manifestation of the
repulsive interaction between the electrons. The comparison with the 2A+R
Yamaguchi potentials overlaps from Table 3 with the Rytova-Keldysh results in
Table 7, showing that the suppression of $\psi_{1}$ is much more dramatic for
the former model. The reason for that is twofold: (i) the relatively smaller
difference between the 3B and 2B binding energies from the Yamaguchi model,
namely $(E_{3B}-E_{2B})/E_{3B}=0.169$ compared to 0.216 from the Rytova-
Keldysh screened model, and (ii) the short-range Yamaguchi potential, while
the Rytova-Keldysh potential has a long-range tail. However, the Rytova-
Keldysh trion has a considerably smaller 3B binding energy than the resulting
one for the screened Rytova-Keldysh electron-electron potential, namely
$(E_{3B}-E_{2B})/E_{3B}=0.061$ obtained from the extrapolated value of -802.9
meV in Fig. 6. Therefore, we expect a more evident clusterization of the
exciton within the trion.

Figure 9: The evolution of the total wave function $\Psi(p_{1},q_{1},\phi_{1})$ with respect to the angle $\phi_{1}$ for a fixed screening parameter $l_{0}=100$ Å in the screening scheme $V_{1}(q)\to(1-e^{-l_{0}q})V_{ee}(q)$. |
---|---
Figure 10: Comparison between trion wave function calculated with Rytova-
Keldysh potential (left panel) and re-scaled Yamaguchi potential model (right
panel) for $\phi_{1}=0$.

In Fig. 7, it is shown the results for the Faddeev components
$\psi_{1}(p_{1},q_{1},\phi_{1}=0)$ (top panel) and
$\psi_{2}(p_{2},q_{2},\phi_{2}=0)=\psi_{3}(p_{3},q_{3},\phi_{3}=0)$ (bottom
panel), where the momenta are defined in terms of $p_{1}$ and $q_{1}$
according to Eq. (18a). The calculations were done for $l_{0}=1,30,50,70,100$
Å (from left to right panels) with the electron-electron screened potential
$V_{1}(q)\to(1-e^{-l_{0}q})V_{ee}(q)$. As $l_{0}$ increases, the trion binding
decreases and the electrons are pushed away from each other. At the same time,
one of them has a hole in the vicinity region of the exciton size as expressed
by the dominant configurations, namely $\psi_{2}\sim$ [3(e)1(h)]—2(e) and
$\psi_{3}\sim$ [1(h)2(e)]—3(e) (the square brackets represent the exciton),
while $\psi_{1}(p_{1},q_{1},0)$ just reflects the short-range repulsion,
spreading $q_{1}$ and $p_{1}$ to larger values, as depicted in the top panels
of Fig. 7. For small values of the screening parameter $l_{0}$, in which the
long-range tail of the repulsive potential is damped, we observe, analyzing
the momentum distribution of $\psi_{1}(p_{1},q_{1},0)$, larger trion bindings
tending to be more symmetrical in $p_{1}$ and $q_{1}$. As a consequence, it
implies a geometric configuration where the hole is equally separated from the
two electrons, disfavoring the cluster structure against a more symmetrical
configuration, schematically like e–h–e. The node appears in $\psi_{1}$ at
larger values of $p_{1}$ for large $l_{0}$ values, as the repulsion is
increased and it is intense at the short range.

In the bottom panel of Fig. 7, the evolution of $\psi_{2}(p_{2}({\bf
p}_{1},{\bf q}_{1}),q_{2}({\bf p}_{1},{\bf q}_{1}),\phi_{2}=0)$ for
$\phi_{1}=0$ with $l_{0}$ is shown in the ($p_{1}\times q_{1}$)–plane. As
already observed in the top panel of Fig. 7 for $\psi_{1}$, as $l_{0}$
increases, the configuration migrates from the symmetrical situation, e–h–e,
to a cluster one, [3(e)1(h)]—2(e). The node line becomes evident together with
the cluster structure when the long-range screening is reduced, and the
electron becomes weakly bounded with respect to the exciton. Higher amplitude
values of the wave function are found for $p_{1}\propto q_{1}$, when the
cluster structure dominates, as already discussed for the 2A+R Yamaguchi 3B
model.

The total wave function is presented in Fig. 8 in the $p_{1}\times q_{1}$ (top
panel) and $p_{2}\times q_{2}$ (bottom panel) planes for several values of
$l_{0}$ from 1 to 100 Å with $\phi_{1}=\phi_{2}=0$. The reduction of the
screening at large distances turns the pattern similar to the one observed for
$\psi_{2}$ in the ($p_{1}\times q_{1}$)–plane [cf. bottom panels in Fig. 7],
more evident due to its dominance over $\psi_{1}$, reminding that
$\psi_{2}\equiv\psi_{3}$ from the symmetric configuration of the two
electrons, which have to be in a singlet spin state or an antisymmetric
combination of different valley states. In the bottom panels of Fig. 8,
results for the total wave function are presented in the ($p_{2}\times
q_{2}$)–plane, where again, the more symmetric e–h–e configuration dominates
at the strong trion binding and weaker repulsion between the electrons. By
reducing the screening of the repulsive potential, the electron turns to be
weakly bound to the exciton, and the system presents an evident cluster
structure with the coherent superposition of the two configurations
[3(e)1(h)]—2(e) and [1(h)2(e)]—3(e). As a function of $p_{2}$ and $q_{2}$, the
total wave functions, demonstrated in the top panels of Fig. 8, exhibit two
branches of higher probability density, namely for $p_{2}\propto q_{2}$ and
for $q_{2}$ small with $p_{2}$ spreading in the region shown in the figure.
This behavior was found in the 2A+R Yamaguchi model and is associated with the
cluster structure, already discussed in detail in Sec. IV.

In Fig. 9, the angular dependence in $\phi_{1}$ is explored for $l_{0}=100$ Å.
Note that, at $\phi_{1}=\theta$ or $\phi_{1}=180^{o}-\theta$, the results are
the same due to the symmetry of the wave function by exchanging the momentum
of the electrons. The configuration space wave function is symmetrical by the
exchange of the electrons once the antisymmetry is ensured by the spin state.
The slope of the node line is deformed when $\phi_{1}$ changes between ${\bf
p}_{1}$ and ${\bf q}_{1}$, and becomes more elongated at 90∘. The node line
format is basically kept regardless of the $\phi_{1}$-parameters, which in
turn reflects the situation where the e-h-e system has a hole with small
momentum with respect to the center of mass of the electron pair. It can be
clearly seen that the zero of the wave function in the ($p_{1}\times
q_{1}$)–plane starts at $q_{1}=0$.

Finally, in Fig. 10, we directly compare the trion wave function computed with
both Rytova-Keldysh (left panel) and 2A+R Yamaguchi (right panel) models. For
that, we performed a re-scaling of the Yamaguchi separable potential model to
physical units of the exciton and trion. The results for the exciton and trion
were obtained for the Yamaguchi model in units of $\hbar=m=1$. Turning to
physical units we have that

$E_{2B}=-\frac{\hbar^{2}}{m_{e}}\lambda^{2}\,0.1\quad\text{and}\quad
E_{3B}=-\frac{\hbar^{2}}{m_{e}}\lambda^{2}\,0.12034\,,$

and to obtain the dimensional constant $\lambda$ we use the exciton binding
energy of $E_{2B}=-753$ meV, and $m=(m_{e}+m_{h})/2=0.505\,m_{0}$, which gives
$\lambda=0.706576$ Å-1. The trion binding energy in the 2A+R Yamaguchi model
in meV units is $E_{t}=153.225$ meV, comparable with 207.26 meV from the
regulated repulsive Rytova-Keldysh potential with $l_{0}=100$ Å given in Table
5. After the re-scaling to physical units, the comparison shows essentially
the same structure of the trion wave function resulting from calculations with
the Rytova-Keldysh and 2A+R Yamaguchi potentials. Remarkably the node line has
the same form, and the region where the wave function attains the highest
values is to the left of the node line with $p_{1}\propto q_{1}$. The node
line is shifted to larger values of $p_{1}$ for the 2A+R Yamaguchi model with
respect to the Rytova-Keldysh potential, which should reflect the short-range
nature of the Yamaguchi potential contrasting with the Rytova-Keldysh
potential.

##  VII Summary

Our work can be helpful for theoretical-computational research on trions on
the following main aspects: (i) we provided a general formulation of the
Faddeev equations to compute the wave function for a 2D trion in momentum
space; (ii) we developed a method for the accurate calculation of the trion
binding energy in freestanding monolayer MoS2 with different regularization
schemes for the repulsive electron-electron Rytova-Keldysh potential, with the
numerical precision checked by computing the expectation value of the
Hamiltonian; and (iii) we analyzed in detail the degree of clusterization of
the weakly bound trion state through the momentum distributions of the total
wave function and its Faddeev components.

The repulsive electron-electron Rytova-Keldysh potential posed a numerical
challenge in solving the Faddeev equations. To overcome that, we have devised
two different regularization functions applied only to the repulsive term,
which turns the trion weakly bound with respect to the exciton. The two
different choices for the regularization functions were chosen as (i)
$e^{-l_{0}q}$ ($7<l_{0}<25$ Å), that acts in the high momentum transfer
region, and (ii) $1-e^{-l_{0}q}$ ($1<l_{0}<100$ Å), that acts in the low
momentum transfer region. The results were then extrapolated to $l_{0}^{-1}\to
0$ in the former case and to $l_{0}\to 0$ in the latter one. The extrapolation
results were good in one part in $10^{4}$, resulting in a prediction of the
trion energy of $-49.3(1)$ meV for monolayer MoS2. Our result lies in the
range of the experimental results for suspended samples 44-80 meV as reported
in Refs. [30, 29]. It can not be compared with the results where the MoS2 is
deposited on a substrate, where the trion binding energy lies in the range of
$20-43$ meV [25, 26, 27, 28], since the interaction becomes weaker due to
screening, affecting both the exciton and trion complexes. The value of
$-49.3(1)$ meV is within the previous numerical results and particularly close
to the converged negatively charged intralayer trion binding energy computed
within an ab-initio many-body theory, which was found to be $58$ meV [38] with
the exciton binding energy of $-0.76$ eV.

Furthermore, we have analyzed the structure of the trion wave function by
decomposing it into its Faddeev components for both the Rytova-Keldysh and
Yamaguchi potentials. Despite the fact that the two interactions have
different large distance tails, we have observed qualitative similarities in
the wave functions. Both trion models showed a remarkable dominance of the
[eh]–e configurations, corresponding to the Faddeev components of the wave
function where the electron is a spectator. We should observe that the
electrons are in a symmetric configuration by exchanging their spatial
coordinates and are considered to be in a spin singlet state or an
antisymmetric combination of different valley states.

We expect that without regulating the repulsive Rytova-Keldysh potential, the
dominance of the strong cluster structure [eh]–e would be more evident, as the
trion will be even more weakly bound, as indicated by the extrapolated
results. Therefore, our study suggests that a realistic Rytova-Keldysh
potential calculation, with the small relative exciton and trion splitting,
can profit from the cluster structure and use it to build more accurate
methods, eventually relying on the dominant exciton structure together with a
simplified potential that contains the low-energy electron-electron continuum
information.

###### Acknowledgements.

We thank Charlotte Elster for the helpful discussion on the Coulomb
interaction screening. This work is a part of the project INCT-FNA proc. No.
464898/2014-5. The work of M. R. H. was supported by the National Science
Foundation under Grant No. NSF-PHY-2000029 with Central State University. M.
R. H. also thanks the Ohio Supercomputer Center (OSC) for the use of their
facilities under Grant No. POS0104. K. M. acknowledges a Ph.D. scholarship
from the Brazilian agency CNPq (Conselho Nacional de Desenvolvimento
Científico e Tecnológico). K. M., T. F., A. J. C., and D. R. C. were supported
by CNPq Grant No. 400789/2019-0, 308486/2015-3, 315408/2021-9, and
313211/2021-3, respectively. A. J. C. and T. F. acknowledge Fundação de Amparo
à Pesquisa do Estado de São Paulo (FAPESP) under Grant No. 2022/08086-0 and
Thematic Projects 2017/05660-0 and 2019/07767-1, respectively.

##  Appendix A Momentum space representation of Faddeev equations

The Faddeev components, i.e., the projections of the coupled Faddeev equations
(3) on the 3B basis states $|{\bf p}_{i}{\bf q}_{i}\rangle$, can be written as

$\langle{\bf p}_{i}{\bf
q}_{i}|\psi_{i}\rangle=G_{0}^{(i,jk)}(E,p_{i},q_{i})\langle{\bf p}_{i}{\bf
q}_{i}|t_{i}\Big{[}|\psi_{j}\rangle+|\psi_{k}\rangle\Big{]}\,.$ (22)

By inserting the completeness relation of Eq. (7) into Eq. (22), it leads to

$\displaystyle\langle{\bf p}_{i}{\bf q}_{i}|\psi_{i}\rangle$ $\displaystyle=$
$\displaystyle G_{0}^{(i,jk)}(E,p_{i},q_{i})\int d^{2}p^{\prime}_{i}\int
d^{2}q^{\prime}_{i}$ (23) $\displaystyle\times$ $\displaystyle\langle{\bf
p}_{i}{\bf q}_{i}|t_{i}|{\bf p}^{\prime}_{i}{\bf
q}^{\prime}_{i}\rangle\bigl{[}\langle{\bf p}^{\prime}_{i}{\bf
q}^{\prime}_{i}|\psi_{j}\rangle+\langle{\bf p}^{\prime}_{i}{\bf
q}^{\prime}_{i}|\psi_{k}\rangle\bigr{]},\quad$ (24)

where

$\displaystyle
G_{0}^{(i,jk)}(E,p_{i},q_{i})=\dfrac{1}{E-\dfrac{p_{i}^{2}}{2\mu_{jk}}-\dfrac{q_{i}^{2}}{2\mu_{i,jk}}},$
(25)

and

$\displaystyle\langle{\bf p}_{i}{\bf q}_{i}|t_{i}|{\bf p}_{i}^{\prime}{\bf
q}_{i}^{\prime}\rangle$ $\displaystyle=$ $\displaystyle\delta({\bf q}_{i}-{\bf
q}^{\prime}_{i})\langle{\bf p}_{i}|t_{i}|{\bf p}_{i}^{\prime}\rangle.$ (26)

In order to evaluate Eq. (24), we need to compute $\langle{\bf
p}^{\prime}_{i}{\bf q}^{\prime}_{i}|\psi_{j}\rangle$ and $\langle{\bf
p}^{\prime}_{i}{\bf q}^{\prime}_{i}|\psi_{k}\rangle$, where by inserting a
completeness relation, one obtains

$\displaystyle\langle{\bf p}^{\prime}_{i}{\bf q}^{\prime}_{i}|\psi_{j}\rangle$
$\displaystyle=$ $\displaystyle\int d^{2}p^{\prime\prime}_{j}\int
d^{2}q^{\prime\prime}_{j}\langle{\bf p}_{i}^{\prime}{\bf q}_{i}^{\prime}|{\bf
p}^{\prime\prime}_{j}{\bf q}^{\prime\prime}_{j}\rangle\langle{\bf
p}^{\prime\prime}_{j}{\bf q}^{\prime\prime}_{j}|\psi_{j}\rangle$ (27)
$\displaystyle=$ $\displaystyle\langle{\bf p}_{j}({\bf p}^{\prime}_{i},{\bf
q}^{\prime}_{i}),{\bf q}_{j}({\bf p}^{\prime}_{i},{\bf
q}^{\prime}_{i})|\psi_{j}\rangle,$ (28) $\displaystyle\langle{\bf
p}^{\prime}_{i}{\bf q}^{\prime}_{i}|\psi_{k}\rangle$ $\displaystyle=$
$\displaystyle\int d^{2}p^{\prime\prime}_{k}\int
d^{2}q^{\prime\prime}_{k}\langle{\bf p}_{i}^{\prime}{\bf q}_{i}^{\prime}|{\bf
p}^{\prime\prime}_{k}{\bf q}^{\prime\prime}_{k}\rangle\langle{\bf
p}^{\prime\prime}_{k}{\bf q}^{\prime\prime}_{k}|\psi_{k}\rangle$ (29)
$\displaystyle=$ $\displaystyle\langle{\bf p}_{k}({\bf p}^{\prime}_{i},{\bf
q}^{\prime}_{i}),{\bf q}_{k}({\bf p}^{\prime}_{i},{\bf
q}^{\prime}_{i})|\psi_{k}\rangle,$ (30)

where the relation between different Jacobi momenta is given by

$\displaystyle{\bf p}_{i}({\bf p}_{j},{\bf q}_{j})$ $\displaystyle\equiv$
$\displaystyle\pmb{\mathcal{P}}_{ij}({\bf p}_{j},{\bf q}_{j})=\alpha_{ij}{\bf
p}_{j}+\beta_{ij}{\bf q}_{j},$ (31) $\displaystyle{\bf q}_{i}({\bf p}_{j},{\bf
q}_{j})$ $\displaystyle\equiv$ $\displaystyle\pmb{\mathcal{Q}}_{ij}({\bf
p}_{j},{\bf q}_{j})=\gamma_{ij}{\bf p}_{j}+\eta_{ij}{\bf q}_{j},$ (32)

with

$\displaystyle\alpha_{ij}=-\dfrac{m_{j}}{m_{jk}},\,\beta_{ij}=\mathcal{E}_{ij}\dfrac{m_{k}\,m_{ijk}}{m_{ik}m_{jk}},$
(33)
$\displaystyle\gamma_{ij}=-\mathcal{E}_{ij},\,\eta_{ij}=-\dfrac{m_{i}}{m_{ik}},$
(34) $\displaystyle m_{ij}=m_{i}+m_{j},\quad m_{ijk}=m_{i}+m_{j}+m_{k},$ (35)
$\displaystyle\mathcal{E}_{ij}=\begin{cases}1&\text{for cyclic permutation}\\\
-1&\text{for anti-cyclic permutation}\end{cases}.$ (36)

By using Eqs. (26) and (30) we can rewrite Eq. (24) as

$\displaystyle\psi_{i}({\bf p}_{i},{\bf q}_{i})$ $\displaystyle=$
$\displaystyle G_{0}^{(i,jk)}(E,p_{i},q_{i})\int
d^{2}p^{\prime}_{i}\,t_{i}({\bf p}_{i},{\bf p}_{i}^{\prime})$ (37)
$\displaystyle\times$
$\displaystyle\Biggl{[}\psi_{j}\biggl{(}\pmb{\mathcal{P}}_{ji}({\bf
p}^{\prime}_{i},{\bf q}_{i}),\pmb{\mathcal{Q}}_{ji}({\bf p}^{\prime}_{i},{\bf
q}_{i})\biggr{)}$ (39)
$\displaystyle+\psi_{k}\biggl{(}\pmb{\mathcal{P}}_{ki}({\bf
p}^{\prime}_{i},{\bf q}_{i}),\pmb{\mathcal{Q}}_{ki}({\bf p}^{\prime}_{i},{\bf
q}_{i})\biggr{)}\Biggr{]}.$

To solve the coupled 2D Faddeev integral equations, i.e. Eq. (39), as shown in
Fig. 11, we choose a coordinate system where ${\bf p}_{i}$ is parallel to the
$x-$axis, ${\bf p}^{\prime}_{i}$ and ${\bf q}_{i}$ are free in the 2D space.

$x$${\bf p}_{i}$$y$${\bf p}^{\prime}_{i}$${\bf
q}_{i}$$\phi^{\prime}_{i}$$\phi_{i}$$\phi_{q_{i},p^{\prime}_{i}}$ Figure 11:
The coordinate system for the solution of the coupled Faddeev integral
equations (39).

So the angle variables will be as

$\displaystyle(\hat{{\bf q}_{i}},\hat{{\bf p}_{i}})$ $\displaystyle=$
$\displaystyle\phi_{i},$ (40) $\displaystyle(\hat{{\bf
p}^{\prime}_{i}},\hat{{\bf p}_{i}})$ $\displaystyle=$
$\displaystyle\phi^{\prime}_{i},$ (41) $\displaystyle(\hat{{\bf
q}_{i}},\hat{{\bf p}^{\prime}_{i}})$ $\displaystyle=$
$\displaystyle\phi_{q_{i},p^{\prime}_{i}}=\phi_{i}-\phi^{\prime}_{i}.$ (42)

The shifted momentum and angle variables are defined as

$\displaystyle\mathcal{P}_{ij}$ $\displaystyle\equiv$
$\displaystyle\mathcal{P}_{ij}({\bf p}^{\prime}_{j},{\bf
q}_{j})=|\alpha_{ij}{\bf p}^{\prime}_{j}+\beta_{ij}{\bf
q}_{j}|=\sqrt{(\mathcal{P}_{ij}^{X})^{2}+(\mathcal{P}_{ij}^{Y})^{2}},$ (43)
$\displaystyle\mathcal{Q}_{ij}$ $\displaystyle\equiv$
$\displaystyle\mathcal{Q}_{ij}({\bf p}^{\prime}_{j},{\bf
q}_{j})=|\gamma_{ij}{\bf p}^{\prime}_{j}+\eta_{ij}{\bf
q}_{j}|=\sqrt{(\mathcal{Q}_{ij}^{X})^{2}+(\mathcal{Q}_{ij}^{Y})^{2}},$ (44)
$\displaystyle\phi_{ij}$ $\displaystyle=$
$\displaystyle(\hat{\mathcal{P}_{ij}},\hat{\mathcal{Q}_{ij}})={\text{atan}2}(det,dot),\quad
0<\phi_{ij}<2\pi,$ (45)

where

$\displaystyle\mathcal{P}_{ij}^{X}$ $\displaystyle=$
$\displaystyle\alpha_{ij}p^{\prime}_{j}\cos(\phi^{\prime}_{j})+\beta_{ij}q_{j}\cos(\phi_{j}),$
(46) $\displaystyle\mathcal{P}_{ij}^{Y}$ $\displaystyle=$
$\displaystyle\alpha_{ij}p^{\prime}_{j}\sin(\phi^{\prime}_{j})+\beta_{ij}q_{j}\sin(\phi_{j}),$
(47) $\displaystyle\mathcal{Q}_{ij}^{X}$ $\displaystyle=$
$\displaystyle\gamma_{ij}p^{\prime}_{j}\cos(\phi^{\prime}_{j})+\eta_{ij}q_{j}\cos(\phi_{j}),$
(48) $\displaystyle\mathcal{Q}_{ij}^{Y}$ $\displaystyle=$
$\displaystyle\gamma_{ij}p^{\prime}_{j}\sin(\phi^{\prime}_{j})+\eta_{ij}q_{j}\sin(\phi_{j}),$
(49)

and

$\displaystyle\begin{cases}det=\mathcal{P}_{ij}^{X}\cdot\mathcal{Q}_{ij}^{Y}-\mathcal{P}_{ij}^{Y}\cdot\mathcal{Q}_{ij}^{X},\\\
dot=\mathcal{P}_{ij}^{X}\cdot\mathcal{Q}_{ij}^{X}+\mathcal{P}_{ij}^{Y}\cdot\mathcal{Q}_{ij}^{Y}.\end{cases}$
(50)

##  Appendix B 3B wave functions in momentum space

3B wave function is given as

$\displaystyle\Psi=\sum_{i=1}^{3}\psi_{i},$ (51)

where

$\langle{\bf p}_{i}{\bf q}_{i}|\Psi\rangle=\langle{\bf p}_{i}{\bf
q}_{i}|\psi_{i}\rangle+\langle{\bf p}_{i}{\bf
q}_{i}|\psi_{j}\rangle+\langle{\bf p}_{i}{\bf q}_{i}|\psi_{k}\rangle.$ (52)

Using Eq. (30) and the coordinate system defined in Fig. 11, one has that

$\displaystyle\Psi(p_{i},q_{i},\phi_{i})=\psi_{i}(p_{i},q_{i},\phi_{i})$ (53)
$\displaystyle+\psi_{j}\biggl{(}\mathcal{P}_{ji}(p_{i},q_{i},\phi_{i}),\mathcal{Q}_{ji}(p_{i},q_{i},\phi_{i}),\phi_{ji}(p_{i},q_{i},\phi_{i})\biggr{)}$
(54)
$\displaystyle+\psi_{k}\biggl{(}\mathcal{P}_{ki}(p_{i},q_{i},\phi_{i}),\mathcal{Q}_{ki}(p_{i},q_{i},\phi_{i}),\phi_{ki}(p_{i},q_{i},\phi_{i})\biggr{)}.\quad$
(55)

The 3B wave function is normalized as

$\displaystyle\langle\Psi|\Psi\rangle$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{3}\langle\Psi|\psi_{i}\rangle$ (56) $\displaystyle=$
$\displaystyle
2\pi\int_{0}^{\infty}dp_{i}\,p_{i}\int_{0}^{\infty}dq_{i}\,q_{i}\int_{0}^{2\pi}d\phi_{i}\,\Psi^{2}(p_{i},q_{i},\phi_{i})$
(57) $\displaystyle=$ $\displaystyle
2\pi\int_{0}^{\infty}dp_{i}\,p_{i}\int_{0}^{\infty}dq_{i}\,q_{i}\int_{0}^{2\pi}d\phi_{i}\,\Psi(p_{i},q_{i},\phi_{i})$
(58) $\displaystyle\times$
$\displaystyle\biggl{[}\psi_{i}(p_{i},q_{i},\phi_{i})$ (61)
$\displaystyle+\psi_{j}(\mathcal{P}_{ji}({\bf p}_{i},{\bf
q}_{i}),\mathcal{Q}_{ji}({\bf p}_{i},{\bf q}_{i}),\phi_{ji})$
$\displaystyle+\psi_{k}(\mathcal{P}_{ki}({\bf p}_{i},{\bf
q}_{i}),\mathcal{Q}_{ki}({\bf p}_{i},{\bf q}_{i}),\phi_{ki})\biggr{]}=1.$

Table 8: The convergence of trion ground state binding energies (in meV) as a function of the number of mesh points for Jacobi momenta $N_{p}=N_{q}$ obtained for different values of screening parameter $l_{0}$ for screening electron-electron interactions $V(q)\to e^{-l_{0}\,q}V_{ee}(q)$ (upper panel) and $V(q)\to(1-e^{-l_{0}q})V_{ee}(q)$ (lower panel). The last row of each panel shows the extrapolation of trion energy eigenvalues to an infinite number of mesh points.  | $l_{0}(\AA)$
---|---
$N_{p}=N_{q}$ | $25$ | $20$ | $17$ | $15$ | $13$ | $11$ | $10$ | $9$ | $8$ | $7$
$200$ | $-1275.6$ | $-1231.1$ | $-1197.7$ | $-1171.8$ | $-1142.5$ | $-1109.0$ | $-1090.5$ | $-1070.8$ | $-1049.8$ | $-1027.6$
$250$ | $-1258.7$ | $-1214.2$ | $-1180.8$ | $-1155.0$ | $-1125.6$ | $-1092.2$ | $-1073.7$ | $-1054.1$ | $-1033.0$ | $-1011.5$
$300$ | $-1248.3$ | $-1203.8$ | $-1170.4$ | $-1144.6$ | $-1115.3$ | $-1081.8$ | $-1063.4$ | $-1043.7$ | $-1022.8$ | $-1000.7$
$350$ | $-1241.3$ | $-1196.8$ | $-1163.5$ | $-1137.7$ | $-1108.3$ | $-1074.9$ | $-1056.5$ | $-1036.8$ | $-1016.0$ | $-993.7$
$N_{p},N_{q}\to\infty$ | $-1195.1$ | $-1150.6$ | $-1117.3$ | $-1091.5$ | $-1062.2$ | $-1028.9$ | $-1010.6$ | $-991.0$ | $-970.3$ | $-948.0$
| $l_{0}(\AA)$
$N_{p}=N_{q}$ | $1$ | $5$ | $10$ | $15$ | $20$ | $30$ | $50$ | $70$ | $90$ | $100$
$200$ | $-1588.7$ | $-1447.8$ | $-1346.2$ | $-1279.4$ | $-1231.0$ | $-1165.0$ | $-1090.8$ | $-1049.3$ | $-1023.3$ | $-1014.8$
$250$ | $-1558.6$ | $-1418.8$ | $-1318.1$ | $-1251.8$ | $-1203.8$ | $-1138.2$ | $-1064.7$ | $-1023.6$ | $-998.3$ | $-988.8$
$300$ | $-1540.0$ | $-1400.9$ | $-1300.7$ | $-1234.8$ | $-1187.0$ | $-1121.7$ | $-1048.4$ | $-1007.1$ | $-981.1$ | $-972.3$
$350$ | $-1527.5$ | $-1388.9$ | $-1289.1$ | $-1223.3$ | $-1175.7$ | $-1110.6$ | $-1037.5$ | $-996.3$ | $-970.0$ | $-960.6$
$N_{p},N_{q}\to\infty$ | $-1444.8$ | $-1309.6$ | $-1212.0$ | $-1147.7$ | $-1101.2$ | $-1037.3$ | $-965.6$ | $-924.8$ | $-898.2$ | $-888.0$

##  Appendix C Expectation value of 3B Hamiltonian

By having the 3B wave function, the expectation value of 3B Hamiltonian
$\langle H\rangle$ can be obtained as

$\langle\Psi|H|\Psi\rangle=\langle\Psi|H_{0}|\Psi\rangle+\langle\Psi|V|\Psi\rangle,$
(62)

where

$\displaystyle\langle\Psi|H_{0}|\Psi\rangle=2\pi\int_{0}^{\infty}dp_{i}\,p_{i}\int_{0}^{\infty}dq_{i}\,q_{i}\left(\dfrac{p_{i}^{2}}{2\mu_{jk}}+\dfrac{q_{i}^{2}}{2\mu_{i,jk}}\right)$
(63)
$\displaystyle\times\int_{0}^{2\pi}d\phi_{i}\,\Psi^{2}(p_{i},q_{i},\phi_{i}).$
(64)

and

$\displaystyle\langle\Psi|V|\Psi\rangle=\sum_{i=1}^{3}\langle\Psi|V_{i}|\Psi\rangle$
(65)
$\displaystyle=2\pi\sum_{i=1}^{3}\int_{0}^{\infty}dp_{i}p_{i}\,\int_{0}^{\infty}dq_{i}q_{i}\,\int_{0}^{2\pi}d\phi_{i}\,\int_{0}^{\infty}dp^{\prime}_{i}\,p^{\prime}_{i}\,\int_{0}^{2\pi}d\phi^{\prime}_{i}$
(66)
$\displaystyle\times\Psi(p_{i},q_{i},\phi_{i})\,V_{i}(p_{i},p^{\prime}_{i},\phi^{\prime}_{i})\,\Psi(p^{\prime}_{i},q_{i},\phi_{i}-\phi^{\prime}_{i}).$
(67)

The matrix elements of non-PW potentials can be obtained from the summation of
PW components as

$V(p_{i},p^{\prime}_{i},\phi^{\prime}_{i})=\frac{1}{2\pi}\,\sum^{\infty}_{m=0}\varepsilon_{m}cos(m\phi^{\prime}_{i})\,V_{m}(p_{i},p^{\prime}_{i}).$
(68)

##  Appendix D Verification of the 3B Schrödinger equation in 2D

The Schrödinger equation for the bound state of three particles is given by

$E_{t}|\Psi\rangle=H|\Psi\rangle=(H_{0}+V_{i}+V_{j}+V_{k})|\Psi\rangle.$ (69)

Using the three different sets of Jacobi momenta in momentum space, we obtain

$\displaystyle E_{t}\Psi({\bf p}_{i},{\bf
q}_{i})=\left[\dfrac{p_{i}^{2}}{2\mu_{jk}}+\dfrac{q_{i}^{2}}{2\mu_{i,jk}}\right]\,\Psi({\bf
p}_{i},{\bf q}_{i})$ (70) $\displaystyle+\int d^{2}p^{\prime}_{i}\,V_{i}({\bf
p}_{i},{\bf p}^{\prime}_{i})\,\Psi({\bf p}^{\prime}_{i},{\bf q}_{i})$ (71)
$\displaystyle+\int d^{2}p^{\prime}_{j}\,V_{j}(\pmb{\mathcal{P}}_{ji}({\bf
p}_{i},{\bf q}_{i}),{\bf p}^{\prime}_{j})\,\Psi({\bf
p}^{\prime}_{j},\pmb{\mathcal{Q}}_{ji}({\bf p}_{i},{\bf q}_{i}))$ (72)
$\displaystyle+\int d^{2}p^{\prime}_{k}\,V_{k}(\pmb{\mathcal{P}}_{ki}({\bf
p}_{i},{\bf q}_{i}),{\bf p}^{\prime}_{k})\,\Psi({\bf
p}^{\prime}_{k},\pmb{\mathcal{Q}}_{ki}({\bf p}_{i},{\bf q}_{i})).\quad$ (73)

By using the coordinate system defined in Fig. 11, Eq. (73) can be written as

$\displaystyle
E_{t}\Psi(p_{i},q_{i},\phi_{i})=\left[\dfrac{p_{i}^{2}}{2\mu_{jk}}+\dfrac{q_{i}^{2}}{2\mu_{i,jk}}\right]\,\Psi(p_{i},q_{i},\phi_{i})$
(74)
$\displaystyle+\int_{0}^{\infty}dp^{\prime}_{i}\,p^{\prime}_{i}\int_{0}^{2\pi}d\phi^{\prime}_{i}\,V_{i}(p_{i},p^{\prime}_{i},\phi^{\prime}_{i})\,\Psi(p^{\prime}_{i},q_{i},\phi_{i}-\phi^{\prime}_{i})$
(75)
$\displaystyle+\int_{0}^{\infty}dp^{\prime}_{j}\,p^{\prime}_{j}\int_{0}^{2\pi}\,d\phi^{\prime}_{j}\,V_{j}(\mathcal{P}_{ji}(p_{i},q_{i},\phi_{i}),p^{\prime}_{j},\phi_{\mathcal{P}_{ji},p^{\prime}_{j}})$
(76)
$\displaystyle\quad\times\Psi(p^{\prime}_{j},\mathcal{Q}_{ji}(p_{i},q_{i},\phi_{i}),\phi_{\mathcal{Q}_{ji},p^{\prime}_{j}})$
(77)
$\displaystyle+\int_{0}^{\infty}dp^{\prime}_{k}\,p^{\prime}_{k}\,\int_{0}^{2\pi}\,d\phi^{\prime}_{k}\,V_{k}(\mathcal{P}_{ki}(p_{i},q_{i},\phi_{i}),p^{\prime}_{k},\phi_{\mathcal{P}_{ki},p^{\prime}_{k}})$
(78)
$\displaystyle\quad\times\Psi(p^{\prime}_{k},\mathcal{Q}_{ki}(p_{i},q_{i},\phi_{i}),\phi_{\mathcal{Q}_{ki},p^{\prime}_{k}}).$
(79)

##  Appendix E Numerical Methods

The coupled Faddeev integral equations (8) have an eigenvalue equation form of
$\lambda\ \psi={\cal K}(E)\cdot\psi$ with the eigenvalue $\lambda=1$ and an
eigenvector composed of three Faddeev components
$\psi=\left(\begin{array}[]{l}\psi_{i}\\\ \psi_{j}\\\
\psi_{k}\end{array}\right).$ We solve the eigenvalue equation with the Lanczos
iterative method, which is successfully implemented in two-, three-, and four-
body bound state calculations [58, 59, 60, 61, 62, 63, 64]. Details of the
implementation of this Lanczos technique are discussed in Appendix C2 of Ref.
[65].

We start the iteration process with an initial Gaussian guess for Faddeev
components and stop it after 10-15 iterations. As the kernel of the eigenvalue
equation ${\cal K}(E)$ is energy dependent, the solution of the eigenvalue
equation can be started with an initial guess for the 3B binding energy, and
the search in the binding energy is stopped when $|\lambda-1|\leq 10^{-6}$. To
discretize the continuous momentum and angle variables, we use the Gauss-
Legendre quadratures with a linear mapping $\phi=\pi(1+x)$ for angle variables
and a hyperbolic mapping $p=\frac{1+x}{1-x}$ for the magnitude of Jacobi
momenta.

A typical number of mesh points for angle variables is 60, and for the
magnitude of Jacobi momenta is 300. The solution of coupled Faddeev integral
equations demands a huge number of 3B interpolations on the Faddeev components
$\psi_{j}(\mathcal{P}_{ji},\mathcal{Q}_{ji},\phi_{ji})$ and
$\psi_{k}(\mathcal{P}_{ki},\mathcal{Q}_{ki},\phi_{ki})$ for shifted momentum
and angle variables in each iteration step. We use the Cubic Hermite spline
interpolation of Ref. [66] for its high computational speed and accuracy. To
avoid extrapolations outside the Gauss-Legendre grids, we add an extra point
$0$ to all Jacobi momenta grids and two extra points $0$ and $2\pi$ to all
angle grids.

##  Appendix F 3B energy eigenvalues

In Table 8, we provide our numerical results for 3B energy eigenvalues
obtained from the solution of the coupled Faddeev integral equations (8), for
the Rytova-Keldysh potential given in Eq. (19), with two screening schemes for
electron-electron interactions shown in Fig. 5 and given in Eq. (21), for
different values of screening parameter $l_{0}$ as a function of the number of
mesh points for Jacobi momenta $N_{p}=N_{q}$.

## References

  * Frenkel [1931] J. Frenkel, On the transformation of light into heat in solids. i, Phys. Rev. 37, 17 (1931).
  * Fox and Hexter [1964] D. Fox and R. M. Hexter, Crystal shape dependence of exciton states in molecular crystals, The Journal of Chemical Physics 41, 1125 (1964), https://doi.org/10.1063/1.1726016 .
  * Gross [1956] E. Gross, Optical spectrum of excitons in the crystal lattice, Il Nuovo Cimento (1955-1965) 3, 672 (1956).
  * Knox [1983] R. Knox, Introduction to exciton physics, in  _Collective Excitations in Solids_ (Springer, 1983) pp. 183–245.
  * Dresselhaus [1957] G. Dresselhaus, Absorption coefficients for exciton absorption lines, Phys. Rev. 106, 76 (1957).
  * Elliott [1957] R. J. Elliott, Intensity of optical absorption by excitons, Phys. Rev. 108, 1384 (1957).
  * Lampert [1958] M. A. Lampert, Mobile and immobile effective-mass-particle complexes in nonmetallic solids, Phys. Rev. Lett. 1, 450 (1958).
  * Thomas and Rice [1977] G. Thomas and T. Rice, Trions, molecules and excitons above the mott density in ge, Solid State Communications 23, 359 (1977).
  * Kheng  _et al._ [1993] K. Kheng, R. T. Cox, M. Y. d’ Aubigné, F. Bassani, K. Saminadayar, and S. Tatarenko, Observation of negatively charged excitons ${\mathit{x}}^{\mathrm{-}}$ in semiconductor quantum wells, Phys. Rev. Lett. 71, 1752 (1993).
  * Stébé and Ainane [1989] B. Stébé and A. Ainane, Ground state energy and optical absorption of excitonic trions in two dimensional semiconductors, Superlattices and Microstructures 5, 545 (1989).
  * Pethick and Smith [2008] C. J. Pethick and H. Smith,  _Bose–Einstein condensation in dilute gases_ (Cambridge university press, 2008).
  * Rosa  _et al._ [2022] D. S. Rosa, T. Frederico, G. Krein, and M. T. Yamashita, $d$-dimensional three-body bound-state problem with zero-range interactions, Phys. Rev. A 106, 023311 (2022).
  * Song  _et al._ [2013] X. Song, J. Hu, and H. Zeng, Two-dimensional semiconductors: recent progress and future perspectives, Journal of Materials Chemistry C 1, 2952 (2013).
  * Ugeda  _et al._ [2014] M. M. Ugeda, A. J. Bradley, S.-F. Shi, F. H. da Jornada, Y. Zhang, D. Y. Qiu, W. Ruan, S.-K. Mo, Z. Hussain, Z.-X. Shen, F. Wang, S. G. Louie, and M. F. Crommie, Giant bandgap renormalization and excitonic effects in a monolayer transition metal dichalcogenide semiconductor, Nature Materials 13, 1091 (2014).
  * Mak  _et al._ [2013] K. F. Mak, K. He, C. Lee, G. H. Lee, J. Hone, T. F. Heinz, and J. Shan, Tightly bound trions in monolayer mos2, Nature Materials 12, 207 (2013).
  * Chernikov  _et al._ [2014] A. Chernikov, T. C. Berkelbach, H. M. Hill, A. Rigosi, Y. Li, B. Aslan, D. R. Reichman, M. S. Hybertsen, and T. F. Heinz, Exciton binding energy and nonhydrogenic rydberg series in monolayer ${\mathrm{ws}}_{2}$, Phys. Rev. Lett. 113, 076802 (2014).
  * Chaves  _et al._ [2020] A. Chaves, J. G. Azadani, H. Alsalman, D. R. da Costa, R. Frisenda, A. J. Chaves, S. H. Song, Y. D. Kim, D. He, J. Zhou, A. Castellanos-Gomez, F. M. Peeters, Z. Liu, C. L. Hinkle, S.-H. Oh, P. D. Ye, S. J. Koester, Y. H. Lee, P. Avouris, X. Wang, and T. Low, Bandgap engineering of two-dimensional semiconductor materials, npj 2D Materials and Applications 4, 29 (2020).
  * Cudazzo  _et al._ [2011] P. Cudazzo, I. V. Tokatly, and A. Rubio, Dielectric screening in two-dimensional insulators: Implications for excitonic and impurity states in graphane, Phys. Rev. B 84, 085406 (2011).
  * Liu  _et al._ [2015] X. Liu, T. Galfsky, Z. Sun, F. Xia, E.-c. Lin, Y.-H. Lee, S. Kéna-Cohen, and V. M. Menon, Strong light–matter coupling in two-dimensional atomic crystals, Nature Photonics 9, 30 (2015).
  * Epstein  _et al._ [2020] I. Epstein, A. J. Chaves, D. A. Rhodes, B. Frank, K. Watanabe, T. Taniguchi, H. Giessen, J. C. Hone, N. M. R. Peres, and F. H. L. Koppens, Highly confined in-plane propagating exciton-polaritons on monolayer semiconductors, 2D Materials 7, 035031 (2020).
  * Emmanuele  _et al._ [2020] R. P. A. Emmanuele, M. Sich, O. Kyriienko, V. Shahnazaryan, F. Withers, A. Catanzaro, P. M. Walker, F. A. Benimetskiy, M. S. Skolnick, A. I. Tartakovskii, I. A. Shelykh, and D. N. Krizhanovskii, Highly nonlinear trion-polaritons in a monolayer semiconductor, Nature Communications 11, 3589 (2020).
  * Chang  _et al._ [2018] Y.-C. Chang, S.-Y. Shiau, and M. Combescot, Crossover from trion-hole complex to exciton-polaron in $n$-doped two-dimensional semiconductor quantum wells, Phys. Rev. B 98, 235203 (2018).
  * Žutić  _et al._ [2019] I. Žutić, A. Matos-Abiague, B. Scharf, H. Dery, and K. Belashchenko, Proximitized materials, Materials Today 22, 85 (2019).
  * Seyler  _et al._ [2018] K. L. Seyler, D. Zhong, B. Huang, X. Linpeng, N. P. Wilson, T. Taniguchi, K. Watanabe, W. Yao, D. Xiao, M. A. McGuire, K.-M. C. Fu, and X. Xu, Valley manipulation by optically tuning the magnetic proximity effect in wse2/cri3 heterostructures, Nano Letters 18, 3823 (2018).
  * Mak  _et al._ [2012] K. F. Mak, K. He, C. Lee, G. H. Lee, J. Hone, T. F. Heinz, and J. Shan, Electrical control of neutral and charged excitons in a monolayer semiconductor, Nature Materials 12, 207 (2012).
  * Ross  _et al._ [2013] J. S. Ross, S. Wu, H. Yu, N. J. Ghimire, A. M. Jones, G. Aivazian, J. Yan, D. G. Mandrus, D. Xiao, W. Yao, and X. Xu, Electrical control of neutral and charged excitons in a monolayer semiconductor, Nature Communications 4, 1474 (2013).
  * Soklaski  _et al._ [2014] R. Soklaski, Y. Liang, and L. Yang, Temperature effect on optical spectra of monolayer molybdenum disulfide, Applied Physics Letters 104, 193110 (2014).
  * Zhang  _et al._ [2015] Y. Zhang, H. Li, H. Wang, R. Liu, S.-L. Zhang, and Z.-J. Qiu, On valence-band splitting in layered mos2, ACS Nano 9, 8514 (2015), pMID: 26222731, https://doi.org/10.1021/acsnano.5b03505 .
  * Lin  _et al._ [2019] T.-Y. Lin, Y.-C. Lee, Y.-W. Lee, S.-W. Chang, D.-L. Ma, B.-C. Lin, and H.-L. Chen, Air gap-based cavities dramatically enhance the true intrinsic spectral signals of suspended and pristine two-dimensional materials, The Journal of Physical Chemistry C 123, 5667 (2019).
  * Lin  _et al._ [2014] Y. Lin, X. Ling, L. Yu, S. Huang, A. L. Hsu, Y.-H. Lee, J. Kong, M. S. Dresselhaus, and T. Palacios, Dielectric screening of excitons and trions in single-layer mos2, Nano letters 14, 5569 (2014).
  * Jadczak  _et al._ [2017] J. Jadczak, A. Delgado, L. Bryja, Y. S. Huang, and P. Hawrylak, Robust high-temperature trion emission in monolayers of $\mathrm{Mo}{({\mathrm{S}}_{y}{\mathrm{Se}}_{1-y})}_{2}$ alloys, Phys. Rev. B 95, 195427 (2017).
  * Berkelbach  _et al._ [2013] T. C. Berkelbach, M. S. Hybertsen, and D. R. Reichman, Theory of neutral and charged excitons in monolayer transition metal dichalcogenides, Phys. Rev. B 88, 045318 (2013).
  * Ramirez-Torres  _et al._ [2014] A. Ramirez-Torres, V. Turkowski, and T. S. Rahman, Time-dependent density-matrix functional theory for trion excitations: Application to monolayer ${\mathrm{mos}}_{2}$ and other transition-metal dichalcogenides, Phys. Rev. B 90, 085419 (2014).
  * Kidd  _et al._ [2016] D. W. Kidd, D. K. Zhang, and K. Varga, Binding energies and structures of two-dimensional excitonic complexes in transition metal dichalcogenides, Phys. Rev. B 93, 125423 (2016).
  * Mayers  _et al._ [2015] M. Z. Mayers, T. C. Berkelbach, M. S. Hybertsen, and D. R. Reichman, Binding energies and spatial structures of small carrier complexes in monolayer transition-metal dichalcogenides via diffusion monte carlo, Phys. Rev. B 92, 161404 (2015).
  * Szyniszewski  _et al._ [2017] M. Szyniszewski, E. Mostaani, N. D. Drummond, and V. Fal’Ko, Binding energies of trions and biexcitons in two-dimensional semiconductors from diffusion quantum monte carlo calculations, Physical Review B 95, 081301 (2017).
  * Van der Donck  _et al._ [2017] M. Van der Donck, M. Zarenia, and F. M. Peeters, Excitons and trions in monolayer transition metal dichalcogenides: A comparative study between the multiband model and the quadratic single-band model, Phys. Rev. B 96, 035131 (2017).
  * Drüppel  _et al._ [2017] M. Drüppel, T. D. Peter, Krüger, and M. Rohlfing, Diversity of trion states and substrate effects in the optical properties of an MoS2 monolayer, Nature Communications 8, 2117 (2017).
  * Kezerashvili and Tsiklauri [2017] R. Y. Kezerashvili and S. M. Tsiklauri, Trion and biexciton in monolayer transition metal dichalcogenides, Few-Body Systems 58, 1 (2017).
  * Filikhin  _et al._ [2018a] I. Filikhin, R. Y. Kezerashvili, and B. Vlahovic, On binding energy of trions in bulk materials, Physics Letters A 382, 787 (2018a).
  * Filikhin  _et al._ [2018b] I. Filikhin, R. Y. Kezerashvili, S. M. Tsiklauri, and B. Vlahovic, Trions in bulk and monolayer materials: Faddeev equations and hyperspherical harmonics, Nanotechnology 29, 124002 (2018b).
  * Filikhin  _et al._ [2018c] I. Filikhin, R. Y. Kezerashvili, and B. Vlahovic, On binding energy of trions in bulk materials, Physics Letters A 382, 787 (2018c).
  * Chang and Chang [2021] Y.-W. Chang and Y.-C. Chang, Variationally optimized orbital approach to trions in two-dimensional materials, The Journal of Chemical Physics 155, 024110 (2021).
  * Cavalcante  _et al._ [2018] L. S. R. Cavalcante, D. R. da Costa, G. A. Farias, D. R. Reichman, and A. Chaves, Stark shift of excitons and trions in two-dimensional materials, Physical Review B 98, 245309 (2018).
  * Yamaguchi [1954] Y. Yamaguchi, Two-nucleon problem when the potential is nonlocal but separable. i, Phys. Rev. 95, 1628 (1954).
  * Glöckle [2012] W. Glöckle,  _The quantum mechanical few-body problem_ (Springer Science & Business Media, 2012).
  * Blume [2012] D. Blume, Few-body physics with ultracold atomic and molecular systems in traps, Reports on Progress in Physics 75, 046401 (2012).
  * Nishida  _et al._ [2013] Y. Nishida, Y. Kato, and C. D. Batista, Efimov effect in quantum magnets, Nature Physics 9, 93 (2013).
  * Adhikari  _et al._ [1988] S. K. Adhikari, A. Delfino, T. Frederico, I. Goldman, and L. Tomio, Efimov and thomas effects and the model dependence of three-particle observables in two and three dimensions, Physical Review A 37, 3666 (1988).
  * Faddeev [1960] L. Faddeev, Scattering theory for a three-particle system, Zhur. Eksptl’. i Teoret. Fiz. 39 (1960).
  * Rytova [1967] N. S. Rytova, The screened potential of a point charge in a thin film, Moscow University Physics Bulletin 3, 18 (1967).
  * Klots  _et al._ [2014] A. Klots, A. Newaz, B. Wang, D. Prasai, H. Krzyzanowska, J. Lin, D. Caudel, N. Ghimire, J. Yan, B. Ivanov, _et al._ , Probing excitonic states in suspended two-dimensional semiconductors by photocurrent spectroscopy, Scientific reports 4, 1 (2014).
  * Zhang  _et al._ [2016] C. Zhang, C. Gong, Y. Nie, K.-A. Min, C. Liang, Y. J. Oh, H. Zhang, W. Wang, S. Hong, L. Colombo, R. M. Wallace, and K. Cho, Systematic study of electronic structure and band alignment of monolayer transition metal dichalcogenides in van der waals heterostructures, 2D Materials 4, 015026 (2016).
  * 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).
  * Deltuva  _et al._ [2005] A. Deltuva, A. Fonseca, and P. Sauer, Momentum-space treatment of the coulomb interaction in three-nucleon reactions with two protons, Physical Review C 71, 054005 (2005).
  * Wang  _et al._ [2018] G. Wang, A. Chernikov, M. M. Glazov, T. F. Heinz, X. Marie, T. Amand, and B. Urbaszek, Colloquium: Excitons in atomically thin transition metal dichalcogenides, Rev. Mod. Phys. 90, 021001 (2018).
  * Schmidt  _et al._ [2003] W. G. Schmidt, S. Glutsch, P. H. Hahn, and F. Bechstedt, Efficient $\mathcal{O}{(N}^{2})$ method to solve the bethe-salpeter equation, Phys. Rev. B 67, 085307 (2003).
  * Mohseni  _et al._ [2021] K. Mohseni, A. J. Chaves, D. R. da Costa, T. Frederico, and M. R. Hadizadeh, Three-boson stability for boosted interactions towards the zero-range limit, Physics Letters B 823, 136773 (2021).
  * Hadizadeh  _et al._ [2020] M. Hadizadeh, M. Radin, and K. Mohseni, A three-dimensional momentum-space calculation of three-body bound state in a relativistic faddeev scheme, Scientific Reports 10, 1 (2020).
  * Ahmadi  _et al._ [2020] M. Ahmadi, M. R. Hadizadeh, M. Radin, and S. Bayegan, Novel regularization scheme for nucleon-nucleon lattice simulations with effective field theory, Phys. Rev. C 102, 044001 (2020).
  * Hadizadeh  _et al._ [2014] M. R. Hadizadeh, C. Elster, and W. N. Polyzou, Relativistic three-body bound state in a 3d formulation, Phys. Rev. C 90, 054002 (2014).
  * Hadizadeh  _et al._ [2011] M. R. Hadizadeh, M. T. Yamashita, L. Tomio, A. Delfino, and T. Frederico, Scaling properties of universal tetramers, Phys. Rev. Lett. 107, 135304 (2011).
  * Bayegan  _et al._ [2008] S. Bayegan, M. R. Hadizadeh, and M. Harzchi, Three-nucleon bound state in a spin-isospin dependent three dimensional approach, Phys. Rev. C 77, 064005 (2008).
  * Hadizadeh and Bayegan [2007] M. Hadizadeh and S. Bayegan, Four-body bound-state calculations in three-dimensional approach, Few-Body Systems 40, 171 (2007).
  * Hadizadeh  _et al._ [2012] M. R. Hadizadeh, M. T. Yamashita, L. Tomio, A. Delfino, and T. Frederico, Binding and structure of tetramers in the scaling limit, Phys. Rev. A 85, 023610 (2012).
  * Hüber  _et al._ [1997] D. Hüber, H. Witała, A. Nogga, W. Glöckle, and H. Kamada, A new look into the partial-wave decomposition of three-nucleon forces, Few-Body Systems 22, 107 (1997).