# A trion in a magnetic field revisited

I. A. Aleksandrov  Department of Physics, Saint Petersburg State University,
Universitetskaya Naberezhnaya 7/9, Saint Petersburg 199034, Russia  Ioffe
Institute, Politekhnicheskaya street 26, Saint Petersburg 194021, Russia
<EMAIL_ADDRESS>A. Kudlis  Abrikosov Center for Theoretical Physics,
MIPT, Dolgoprudnyi, Moscow Region 141701, Russia<EMAIL_ADDRESS>I. A.
Shelykh  Science Institute, University of Iceland, Dunhagi-3, IS-107
Reykjavik, Iceland  Abrikosov Center for Theoretical Physics, MIPT,
Dolgoprudnyi, Moscow Region 141701, Russia<EMAIL_ADDRESS>
###### Abstract

We revisit the problem of a two dimensional trion in an external magnetic
field. We demonstrate that the approximations used previously for finding the
energy spectrum of this system break down in the experimentally accessible
range of magnetic fields. It is shown that the neglect of the Coulomb-induced
mixing of different Landau levels corresponding to non-interacting particles
leads to a strong underestimation of the trion binding energies even at
extremely high magnetic fields (hundreds of Tesla). Moreover, proper account
of the Coulomb effects for certain values of the parameters can lead to the
appearance of additional discrete trion states which were overlooked
previously. Finally, we provide a database of the matrix elements necessary
for calculation of the magnetotrion spectra for a wide class of materials.

##  I Introduction

Understanding of magneto-optical spectra of various quasi-two-dimensional (2D)
electron-hole ($e$–$h$) complexes in external magnetic fields is an important
problem in condensed matter theory Miura (2008); Arora (2021). The resonances
corresponding to such composite quasiparticles in magneto-optical spectra are
routinely observed in semiconductor quantum wells and 2D materials. These
experimental observations give strong motivation for theorists to provide a
corresponding quantitative description together with an unambiguous
classification of the spectra. The developed theoretical techniques should be
sufficiently flexible and allow one to examine a broad variety of materials
with different effective masses of carriers and types of the screening in
experimentally relevant diapasons of magnetic fields.

The fundamental basis for theoretical description of magneto-optic response
was laid by L. D. Landau decades ago in his seminal work introducing the
concept of Landau levels (LLs) for a single charged particle in an external
homogeneous magnetic field Landau (1930). The next step was consideration of
an exciton, a neutral composite particle consisting of an electron and a hole
bound together by the Coulomb attraction Akimoto and Hasegawa (1967);
Edelstein  _et al._ (1989); MacDonald and Ritchie (1986); Potemski  _et al._
(1991); Yang and Sham (1987); Bodnar  _et al._ (2017). In this case, it was
demonstrated that center-of-mass coordinates of an exciton can be separated
from the internal degrees of freedom Gor’kov and Dzyaloshinskii (1968); Avron
and Herbst (1978). This allowed one to build the theory of a magnetoexciton
using a perturbative approach Lerner and Lozovik (1980), WKB approximation
Akimoto and Hasegawa (1967), and brute force numerical modelling Edelstein
_et al._ (1989), and provide the classification of the corresponding spectra
MacDonald and Ritchie (1986).

Besides excitons, electron-hole complexes containing more than two particles,
such as biexcitons He  _et al._ (2016); Stevens  _et al._ (2018); Torche and
Bester (2021) and trions Finkelstein and Bar-Joseph (1995); Shields  _et al._
(1995); Sanvitto  _et al._ (2002); Smoleński  _et al._ (2019); Munteanu  _et
al._ (2000); Andronikov  _et al._ (2005); Klein  _et al._ (2022); Van Tuan
_et al._ (2022); Wang  _et al._ (2020); Li  _et al._ (2020) can be formed. The
latter case is of particular interest as a trion (T), being a composite
particle consisting of two holes and one electron ($\textup{X}^{+}$) or two
electrons and one hole ($\textup{X}^{-}$), possesses a net electric charge,
and its center-of-mass dynamics should be thus strongly affected by an
external magnetic field Stébé  _et al._ (1987); Chapman  _et al._ (1998);
Finkelstein  _et al._ (1995); Wojs and Hawrylak (1995); Palacios  _et al._
(1996); Whittaker and Shields (1997); Stébé  _et al._ (1998); Kozhevnikov  _et
al._ (1997); Hayne  _et al._ (1999); Riva  _et al._ (2002, 2001); Kuznetsov
_et al._ (2018); Filikhin  _et al._ (2018); Chang and Chang (2022, 2021);
Stébé and Moradi (2000); Kaur  _et al._ (2019). Although for the case of a
magnetotrion the center-of-mass coordinates cannot be straightforwardly
separated as in the case of an exciton, one can take advantage of an another
exact symmetry, namely magnetic translations, to reduce the number of the
degrees of freedom Schmelcher and Cederbaum (1991). An additional quantum
number characterizing trion states can be thus introduced, together with the
total orbital angular momentum $M_{z}$ along the magnetic-field direction, and
the spins of electrons $S_{e}$ and holes $S_{h}$ Dzyubenko and Sivachenko
(2000).

A number of different techniques were employed to analyze magnetotrion spectra
quantitatively, such as variational methods with Slater-type orbitals as trial
functions Chang and Chang (2021, 2022) and a stochastic variational approach
Riva  _et al._ (2002, 2001). These methods, however, possess an important
drawback as they do not present a clear classification of the energy levels.
Another possible approach is based on the perturbative treatment where the
mixing of the LLs due to the the Coulomb interaction is supposed to be small
Dzyubenko and Sivachenko (2000); Dzyubenko (2000, 2001); Kuznetsov  _et al._
(2018). However, as we will show below, this approach neglects an entire class
of the important matrix elements, and not only gives wrong values of trion
energies, but even fails to reproduce the correct number of the discrete
states in the experimentally relevant range of magnetic fields. Revisiting the
problem is thus a relevant task.

When considering a trion in an external magnetic field, one should analyze the
regimes of strong and weak fields separately. Indeed, on the one hand, in the
case of a weak magnetic field, the simplest natural approximation would be to
consider a trion as a single particle with the net charge $Q_{\textup{T}}=e$
and a mass equal to the total mass of the particles forming it,
$m_{\textup{T}}=2m_{e}+m_{h}$. In this limit, the Coulomb interaction binding
the particles together dominates over the magnetic field effects, so that

$E_{\textup{T}}\gg\hbar\omega_{\textup{T}}=\frac{\hbar|e|B}{m_{\textup{T}}c},$
(1)

where $E_{\textup{T}}$ is the binding energy of a trion and
$\omega_{\textup{T}}=|e|B/(m_{\textup{T}}c)$ is the trion cyclotron frequency.

On the other hand, in the case of strong magnetic fields, one can consider
three particles interacting with a magnetic field as almost independent ones.
It was argued that the Coulomb interaction in this case can be treated as a
perturbation leading to a weak mixing of the LLs corresponding to the
individual particles. This approach is expected to be valid if Dzyubenko and
Sivachenko (2000); Dzyubenko (2000, 2001)

$\hbar\omega_{e},\,\hbar\omega_{h},\,|\hbar\omega_{e}-\hbar\omega_{h}|\gg\frac{e^{2}}{\varepsilon\lambda}=\frac{e^{2}}{\varepsilon}\sqrt{\frac{|e|B}{\hbar
c}},$ (2)

where $\omega_{e}$ and $\omega_{h}$ are the cyclotron frequencies of an
electron and a hole, respectively, $\varepsilon$ is a permittivity of a
medium, and $\lambda=\sqrt{\hbar c/(|e|B)}$, with is a magnetic length. This
corresponds to fields

$B\gg\frac{|e|^{3}m_{e,h}^{2}c}{\varepsilon^{2}\hbar^{2}}.$ (3)

Here we employ the Gaussian-cgs units (the SI expressions can be obtained by
substituting $\varepsilon\to 4\pi\varepsilon\varepsilon_{0}$ and omitting
$c$). The condition (2) seems to have a clear physical meaning, being
equivalent to the requirement that the characteristic mixing energy between
the LLs were much smaller than the distance between them.

However, one can note that there are certain problems with this criterion.
First, Eq. (3) imposes a very strong restriction on the values of the magnetic
field strength. For instance, in the case of GaAs, it corresponds to $B\sim
1000~{}\textup{T}$ or stronger, which is out of the experimentally accessible
range.

Second, even if the condition (3) is satisfied, one can argue that the mixing
of the LLs should be taken into account anyway. Indeed, from the viewpoint of
the perturbation theory, the corresponding second-order energy correction is
proportional to the square of the Coulomb matrix element (which scales as
$\sqrt{B}$), divided by the distance between the LLs (which is proportional to
$B$). This combination is thus independent on the magnetic field, which means
that the increase of the latter does not reduce the role of the mixing terms.

These two considerations motivated us to go beyond the approximation used in
Refs. Dzyubenko and Sivachenko (2000); Dzyubenko (2000, 2001) and treat the
Coulomb interaction terms in the magnetotrion problem exactly. In this paper,
we present the corresponding development of the operator method introduced in
Refs. Dzyubenko and Sivachenko (2000); Dzyubenko (2000, 2001) generalizing it
for the case when mixing of the LLs is accounted for, and provide a reader
with all necessary formulae for all relevant matrix elements. We show that the
approach used previously considerably underestimates the trion binding energy
and, moreover, fails to reproduce completely the discrete part of the
spectrum.

The paper has the following structure. In Sec. II, we briefly outline our
theoretical approach (detailed derivations are given in the Appendices). In
Sec. III, we provide a quantitative analysis of trion energy spectra in the
cases of GaAs and CdTe quantum wells (QWs) for various values of the magnetic
field strength. Conclusions are summarized in Sec. IV.

Throughout the text, we denote the electron (hole) charge and mass by $e<0$
and $m_{e}$ ($-e$ and $m_{h})$.

##  II Theoretical description

In this section, we discuss the main theoretical aspects of the calculations
of trion energy spectra. We briefly describe the symmetries of the problem,
the basis used, and present the final expressions for the matrix elements of
the full trion Hamiltonian. In the main text we present only the most
important steps of the derivation, whereas detailed description of the
calculation procedures is given in the Appendices. We tried to make the paper
self-consistent, so Appendices contain all relevant expressions for the cases
of individual electrons, excitons, trions, and unbound electron-exciton pairs
placed in external magnetic fields (see Appendix A, B, C, and D,
respectively).

The Hamiltonian of a three-particle system involving a pair of two dimensional
(2D) electrons with coordinates $\mathbf{r}_{1}$ and $\mathbf{r}_{2}$ and a
hole with coordinates $\mathbf{r}_{h}$ reads

$\hat{H}^{(\textup{T})}=\hat{H}_{0}+\hat{V}_{eh}+\hat{V}_{ee},$ (4)

where

$\displaystyle\hat{H}_{0}$
$\displaystyle=\sum_{j=1,2}\frac{\hat{\boldsymbol{\pi}}_{j}^{2}}{2m_{e}}+\frac{\hat{\boldsymbol{\pi}}_{h}^{2}}{2m_{h}},$
(5) $\displaystyle\hat{V}_{eh}$
$\displaystyle=-\sum_{j=1,2}\frac{e^{2}}{\lambda}V(|\mathbf{r}_{j}-\mathbf{r}_{h}|),$
(6) $\displaystyle\hat{V}_{ee}$
$\displaystyle=\frac{e^{2}}{\lambda}V(|\mathbf{r}_{1}-\mathbf{r}_{2}|).$ (7)

Here the momenta operators are obtained by the canonical substitution
$\hat{\boldsymbol{\pi}}_{j}=-i\hbar\boldsymbol{\nabla}_{j}-(e/c)\boldsymbol{A}(\boldsymbol{r}_{j})$,
$\hat{\boldsymbol{\pi}}_{h}=-i\hbar\boldsymbol{\nabla}_{h}+(e/c)\mathbf{A}(\mathbf{r}_{h})$,
$\lambda=\sqrt{\hbar c/(|e|B)}$, and $V(r)$ is a dimensionless function
describing the interaction potential. In the current paper, we focus on the
case of conventional semiconductor QWs and thus employ the standard expression
for a Coulomb potential,

$V(r)=\frac{\lambda}{\varepsilon r}.$ (8)

We did not wright explicitly the terms, correspopnding to the direct action of
the magnetic field on spins of individual electrons and holes, as their
account is trivial.

In contrast to the case of a single charged particle (an electron or a hole),
where the spectrum is given by the discrete harmonic oscillator series (LLs),
the spectrum of a multi-particle complex is substantially modified by the
Coulomb interaction, and the problem cannot be solved analytically. Therefore,
we first construct the basis consisting of the eigenfunctions of the free
particle Hamiltonian $\hat{H}_{0}$ and then diagonalize the matrix of the
total Hamiltonian $\hat{H}^{(\textup{T})}$ numerically. When choosing the
basis, we fully take into account the complete symmetries of the system, which
allows us to substantially reduce the calculation costs.

In the case of the three interacting particles, the operator of the magnetic
translations defined as Lerner and Lozovik (1980); Dzyubenko and Sivachenko
(2000); Dzyubenko (2000, 2001)

$\hat{\mathbf{K}}=\sum_{j=1,2}\hat{\boldsymbol{\pi}}_{j}+\hat{\boldsymbol{\pi}}_{h}+\frac{e}{c}\sum_{j=1,2}\mathbf{B}\times\mathbf{r}_{j}-\frac{e}{c}\,\mathbf{B}\times\mathbf{r}_{h}$
(9)

commutes with the Hamiltonian $\hat{H}^{(\textup{T})}$. Note that this
commutator is proportional to the net electric charge of the system, so the
components of the corresponding operator for an exciton do commute (see
Appendix B). However, the $x$ and $y$ components of $\hat{\mathbf{K}}$ do not
commute with each other, $[\hat{K}_{x},\,\hat{K}_{y}]=i|e|\hbar B/c$, which
means that in the case of a trion, when constructing the basis, one can take
advantage of only one integral of motion (for instance, either $\hat{K}_{x}$
or $\hat{K}_{y}$). Instead of choosing a specific component of the vector
$\hat{\mathbf{K}}$ (either $\hat{K}_{x}$ or $\hat{K}_{y}$), one can use the
operator $\hat{\mathbf{K}}^{2}$ Dzyubenko (2000), which commutes with the
Hamiltonian as well.

Let us perform the following coordinate transformation:

$\displaystyle\mathbf{R}$
$\displaystyle=\frac{1}{\sqrt{2}}\,\big{(}\mathbf{r}_{1}+\mathbf{r}_{2}\big{)},\quad\mathbf{r}=\frac{1}{\sqrt{2}}\,\big{(}\mathbf{r}_{1}-\mathbf{r}_{2}\big{)}.$
(10)

In the most straightforward approach the basis functions corresponding to the
case of non-interacting particles can be chosen as:

$\varphi^{(\textup{T})}_{n_{1}m_{1}n_{2}m_{2}n_{h}m_{h}}(\mathbf{r},\mathbf{R},\mathbf{r}_{h})\\\
=\varphi^{(e)}_{n_{1}m_{1}}(\mathbf{r})\varphi^{(e)}_{n_{2}m_{2}}(\mathbf{R})\varphi^{(h)}_{n_{h}m_{h}}(\mathbf{r}_{h}),$
(11)

where $\varphi^{(h)}_{nm}(\mathbf{r})=[\varphi^{(e)}_{nm}(\mathbf{r})]^{*}$
correspond to single particle wavefunctions in a magnetic field (see Appendix
A). The quantum number $n=0,1,2,\dots$ is the LL number, and $m=0,1,2,\dots$
defines the angular momentum projection. However, this basis set completely
disregards the translational symmetry of the system as the corresponding
functions are characterized only by a well-defined projection of the total
angular momentum $M_{z}=(n_{1}-m_{1})+(n_{2}-m_{2})-(n_{h}-m_{h})$ and total
electron spin $S_{e}$ [$S_{e}=0$ (spin-singlet state) for even $n_{1}-m_{1}$
and $S_{e}=1$ (spin-triplet state) for odd $n_{1}-m_{1}$; see Appendix C].

To account for the translational symmetry, one can perform a unitary
transformation making the basis functions (11) also eigenfunctions of
$\hat{\mathbf{K}}^{2}$ Dzyubenko (2000) (see Appendix C for details), and get
the following set of the basis functions

$\displaystyle\psi_{n_{1}n_{2}n_{h}ml}(\mathbf{r},\mathbf{R},\mathbf{r}_{h})$
$\displaystyle=\frac{1}{\sqrt{N_{n_{1}n_{2}n_{h}ml}}}\left[\dfrac{1}{\sqrt{2}}\left(\xi-\partialderivative{\xi^{*}}\right)\right]^{n_{1}}\left[\dfrac{1}{\sqrt{2}}\left(\xi_{R}-\partialderivative{\xi^{*}_{R}}\right)\right]^{n_{2}}\left[\dfrac{1}{\sqrt{2}}\left(\xi^{*}_{h}-\partialderivative{\xi_{h}}\right)\right]^{n_{h}}\psi_{ml}(\mathbf{r},\mathbf{R},\mathbf{r}_{h}),$
(12)

where

$\displaystyle\psi_{ml}(\mathbf{r},\mathbf{R},\mathbf{r}_{h})=\left(\xi^{*}\right)^{m}\left(\xi_{h}\right)^{l}\tilde{\Phi}^{(\textup{T})}_{0}\\!(\mathbf{r},\\!\mathbf{R},\mathbf{r}_{h})$
(13)

and

$\displaystyle\tilde{\Phi}^{(\textup{T})}_{0}\\!(\mathbf{r},\\!\mathbf{R},\mathbf{r}_{h})\\!=\\!\dfrac{\exp{-\xi\xi^{*}\\!-\\!\xi_{R}\xi^{*}_{R}\\!-\\!\xi_{h}\xi^{*}_{h}\\!+\\!\sqrt{2}\xi_{h}\xi^{*}_{R}}}{\sqrt{2}\left(2\pi\right)^{3/2}\lambda^{3}}.$
(14)

is a normalized wave function of the ground state. The normalization factors
$N_{n_{1}n_{2}n_{h}ml}$ are calculated in Appendix C. The dimensionless
coordinates read: $\xi=z/(2\lambda)$, $\xi_{R}=z_{R}/(2\lambda)$, and
$\xi_{h}=z_{h}/(2\lambda)$, where $z=x+{i\mkern 1.0mu}y$. The angular momentum
projection of the state (12) reads

$M_{z}=n_{1}-m+n_{2}-(n_{h}-l),$ (15)

while the spin quantum number is $S_{e}=0$ ($S_{e}=1$) if $n_{1}-m$ is even
(odd). Note that the energies are degenerate with respect to the quantum
number $k$ associated with the operator $\mathbf{K}^{2}$, so we have set
$k=0$.

The matrix elements of the full Hamiltonian (4) are defined by the following
expression:

$\displaystyle
H^{(\textup{T})}_{n_{1}n_{2}n_{h}ml;n_{1}^{\prime}n_{2}^{\prime}n_{h}^{\prime}m^{\prime}l^{\prime}}$
$\displaystyle=E^{(\textup{T},0)}_{n_{1}+n_{2},n_{h}}\delta_{n_{1}n_{1}^{\prime}}\delta_{n_{2}n_{2}^{\prime}}\delta_{n_{h}n_{h}^{\prime}}\delta_{mm^{\prime}}\delta_{ll^{\prime}}$
$\displaystyle+V^{(eh)}_{n_{1}n_{2}n_{h}ml;n_{1}^{\prime}n_{2}^{\prime}n_{h}^{\prime}m^{\prime}l^{\prime}}$
$\displaystyle+V^{(ee)}_{n_{1}n_{2}n_{h}ml;n_{1}^{\prime}n_{2}^{\prime}n_{h}^{\prime}m^{\prime}l^{\prime}},$
(16)

where

$E^{(\textup{T},0)}_{n_{e},n_{h}}=\hbar\omega_{e}\bigg{(}1+n_{e}\bigg{)}+\hbar\omega_{h}\bigg{(}\frac{1}{2}+n_{h}\bigg{)}$
(17)

corresponds to the LLs of non-interacting particles, $n_{e}=n_{1}+n_{2}$ and
$\omega_{e,h}=|e|B/(m_{e,h}c)$. The Coulomb interaction matrix elements
$V^{(eh)}_{n_{1}n_{h}ml;n_{1}^{\prime}n_{2}^{\prime}n_{h}^{\prime}m^{\prime}l^{\prime}}$
and
$V^{(ee)}_{n_{1}n_{2}n_{h}ml;n_{1}^{\prime}n_{2}^{\prime}n_{h}^{\prime}m^{\prime}l^{\prime}}$
are bulky and given in Appendic C ( Eqs. (96) and (108), respectively).

In Ref. Dzyubenko (2000) the Coulomb mixing was taken into account only within
the subspaces corresponding to $n_{e}=n_{e}\textasciiacute$ and
$n_{h}=n_{h}\textasciiacute$, i.e. it was thus assumed, that $n_{e},n_{h}$
remain good quantum numbers of the full interaction problem. It was argued
that this approximation should be valid in the case of very strong magnetic
fields, where the magnetic energy dominates over the Coulomb energy. However,
as was discussed in the Introduction [see Eq. (3) and below], this approach
can hardly be justified in any physically relevant situation and exact
treatment of the Coulomb mixing is necessary.

The resulting spectrum contains discrete levels and continuous regions,
corresponding to the states of unbound electron-exciton pairs. To identify the
onset of the energy continuum, we additionally perform the independent
calculations for the case of an uncoupled electron and an exciton (see
Appendix D).

Figure 1: Lower part of the energy spectrum for a trion in a GaAs QW for
$B=30~{}\textup{T}$. The first vertical panel corresponds to the-single-
composite particle (SCP) approximation, where the energies are calculated by
means of Eq. (18). The second vertical panel displays the free LLs obtained
via Eq. (17). Near each horizontal line, we indicate the values of the pair
$(n_{e}=n_{1}+n_{2},n_{h})$. In the third panel, we present the results
obtained with use of the technique proposed in Ref. Dzyubenko and Sivachenko
(2000), where Coulomb mixing of different LLs was neglected. The fourth panel
displays the spectrum obtained in this work by diagonalization of the matrix
of the full Hamiltonian, (16), accounting for the LLs mixing. In the third and
fourth panels thin horizontal lines correspond to the discrete energy state of
magnetotrions, and color bars display energy continuums corresponding to
unbound electron-hole pairs. As one can see, even for very strong value of the
magnetic field the account of the Coulomb mixing is crucial. It affects
substantially both the boundary of the energy continuum and magnetotrion
binding energy, defined as the distance between the discrete trion state and
the boundary of energy continuum. Moreover, it also leads to the appearance of
a second discrete level. To make the figure less cumbersome, we have not
plotted all regions of the energy continuums in the third panel, which should
start down from every LL of non-interacting particles presented in the second
panel.

##  III Results and discussion

In this section, we examine the trion spectra for GaAs and CdTe QWs for
different values of the external magnetic field strengths. For these
conventional semiconductors, the interaction among charge carriers is well
described by the standard Coulomb potential (8).

For the case of GaAs, we use the following parameters: effective masses of
carriers $m_{e}=0.063m_{0}$, $m_{h}=0.51m_{0}$ ($m_{0}$ is the free electron
mass), static dielectric permittivity $\varepsilon=12.9$.

Our results for $B=30~{}\textup{T}$ are shown in Fig. 1. The corresponding
spectrum is presented in the fourth vertical panel from the left. Constructing
the matrix of the Hamiltonian, we accounted for the LLs with $n_{h}$ and
$n_{e}=n_{1}+n_{2}$ spanning from $0$ to $4$ (see the discussion of the
convergence of the procedure in Appendix G). We compared our results with
those obtained in Refs. Dzyubenko and Sivachenko (2000); Dzyubenko (2000,
2001), where the Coulomb mixing of the different LLs was neglected, see third
vertical panel from the left. The thin horizontal lines correspond to the
discrete energy state of magnetotrions, and the color bars display the energy
continua corresponding to unbound electron-hole pairs. For comparison we also
illustated the simple cases of a trion as a single composite particle (first
panel from the right)with the spectrum given by:

$\displaystyle
E^{\textup{SCP}}_{n}=\hbar\omega_{\textup{T}}\left(n+\dfrac{1}{2}\right),$
(18)

where the trion cyclotron frequency $\omega_{T}$ is defined in Eq. (1), and
LLs corresponding to three non-interacting particles with spectrum given by
Eq. (17) (second panel from the left). Both of these approximations yield
obviously only discrete energy levels.

Comparing third and fourth panels of Fig. 1 one can note that even for the
value of a magnetic field equal to several tens of Tesla, the account of the
Coulomb mixing is crucial. It substantially affects the boundary of the energy
continuum, which shifts down by $5.688~{}\textup{meV}$). The magnetotrion
binding energy, defined as the distance between the discrete trion state and
the boundary of the energy continuum, is increased from $1.298~{}\textup{meV}$
($M_{z}=-1$ and $S_{e}=1$) to $4.680~{}\textup{meV}$ ($M_{z}=0$ and
$S_{e}=0$). For details, see Appendices E and F.

Figure 2: Lower part of the energy spectrum for a trion in a GaAs QW for
$B=5~{}\textup{T}$. The first vertical panel corresponds to the-single-
composite particle (SCP) approximation, where the energies are calculated by
means of Eq. (18). The second vertical panel displays the free LLs obtained
via Eq. (17). Near each horizontal line, we indicate the values of the pair
$(n_{e}=n_{1}+n_{2},n_{h})$. In the third panel, we present the results
obtained with use of the technique proposed in Ref. Dzyubenko and Sivachenko
(2000), where Coulomb mixing of different LLs was neglected. The fourth panel
displays the spectrum obtained in this work by diagonalizing the matrix of the
full Hamiltonian, (16), accounting for the LLs mixing. In the third and fourth
panels thin horizontal lines correspond to the discrete energy state of
magnetotrions, and color bars display energy continuums corresponding to
unbound electron-hole pairs. One clearly sees the substantial shift of both
the energy continuum and the trion binding energy.

Moreover, exact treatment of the Coulomb mixing of LLs modifies the spectrum
quantitatively, leading to the appearence of an additonal discrete trion state
below the energetic continuum with binding energy of $1.043~{}\textup{meV}$.
This state corresponds to $M_{z}=-1$ and $S_{e}=1$.

We also considered the the case of $B=5~{}\textup{T}$, which is illustrated in
Fig. 2. Here the energy continuum moves down by $5.329~{}\textup{meV}$, and
the trion binding energy increases from $0.530~{}\textup{meV}$ to
$1.604~{}\textup{meV}$. Interestingly, for this value of the magnetic field,
an additional trion level does not appear.

For comparison, we also present the results for another widely used
semiconductor material, CdTe. It has effective masses $m_{e}=0.11m_{0}$,
$m_{h}=0.40m_{0}$ and relative permittivity $\varepsilon=11.0$. We show the
spectrum for $B=30~{}\textup{T}$ only (see Fig. 3). As in the case of GaAs,
our findings demonstrate that proper account for the Coulomb mixing of the LLs
is important and leads to the shift of the energy continuum by
$11.796~{}\textup{meV}$ and of the trion binding energy by
$2.735~{}\textup{meV}$.

Figure 3: Lower part of the energy spectrum for a trion in a CdTe QW for
$B=30~{}\textup{T}$. The first vertical panel corresponds to the-single-
composite particle (SCP) approximation, where the energies are calculated by
means of Eq. (18). The second vertical panel displays the free LLs obtained
via Eq. (17). Near each horizontal line, we indicate the values of the pair
$(n_{e}=n_{1}+n_{2},n_{h})$. In the third panel, we present the results
obtained by means of the technique proposed in Ref. Dzyubenko and Sivachenko
(2000), where Coulomb mixing of different LLs was neglected. The fourth
displays the spectrum obtained in this work by diagonalizing the matrix of the
full Hamiltonian, (16), accounting for the LLs mixing. In the third and fourth
panels thin horizontal lines correspond to the discrete energy state of
magnetotrions, and color bars display energy continuums corresponding to
unbound electron-hole pairs. One clearly sees the substantial shift of both
the energy continuum and the trion binding energy. Figure 4: The energy
spectrum of GaAs based system as function of the magnetic field $B$. Solid
green lines correspond to the two discrete trion states obtained in this work,
the dashed green line corresponds to the lower boundary of the energy
continuum, solid blue line corresponds to the trion bound state reported in
Ref. Dzyubenko (2000), dashed blue line indicates the corresponding energetic
continuum boundary. We also plot the position of the lowest LL corresponding
to non-interacting particles (solid black line). The increase of magnetic
field does not lead to the merging of green and blue lines, which means that
the role of the Coulomb mixing between different LLs does not decrease with
magnetic field. We accounted for LLs with $n_{h}$ and $n_{e}=n_{1}+n_{2}$
spanning from $0$ to $4$. In Appendix H, we also present the similar
dependence extracted on the basis of extrapolating procedure described there
(see Fig. 12).

In Fig. 4 we show the energy spectrum of the system as a function of the
magnetic field $B$. We display the positions of the two discrete trion states,
which were identified in our calculations (thin and thick solid green lines),
the lower boundary of the energy continuum (dashed green line) and compare
them to the results, presented in Ref. Dzyubenko (2000) shown by the solid
blue line (bound trion) and dashed line (continuum boundary). The position of
the lowest LL corresponding to non-interacting particles is represented by the
solid black line.

As one can clearly see, the increase of the magnetic field does not lead to
the merging of green and blue lines, which means that the role of the Coulomb
mixing between the different LLs does not decrease with increasing $B$. This
is in agreement with a qualitative argument that although the increase of the
magnetic field leads to LL separation, it also enhances the Coulomb matrix
elements due to the reduction of cyclotron radii ($\lambda\sim 1/\sqrt{B}$),
and these two effects compensate each other (see the discussion after Eq. (3)
in the Introduction).

We point out once more that in our calculations the spin structure of the wave
functions was not presented explicitly. The electron spin state $S_{e}$ only
defines the parity of the coordinate part of the wave function. The account of
the interaction of spin with magnetic field is trivial, and will lead to the
splitting of a state with $S_{e}$ into $(2S_{h}+1)(2S_{e}+1)=2(2S_{e}+1)$
sublevels.

##  IV Conclusion

In this study, we revisited the problem of a 2D trion in an external magnetic
field. We demonstrated, that the mixing between different Landau levels
induced by Coulomb interaction plays crucial role in the determination of the
energy spectrum of the problem and can not be neglected even in the limit of
very strong magnetic fields. This is because magnetic field increases not only
the distance between the Landau levels, but Coulomb matrix elements as well
via decrease of the magnetic length. We have shown that proper account of the
Coulomb mixing not only leads to the substantial increase of magnetotrion
binding energy, but can also lead to qualitative reshaping of the spectrum,
resulting in the appearance of additional discrete energy states.

## Acknowledgments

I.A.S. acknowledge financial support from Icelandic Research Fund (Rannis, the
project “Hybrid polaritonics”), “Priority 2030 Academic Leadership Program”
and “Goszadanie no. 2019-1246”.

##  Appendix A Single electron in a magnetic field

Here and in Appendices B and C, we first briefly outline how the energy
spectrum can be computed and then provide a detailed theoretical description
of the corresponding problems. In Appendix D, we give a direct recipe for
calculating the spectrum for an unbound electron-exciton pair, since we
presented all the necessary theoretical calculations in the previous sections.

###  A.1 Energy spectrum

In the case of a single electron (hole) in a magnetic field, the energy is
given by a simple formula describing the spectrum of a harmonic oscillator.
The corresponding wave functions $\varphi^{(e)}_{nm}(\mathbf{r})$ are
constructed explicitly and form a basis in which the Hamiltonian is diagonal:

$\displaystyle H^{(e)}_{n^{\prime};n}=E^{(e)}_{n}\delta_{nn^{\prime}},$ (19)

where $n,n^{\prime}=0,1,2,\dots$ and the eigenvalues correspond to the Landau
levels (LLs),

$E_{n}^{(e)}=\hbar\omega_{e}\bigg{(}\frac{1}{2}+n\bigg{)}.$ (20)

Here $\omega_{e}=|e|B/(m_{e}c)$ is the electron cyclotron frequency. The
energy is degenerate with respect to the second quantum number $m$. The $z$
projection of the orbital angular momentum is $m_{z}=n-m$.

###  A.2 Details

We start with the Hamiltonian, which in the case of a single electron
interacting with a constant uniform magnetic field reads

$\hat{H}^{(e)}=\frac{\hat{\boldsymbol{\pi}}_{e}^{2}}{2m_{e}},$ (21)

where
$\hat{\boldsymbol{\pi}}_{e}=-i\hbar\boldsymbol{\nabla}-(e/c)\mathbf{A}(\mathbf{r})$
and we employ a symmetric gauge,
$\mathbf{A}(\mathbf{r})=(\mathbf{B}\times\mathbf{r})/2$ with $\mathbf{B}$
directed along the $z$ axis. We assume that the particle is confined within
the $xy$ plane and its quantum state is described by a wave function depending
solely on $x$ and $y$ (or polar coordinates $r$ and $\varphi$).

In the external field specified above, the Hamiltonian (21) obviously commutes
with the operator of the angular momentum projection
$\hat{\ell}_{z}=-i\hbar\partial_{\varphi}$, so the eigenfunctions of $H^{(e)}$
will also be characterized by the corresponding quantum number $m_{z}$, i.e.,
one can utilize the ansatz $\psi(r,\varphi)=\mathrm{e}^{im_{z}\varphi}R(r)$.
The ordinary differential equation for $R(r)$ can then be solved analytically
(see, e.g., Ref. Landau and Lifshitz (1977)). Here we will use another
approach based on the formalism of creation and annihilation operators (see,
e.g., Ref. Dzyubenko (2000)).

First, let us introduce the one-particle normalized vacuum state in the
coordinate representation:

$\Phi^{(e)}_{0}(\mathbf{r})=\langle\mathbf{r}|0\rangle=\frac{1}{\lambda\sqrt{2\pi}}\,\mathrm{e}^{-r^{2}/(4\lambda^{2})},$
(22)

where

$\lambda=\sqrt{\frac{\hbar c}{|e|B}}$ (23)

is the so-called magnetic length, which was also introduced in the main text.
It is convenient to use a complex variable $z=r\mathrm{e}^{i\varphi}=x+iy$ and
its complex conjugate $z^{*}$. We will also need the following creation
operators in the coordinate representation:

$\displaystyle\hat{A}_{e}^{\dagger}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}\bigg{(}\frac{z}{2\lambda}-2\lambda\,\frac{\partial}{\partial
z^{*}}\bigg{)},$ (24) $\displaystyle\hat{B}_{e}^{\dagger}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}\bigg{(}\frac{z^{*}}{2\lambda}-2\lambda\,\frac{\partial}{\partial
z}\bigg{)},$ (25)

Together with the annihilation operators $\hat{A}_{e}$ and $\hat{B}_{e}$, they
obey the relations $[\hat{A}_{e},\,\hat{A}_{e}^{\dagger}]=1$,
$[\hat{B}_{e},\,\hat{B}_{e}^{\dagger}]=1$, and
$[\hat{A}_{e},\,\hat{B}_{e}]=[\hat{A}_{e},\,\hat{B}_{e}^{\dagger}]=0$. These
operators allow one to generate the following wave functions:

$\varphi^{(e)}_{nm}(\mathbf{r})=\frac{1}{\sqrt{n!m!}}\,(\hat{A}_{e}^{\dagger})^{n}(\hat{B}_{e}^{\dagger})^{m}\Phi^{(e)}_{0}(\mathbf{r}),$
(26)

where $n$, $m=0$, $1$, $2$, … One can explicitly verify that these functions
form an orthonormal and complete set of the Hamiltonian eigenfunctions and the
corresponding energies are given by Eq. (20). Instead of the non-negative
integers $n$ and $m$, one can also label the solutions by the following
numbers:

$\displaystyle\tilde{n}$ $\displaystyle=$
$\displaystyle\mathrm{min}\,(n,\,m),$ (27) $\displaystyle\tilde{m}$
$\displaystyle=$ $\displaystyle m-n.$ (28)

These allow one to provide the following explicit form of the functions (26):

$\displaystyle\varphi^{(e)}_{nm}(\mathbf{r})$ $\displaystyle=$
$\displaystyle\frac{(-1)^{\tilde{n}}}{\lambda}\sqrt{\frac{2^{|\tilde{m}|}\tilde{n}!}{2\pi(\tilde{n}+|\tilde{m}|)!}}\bigg{(}\frac{r}{2\lambda}\bigg{)}^{|\tilde{m}|}$
(29) $\displaystyle\times$
$\displaystyle\mathrm{e}^{-i\tilde{m}\varphi}L^{|\tilde{m}|}_{\tilde{n}}\bigg{(}\frac{r^{2}}{2\lambda^{2}}\bigg{)}\,\mathrm{e}^{-r^{2}/(4\lambda^{2})},$

where $L^{|\tilde{m}|}_{\tilde{n}}$ are the Laguerre polynomials. In terms of
$n$ and $m$, the explicit form (29) strongly depends on whether $n>m$ or not,
so using $\tilde{n}$ and $\tilde{m}$ is often more convenient. The angular
momentum projection reads $m_{z}=-\tilde{m}=n-m$.

##  Appendix B Exciton $(\textup{X})$ in a magnetic field

In the case of two- and theree-particle systems considered in what follows,
the free LLs are modified due to the Coulomb interaction and the Hamiltonian
eigenfunctions cannot be found analytically. Next, we will discuss how the
energy eigenvalues can be obtained by diagonalizing the corresponding matrix.

###  B.1 Spectrum calculation recipe

The exciton energy spectrum turns out to be striped — it contains continuous
regions separated by empty areas. To evaluate the energies, one has first to
separate the center-of-mass motion by introducing a conserved vector
$\mathbf{K}$ (see the next section for details). Using then the basis set
(39), one obtains the matrix elements in the form

$H^{(\textup{X})}_{\mathbf{K},n^{\prime}m^{\prime};nm}=E^{(\textup{X},0)}_{nm}\delta_{nn^{\prime}}\delta_{mm^{\prime}}+V_{\mathbf{K},n^{\prime}m^{\prime};nm},$
(30)

where $n$ and $m$ are non-negative integer numbers, $E^{(\textup{X},0)}_{nm}$
is defined in Eq. (40), while the expression for
$V_{\mathbf{K},n^{\prime}m^{\prime};nm}$ is presented in Eq. (52). The
Hamiltonian is block-diagonal in $M_{z}=n-m$, due to the angular momentum
conservation. Having constructed the matrix (30), we should diagonalize it
varying $\mathbf{K}$ and $M_{z}$. In fact, the matrix elements depend only on
the magnitude $|\mathbf{K}|$ and $M_{z}$. Within numerical computations, one
has to make sure that for given $|\mathbf{K}|$ and $M_{z}$, the quantum
numbers $n$ and $m$ change in a sufficiently large range, so that the results
converge to a desired accuracy.

###  B.2 Details

As will be seen below, in the case of a two-dimensional electron-hole
($e$–$h$) system, it is possible to separate the center-of-mass coordinates
and make use of the one-particle solutions (29). The two-particle Hamiltonian
has the form

$\hat{H}^{(\textup{X})}=\frac{\hat{\boldsymbol{\pi}}_{e}^{2}}{2m_{e}}+\frac{\hat{\boldsymbol{\pi}}_{h}^{2}}{2m_{h}}-\frac{e^{2}}{\lambda}V(|\mathbf{r}_{e}-\mathbf{r}_{h}|),$
(31)

where
$\hat{\boldsymbol{\pi}}_{e}=-i\hbar\boldsymbol{\nabla}_{e}-(e/c)\mathbf{A}(\mathbf{r}_{e})$,
$\hat{\boldsymbol{\pi}}_{h}=-i\hbar\boldsymbol{\nabla}_{h}+(e/c)\mathbf{A}(\mathbf{r}_{h})$,
and $V(r)$ is a dimensionless function governing the interaction potential
[for the usual Coulomb interaction without the screening effects, it reads
$V(r)=\lambda/(\varepsilon r)$]. Although the interparticle-interaction law
has nothing to do with the magnetic length, its typical values does involve
$\lambda$ due to the localization of the wave functions in the external
magnetic field [for instance, the matrix elements of
$V(r)=\lambda/(\varepsilon r)$ computed in the basis of the Landau wave
functions are $\lambda$-independent].

Let us now take advantage of the translational symmetry of the problem and
separate the center-of-mass motion. As was shown in Ref. Gor’kov and
Dzyaloshinskii (1968), the Hamiltonian (31) commutes with the following vector
operator:

$\displaystyle\hat{\mathbf{K}}$ $\displaystyle=$
$\displaystyle\hat{\boldsymbol{\pi}}_{e}+\hat{\boldsymbol{\pi}}_{h}+\frac{e}{c}\,\mathbf{B}\times(\mathbf{r}_{e}-\mathbf{r}_{h})$
(32) $\displaystyle=$
$\displaystyle-i\hbar\boldsymbol{\nabla}_{e}-i\hbar\boldsymbol{\nabla}_{h}+\frac{e}{c}\,\big{[}\mathbf{A}(\mathbf{r}_{e})-\mathbf{A}(\mathbf{r}_{h})\big{]}.$

Since $\hat{K}_{x}$ and $\hat{K}_{y}$ also commute with each other, the
Hamiltonian eigenfunctions can be characterized by the corresponding vectors
$\mathbf{K}$ consisting of two quantum numbers $K_{x}$ and $K_{y}$. Note that
in the absence of the external field, $\mathbf{K}$ is the total momentum of
the particles. Let us introduce the center-of-mass position $\mathbf{R}$ and
relative position $\mathbf{r}$:

$\displaystyle\mathbf{R}$ $\displaystyle=$
$\displaystyle\frac{m_{e}\mathbf{r}_{e}+m_{h}\mathbf{r}_{h}}{m_{e}+m_{h}},$
(33) $\displaystyle\mathbf{r}$ $\displaystyle=$
$\displaystyle\mathbf{r}_{e}-\mathbf{r}_{h}.$ (34)

The operator (32) takes the form

$\hat{\mathbf{K}}=-i\hbar\boldsymbol{\nabla}_{R}+\frac{e}{2c}\,\mathbf{B}\times\mathbf{r}.$
(35)

The Hamiltonian (31) now reads

$\displaystyle\hat{H}^{(\textup{X})}$ $\displaystyle=$
$\displaystyle-\frac{\hbar^{2}}{2m_{eh}}\,\Delta_{R}-\frac{\hbar^{2}}{2\mu}\,\Delta_{r}$
(36) $\displaystyle+$
$\displaystyle\frac{ie\hbar}{2c}\Bigg{\\{}\frac{1}{m_{eh}}\big{(}\mathbf{B}\times\mathbf{r}\big{)}\boldsymbol{\nabla}_{R}+\frac{1}{\mu}\Big{[}\mathbf{B}\times\big{(}\mathbf{R}+\gamma\mathbf{r}\big{)}\Big{]}\boldsymbol{\nabla}_{r}\Bigg{\\}}$
$\displaystyle+$ $\displaystyle\frac{e^{2}B^{2}}{8\mu
c^{2}}\bigg{(}R^{2}+\beta
r^{2}+2\gamma\mathbf{r}\mathbf{R}\bigg{)}-\frac{e^{2}}{\lambda}V(r),$

where $m_{eh}=m_{e}+m_{h}$, $\mu=m_{e}m_{h}/m_{eh}$,
$\gamma=(m_{h}-m_{e})/m_{eh}$, and $\beta=(m_{e}^{3}+m_{h}^{3})/m_{eh}^{3}$.
The wave function can then be represented as Gor’kov and Dzyaloshinskii (1968)

$\varphi^{(\textup{X})}_{\mathbf{K}}(\mathbf{r},\mathbf{R})=\mathrm{exp}\bigg{[}\frac{i}{\hbar}\bigg{(}\mathbf{K}-\frac{e}{2c}\,\mathbf{B}\times\mathbf{r}\bigg{)}\mathbf{R}\bigg{]}\\\
{}\times\mathrm{exp}\bigg{(}\frac{i\gamma}{2\hbar}\,\mathbf{K}\mathbf{r}\bigg{)}\varphi^{(\textup{X})}(\mathbf{r}-\mathbf{r}_{0}).$
(37)

It is convenient to choose
$\mathbf{r}_{0}=-c(\mathbf{B}\times\mathbf{K})/(eB^{2})$ in order to obtain
the following equation for $\varphi^{(\textup{X})}(\mathbf{r})$:

$\bigg{[}-\frac{\hbar^{2}}{2\mu}\,\Delta+\frac{ie\hbar\gamma}{2\mu
c}\,(\mathbf{B}\times\mathbf{r})\boldsymbol{\nabla}+\frac{e^{2}B^{2}}{8\mu
c^{2}}\,r^{2}\\\
{}-\frac{e^{2}}{\lambda}V(|\mathbf{r}+\mathbf{r}_{0}|)-E\bigg{]}\varphi^{(\textup{X})}(\mathbf{r})=0.$
(38)

Note that the wave function $\varphi^{(\textup{X})}(\mathbf{r})$ does not
involve $\mathbf{K}$.

First, we will neglect the Coulomb interaction $V$ and find the corresponding
two-particle wave functions and energy eigenvalues (free LLs). The total
Hamiltonian $\hat{H}^{(\textup{X})}$ will be then considered in this basis,
and the corresponding matrix elements are given by Eq. (30). The solutions of
Eq. (38) without the $V$ term are exactly those given in Eq. (29), so the
normalized two-particle wave functions within the zeroth order with respect to
the Coulomb interaction have the form

$\varphi^{(\textup{X})}_{\mathbf{K},nm}(\mathbf{r},\mathbf{R})=\frac{1}{\sqrt{2\pi\hbar}}\mathrm{exp}\bigg{[}\frac{i}{\hbar}\bigg{(}\mathbf{K}-\frac{e}{2c}\,\mathbf{B}\times\mathbf{r}\bigg{)}\mathbf{R}\bigg{]}\\\
{}\times\mathrm{exp}\bigg{(}\frac{i\gamma}{2\hbar}\,\mathbf{K}\mathbf{r}\bigg{)}\varphi^{(e)}_{nm}(\mathbf{r}-\mathbf{r}_{0}).$
(39)

The corresponding energies involve the cyclotron frequency with the reduced
mass, $\omega_{0}=|e|B/(\mu c)$, and read

$\displaystyle E^{(\textup{X},0)}_{nm}$ $\displaystyle=$
$\displaystyle\hbar\omega_{0}\bigg{(}\frac{1}{2}+\frac{1+\gamma}{2}\,n+\frac{1-\gamma}{2}\,m\bigg{)}$
(40) $\displaystyle=$
$\displaystyle\hbar\omega_{0}\bigg{(}\frac{1}{2}+\tilde{n}+\frac{|\tilde{m}|-\gamma\tilde{m}}{2}\bigg{)}.$

Here the superscript “0” indicates that the corresponding energies are
obtained to zeroth order in the Coulomb interaction. Note that
$\omega_{0}=\omega_{e}+\omega_{h}$, where $\omega_{h}=|e|B/(m_{h}c)$. In the
heavy-hole limit $m_{h}=\infty$, one obtains $\gamma=1$, $\mu=m_{e}$,
$\omega_{h}=0$ and recovers the electron LLs (20). The angular momentum
projection of the two-particle system is again $M_{z}=-\tilde{m}=n-m$.

Next, one has to compute the matrix elements of the $V$ operator in Eq. (31)
using the basis functions (39). The integration over $\mathbf{R}$ is trivial
and the matrix is diagonal with respect to $\mathbf{K}$:

$V_{\mathbf{K}^{\prime},n^{\prime}m^{\prime};\mathbf{K},nm}\equiv-(e^{2}/\lambda)\big{\langle}\varphi^{(\textup{X})}_{\mathbf{K}^{\prime},n^{\prime}m^{\prime}}\big{|}V(r)\big{|}\varphi^{(\textup{X})}_{\mathbf{K},nm}\big{\rangle}\\\
=-(e^{2}/\lambda)\delta(\mathbf{K}-\mathbf{K}^{\prime})\int\differential\mathbf{r}\big{[}\varphi^{(e)}_{n^{\prime}m^{\prime}}(\mathbf{r})\big{]}^{*}V(|\mathbf{r}+\mathbf{r}_{0}|)\varphi^{(e)}_{nm}(\mathbf{r})\\\
\equiv\delta(\mathbf{K}-\mathbf{K}^{\prime})V_{\mathbf{K},n^{\prime}m^{\prime};nm}.$
(41)

There are several technical difficulties which one encounters in evaluating
the integral over $\mathbf{r}$. First, note that the function
$V(|\mathbf{r}+\mathbf{r}_{0}|)$ is symmetric about
$\mathbf{r}=-\mathbf{r}_{0}$, whereas the wave functions are symmetric about
the origin (apart from the angular part $\mathrm{e}^{-i\tilde{m}\varphi}$).
Second, the radial parts of the wave functions involve special functions
(Laguerre polynomials), which complicates the calculations. Moreover, the
explicit form of the electrostatic potential may also be highly complex.

In the present study, we compute integrals containing monomials instead of the
Laguerre polynomials and then combine the results taking into account the
corresponding coefficients. In order to do this, it is necessary to inspect
the structure of the wave functions. First, it is convenient to introduce non-
normalized wave functions, which will be denoted by $\psi$ (normalized
functions are denoted by $\varphi$). Second, we rewrite the creation operators
(24) in terms of a new dimensionless variable $\xi=z/(2\lambda)$:

$\displaystyle\hat{A}^{\dagger}_{e}=\frac{1}{\sqrt{2}}\left(\xi-\partialderivative{\xi^{*}}\right),\quad\hat{B}^{\dagger}_{e}=\frac{1}{\sqrt{2}}\left(\xi^{*}-\partialderivative{\xi}\right).$
(42)

The vacuum state (22) in terms of $\xi$ reads

$\displaystyle\Phi^{(e)}_{0}(\mathbf{r})=\frac{1}{\lambda\sqrt{2\pi}}\,\mathrm{e}^{-\xi\xi^{*}}.$
(43)

The arbitrary wave function can be obtained via

$\displaystyle\varphi^{(\textup{X})}_{\mathbf{K},nm}(\mathbf{r},\mathbf{R})=\phi^{(\textup{X})}_{\mathbf{K}}(\mathbf{r},\mathbf{R})\dfrac{\psi_{nm}(\mathbf{r})}{\sqrt{N_{nm}}},$
(44)

where

$\phi^{(\textup{X})}_{\mathbf{K}}(\mathbf{r},\mathbf{R})=\frac{1}{\sqrt{2\pi\hbar}}\mathrm{exp}\bigg{(}\frac{i\gamma}{2\hbar}\,\mathbf{K}\mathbf{r}\bigg{)}\\\
\times\mathrm{exp}\bigg{[}\frac{i}{\hbar}\bigg{(}\mathbf{K}-\frac{e}{2c}\,\mathbf{B}\times\mathbf{r}\bigg{)}\mathbf{R}\bigg{]},$
(45)

and the relative-motion wave function has the form

$\psi_{nm}(\mathbf{r})=\frac{1}{2^{(n+m)/2}}\bigg{(}\xi-\partialderivative{\xi^{*}}\bigg{)}^{n}\bigg{(}\xi^{*}-\partialderivative{\xi}\bigg{)}^{m}\Phi^{(e)}_{0}(\mathbf{r}).$
(46)

The part of Eq. (45) concerning the center-of-mass motion already includes the
necessary normalization factor, so $N_{nm}$ should be evaluated as

$\displaystyle
N_{nm}=\int\\!\differential{\mathbf{r}}\psi_{nm}^{*}(\mathbf{r})\psi_{nm}(\mathbf{r}).$
(47)

Applying operators $\hat{A}^{\dagger}_{e}$ and $\hat{B}^{\dagger}_{e}$ to the
vacuum state according to Eq. (46), we obtain the following structure of the
wave function $\psi_{nm}$:

$\displaystyle\psi_{nm}(\mathbf{r})=\sum_{a_{1},a_{2}}B^{nm}_{a_{1}a_{2}}(\xi^{*})^{a_{1}}(\xi)^{a_{2}}\Phi^{(e)}_{0}(\mathbf{r}).$
(48)

The coefficients $B^{nm}_{a_{1}a_{2}}$ can be found directly by means of Eq.
(46). We have implemented the corresponding procedure within the Wolfram
Mathematica subroutine, which is included in the Supplemental Material (for
the details see Appendix I).

Taking into account Eq. (48), the normalization factor can be calculated as
follows (we do not display the summation indices explicitly):

$N_{nm}=\sum\sum B^{nm*}_{a_{1}^{\prime}a_{2}^{\prime}}B^{nm}_{a_{1}a_{2}}\\\
\times\int\\!\differential{\mathbf{r}}\big{[}\Phi^{(e)}_{0}(\mathbf{r})\big{]}^{*}\Phi^{(e)}_{0}(\mathbf{r})(\xi^{*})^{a_{1}+a_{2}^{\prime}}(\xi)^{a_{2}+a_{1}^{\prime}}\\\
=\sum\sum
B^{nm*}_{a_{1}^{\prime}a_{2}^{\prime}}B^{nm}_{a_{1}a_{2}}I^{(\textup{X},\textup{N})}_{a_{1}+a_{2}^{\prime},a_{2}+a_{1}^{\prime}}.$
(49)

Here we have defined the normalization master integral (MI)
$I^{(\textup{X},\textup{N})}_{p_{1},p_{2}}$ involving a monomial:

$\displaystyle
I^{(\textup{X},\textup{N})}_{p_{1},p_{2}}=\int\\!\differential{\mathbf{r}}\big{[}\Phi^{(e)}_{0}(\mathbf{r})\big{]}^{*}\Phi^{(e)}_{0}(\mathbf{r})(\xi^{*})^{p_{1}}(\xi)^{p_{2}}=\dfrac{1}{2\pi\lambda^{2}}\int\limits_{0}^{\infty}\differential{r}r\mathrm{e}^{-r^{2}/(2\lambda^{2})}\left(\dfrac{r}{2\lambda}\right)^{p_{1}+p_{2}}\int\limits_{0}^{2\pi}\differential{\varphi}\mathrm{e}^{{i\mkern
1.0mu}\varphi(p_{2}-p_{1})}=\dfrac{p_{1}!}{2^{p_{1}}}\,\delta_{p_{1}p_{2}}.$
(50)

With the aid of the normalization factor (49), one can construct now the
normalized wave functions

$\displaystyle\varphi^{(e)}_{nm}(\mathbf{r})=\dfrac{\psi_{nm}(\mathbf{r})}{\sqrt{N_{nm}}},$
(51)

which coincide with the expression (29), and also obtain the full two-particle
wave function (44).

Having described how to calculate the normalization factor $N_{nm}$ using the
normalization MI (50), we can proceed in the same way and compute the matrix
element of the electron-hole interaction
$V_{\mathbf{K},n^{\prime}m^{\prime};nm}$, defined in Eq. (41):

$V_{\mathbf{K},n^{\prime}m^{\prime};nm}=-\frac{e^{2}}{\lambda}\int\differential{\mathbf{r}}\big{[}\varphi^{(e)}_{n^{\prime}m^{\prime}}(\mathbf{r})\big{]}^{*}V(|\mathbf{r}+\mathbf{r}_{0}|)\varphi^{(e)}_{nm}(\mathbf{r})\\\
=\dfrac{\sum\sum
B^{nm*}_{a_{1}^{\prime}a_{2}^{\prime}}B^{nm}_{a_{1}a_{2}}I^{(\textup{X},eh)}_{a_{1}+a_{2}^{\prime},a_{2}+a_{1}^{\prime}}}{\sqrt{N_{n^{\prime}m^{\prime}}}\sqrt{N_{nm}}}.$
(52)

Here we have introduced the electron-hole interaction MI
$I^{(\textup{X},eh)}_{p_{1}p_{2}}$:

$I^{(\textup{X},eh)}_{p_{1}p_{2}}=-\frac{e^{2}}{\lambda}\int\\!\differential{\mathbf{r}}\big{[}\Phi^{(e)}_{0}(\mathbf{r})\big{]}^{*}\Phi^{(e)}_{0}(\mathbf{r})V(|\mathbf{r}+\mathbf{r}_{0}|)(\xi^{*})^{p_{1}}(\xi)^{p_{2}}=-\dfrac{1}{2\pi}\dfrac{1}{2^{p_{1}+p_{2}}}\frac{e^{2}}{\lambda}\left(\dfrac{r_{0}}{\lambda}\right)^{2+p_{1}+p_{2}}\\\
\times\int\limits_{0}^{\infty}\differential{x}x^{p_{1}+p_{2}+1}\mathrm{exp}\bigg{[}{-\left(\dfrac{r_{0}}{\lambda}\right)^{2}\dfrac{x^{2}}{2}}\bigg{]}\int\limits_{0}^{2\pi}\differential{\varphi}\mathrm{e}^{{i\mkern
1.0mu}\varphi(p_{2}-p_{1})}V(r_{0}\sqrt{1+x^{2}+2x\cos{\varphi}}),$ (53)

which in the case of the Coulomb potential $V(r)=\lambda/(\varepsilon r)$
reads

$\displaystyle
I^{(\textup{X},eh,\textup{C})}_{p_{1}p_{2}}=-\dfrac{E_{0}}{\pi\sqrt{2\pi}}\dfrac{1}{2^{p_{1}+p_{2}}}\left(\dfrac{r_{0}}{\lambda}\right)^{1+p_{1}+p_{2}}\int\limits_{0}^{\infty}\differential{x}x^{p_{1}+p_{2}+1}\mathrm{exp}\bigg{[}{-\left(\dfrac{r_{0}}{\lambda}\right)^{2}\dfrac{x^{2}}{2}}\bigg{]}\int\limits_{0}^{2\pi}\differential{\varphi}\dfrac{\exp[{i\mkern
1.0mu}\varphi(p_{2}-p_{1})]}{\sqrt{1+x^{2}+2x\cos{\varphi}}}.$ (54)

Here $E_{0}=\sqrt{\pi/2}(e^{2}/\varepsilon\lambda)$. With this expression in
hands, we are now able to calculate an arbitrary matrix element (30).

##  Appendix C Trion (T) in a magnetic field

###  C.1 Spectrum calculation recipe

In contrast to exciton, this three-particle system is not electrically
neutral, so one can only partially separate the center-of-mass degrees of
freedom (for the details see the next section). Instead of a conserved vector
$\mathbf{K}$, one has in this case an integer quantum number $k$ which can be
set to zero due to the degeneracy of the energy levels. The matrix elements
read

$\displaystyle
H^{(\textup{T})}_{n_{1}n_{2}n_{h}ml;n_{1}^{\prime}n_{2}^{\prime}n_{h}^{\prime}m^{\prime}l^{\prime}}$
$\displaystyle=E^{(\textup{T},0)}_{n_{1}+n_{2},n_{h}}\delta_{n_{1}n_{1}^{\prime}}\delta_{n_{2}n_{2}^{\prime}}\delta_{n_{h}n_{h}^{\prime}}\delta_{mm^{\prime}}\delta_{ll^{\prime}}$
$\displaystyle+V^{(eh)}_{n_{1}n_{2}n_{h}ml;n_{1}^{\prime}n_{2}^{\prime}n_{h}^{\prime}m^{\prime}l^{\prime}}$
$\displaystyle+V^{(ee)}_{n_{1}n_{2}n_{h}ml;n_{1}^{\prime}n_{2}^{\prime}n_{h}^{\prime}m^{\prime}l^{\prime}},$
(55)

where each index is a non-negative integer, $E^{(\textup{T},0)}_{n_{e},n_{h}}$
is given in Eq. (95). The matrix elements of the interparticle interaction are
displayed in Eqs. (96) and (108). Note that in what follows, we will compute
the whole matrix by independently varying all of the indices and taking into
account the Coulomb mixing of the free LLs, i.e. incorporating the $V$ terms
in Eq. (16) with $n_{1}\neq n_{1}^{\prime}$, $n_{2}\neq n_{2}^{\prime}$, and
$n_{h}\neq n_{h}^{\prime}$.

###  C.2 Details

A three-particle system involving two electrons with positions
$\mathbf{r}_{1}$ and $\mathbf{r}_{2}$ and a hole with position
$\mathbf{r}_{h}$ is described by the following Hamiltonian:

$\hat{H}^{(\textup{T})}=\hat{H}_{0}+\hat{V}_{eh}+\hat{V}_{ee},$ (56)

where

$\displaystyle\hat{H}_{0}$
$\displaystyle=\sum_{j=1,2}\frac{\hat{\boldsymbol{\pi}}_{j}^{2}}{2m_{e}}+\frac{\hat{\boldsymbol{\pi}}_{h}^{2}}{2m_{h}},$
(57) $\displaystyle\hat{V}_{eh}$
$\displaystyle=-\sum_{j=1,2}\frac{e^{2}}{\lambda}V(|\mathbf{r}_{j}-\mathbf{r}_{h}|),$
(58) $\displaystyle\hat{V}_{ee}$
$\displaystyle=\frac{e^{2}}{\lambda}V(|\mathbf{r}_{1}-\mathbf{r}_{2}|).$ (59)

Here
$\hat{\boldsymbol{\pi}}_{j}=-i\hbar\boldsymbol{\nabla}_{j}-(e/c)\boldsymbol{A}(\boldsymbol{r}_{j})$.
In what follows, we will again construct first the eigenfunctions of the
Hamiltonian $\hat{H}_{0}$ without the interparticle interaction and then use
the corresponding wave functions as a basis to subsequently diagonalize the
total Hamiltonian $\hat{H}^{(\textup{T})}$.

The most straightforward approach is to employ the one-particle solutions (26)
and simply multiply the functions depending on $\mathbf{r}_{1}$,
$\mathbf{r}_{2}$, and $\mathbf{r}_{h}$, respectively. However, this basis set
completely disregards the translational symmetry of the system as the
functions have only a well-defined projection of the total angular momentum
$M_{z}$. First, let us inspect a three-particle generalization of the vector
operator (32):

$\hat{\mathbf{K}}=\sum_{j=1,2}\hat{\boldsymbol{\pi}}_{j}+\hat{\boldsymbol{\pi}}_{h}+\frac{e}{c}\sum_{j=1,2}\mathbf{B}\times\mathbf{r}_{j}-\frac{e}{c}\,\mathbf{B}\times\mathbf{r}_{h}.$
(60)

This operator (so-called operator of the magnetic translations Lerner and
Lozovik (1980); Dzyubenko and Sivachenko (2000); Dzyubenko (2000, 2001))
commutes with both $\hat{H}_{0}$ and $\hat{H}^{(\textup{T})}$, but the crucial
difference between an exciton and trion is that in the latter case the $x$ and
$y$ components of $\hat{\mathbf{K}}$ do not commute with each other. It turns
out that $[\hat{K}_{x},\,\hat{K}_{y}]=i|e|\hbar B/c$ and generally this
commutator is proportional to the total electric charge of the system (an
exciton is neutral). This means that one can take advantage only of one
integral of motion (for instance, either $K_{x}$ or $K_{y}$). In fact, it is
convenient to consider the operator $\hat{\mathbf{K}}^{2}$ Dzyubenko (2000).
It obviously commutes with the Hamiltonian. Let us normalize the vector
operator according to

$\hat{\mathbf{k}}=\sqrt{\frac{c}{|e|\hbar B}}\,\hat{\mathbf{K}},$ (61)

so that $[\hat{k}_{x},\,\hat{k}_{y}]=i$. Introducing then operator
$\hat{\varkappa}=(i\hat{k}_{x}-\hat{k}_{y})/\sqrt{2}$, one obtains

$\displaystyle[\hat{\varkappa},\,\hat{\varkappa}^{\dagger}]=1,$ (62)
$\displaystyle\hat{\mathbf{K}}^{2}=2\hat{\varkappa}^{\dagger}\hat{\varkappa}+1.$
(63)

Therefore, the eigenvalues of $\hat{\mathbf{K}}^{2}$ are $2k+1$, where $k=0$,
$1$, … The solutions of the Schrödinger equation can be characterized by the
quantum number $k$, which allows one to partially separate the center-of-mass
motion. The question now is how to use this symmetry in practical
calculations, i.e., how to construct the eigenfunctions of
$\hat{\mathbf{K}}^{2}$.

Let us briefly outline the procedure proposed in Ref. Dzyubenko (2000). First,
we perform the following coordinate transformation:

$\displaystyle\mathbf{R}$
$\displaystyle=\frac{1}{\sqrt{2}}\,\big{(}\mathbf{r}_{1}+\mathbf{r}_{2}\big{)},\quad\mathbf{r}=\frac{1}{\sqrt{2}}\,\big{(}\mathbf{r}_{1}-\mathbf{r}_{2}\big{)}.$
(64)

Introducing then the three-particle vacuum state

$\Phi^{(\textup{T})}_{0}(\mathbf{r},\mathbf{R},\mathbf{r}_{h})=\Phi^{(e)}_{0}(\mathbf{r})\Phi^{(e)}_{0}(\mathbf{R})\Phi^{(e)}_{0}(\mathbf{r}_{h}),$
(65)

one can generate naive solutions by the corresponding creation operators of
the form (24) and (25) and obtain:

$\varphi^{(\textup{T})}_{n_{1}m_{1}n_{2}m_{2}n_{h}m_{h}}(\mathbf{r},\mathbf{R},\mathbf{r}_{h})\\\
=\varphi^{(e)}_{n_{1}m_{1}}(\mathbf{r})\varphi^{(e)}_{n_{2}m_{2}}(\mathbf{R})\varphi^{(h)}_{n_{h}m_{h}}(\mathbf{r}_{h}),$
(66)

where $\varphi^{(h)}_{nm}(\mathbf{r})=[\varphi^{(e)}_{nm}(\mathbf{r})]^{*}$.
The $z$ projection of the angular momentum equals
$M_{z}=(n_{1}-m_{1})+(n_{2}-m_{2})-(n_{h}-m_{h})$. Moreover, the reflection
$\mathbf{r}\to-\mathbf{r}$, which is equivalent to exchanging the positions of
the electrons, yields the factor $(-1)^{n_{1}-m_{1}}$, so the spin wave
function of the electron subsystem (not displayed here explicitly) is even
(odd) if $n_{1}-m_{1}$ is odd (even). In other words, the total electron spin
is $S_{e}=0$ (spin-singlet state) for even $n_{1}-m_{1}$ and $S_{e}=1$ (spin-
triplet state) for odd $n_{1}-m_{1}$. As was shown in Ref. Dzyubenko (2000),
the functions (66) can be unitary transformed, so that they become
eigenfunctions of $\hat{\mathbf{K}}^{2}$, i.e., they gain additional quantum
number $k$ preserving $M_{z}$ and $S_{e}$.

To perform this unitary transformation, let us again introduce dimensionless
variables $\xi=z/(2\lambda)$, $\xi_{R}=z_{R}/(2\lambda)$, and
$\xi_{h}=z_{h}/(2\lambda)$ and consider the following creation operators:

$\displaystyle\hat{A}^{\dagger}_{e}(\mathbf{r})$
$\displaystyle=\dfrac{\left(\xi-\partialderivative{\xi^{*}}\right)}{\sqrt{2}},\
\
\hat{B}^{\dagger}_{e}(\mathbf{r})=\dfrac{\left(\xi^{*}-\partialderivative{\xi}\right)}{\sqrt{2}},$
(67) $\displaystyle\hat{A}^{\dagger}_{e}(\mathbf{R})$
$\displaystyle=\dfrac{\left(\xi_{R}-\partialderivative{\xi_{R}^{*}}\right)}{\sqrt{2}},\
\
\hat{B}^{\dagger}_{e}(\mathbf{R})=\dfrac{\left(\xi_{R}^{*}-\partialderivative{\xi_{R}}\right)}{\sqrt{2}},$
(68) $\displaystyle\hat{A}^{\dagger}_{h}(\mathbf{r}_{h})$
$\displaystyle=\dfrac{\left(\xi_{h}^{*}-\partialderivative{\xi_{h}}\right)}{\sqrt{2}},\
\
\hat{B}^{\dagger}_{h}(\mathbf{r}_{h})=\dfrac{\left(\xi_{h}-\partialderivative{\xi_{h}^{*}}\right)}{\sqrt{2}}.$
(69)

Let us note that

$\hat{\varkappa}^{\dagger}=\sqrt{2}\hat{B}^{\dagger}_{e}(\mathbf{R})-\hat{B}_{h}(\mathbf{r}_{h}).$
(70)

This relation can be interpreted as a Bogolyubov transformation of
$B^{\dagger}_{e}(\mathbf{R})$:

$\hat{\widetilde{B}}^{\dagger}_{e}(\mathbf{R})=\hat{\varkappa}^{\dagger}=\sqrt{2}\hat{B}^{\dagger}_{e}(\mathbf{R})-\hat{B}_{h}(\mathbf{r}_{h})=\hat{S}\hat{B}^{\dagger}_{e}(\mathbf{R})\hat{S}^{\dagger},$
(71)

The corresponding unitary operator is

$\displaystyle\hat{S}=\exp{-\Theta\left[\hat{B}_{e}(\mathbf{R})\hat{B}_{h}(\mathbf{r}_{h})-\hat{B}^{\dagger}_{h}(\mathbf{r}_{h})\hat{B}^{\dagger}_{e}(\mathbf{R})\right]},$
(72)

where $\sinh\Theta=1$ and $\cosh\Theta=\sqrt{2}$. Another operator that does
not commute with $\hat{S}$ and should be transformed is
$\hat{B}^{\dagger}_{h}(\mathbf{r}_{h})$:

$\hat{\widetilde{B}}^{\dagger}_{h}(\mathbf{r}_{h})=\sqrt{2}\hat{B}^{\dagger}_{h}(\mathbf{r}_{h})-\hat{B}_{e}(\mathbf{R}).$
(73)

In addition to the operators, the wave functions themselves also have to be
transformed. Thus, up to the normalization factor the new basis is

$\displaystyle\\!\\!\\!\\!\hat{A}^{{\dagger}}_{e}(\mathbf{r})^{n_{1}}\\!\hat{A}^{\dagger}_{e}(\mathbf{R})^{n_{2}}\\!\hat{A}^{\dagger}_{h}(\mathbf{r}_{h})^{n_{h}}\\!\hat{\widetilde{B}}^{\dagger}_{e}(\mathbf{R})^{k}\\!\hat{B}^{\dagger}_{e}(\mathbf{r})^{m}\\!\hat{\widetilde{B}}_{h}^{\dagger}(\mathbf{r}_{h})^{l}\Psi^{(\textup{T})}_{0},\\!\\!\\!$
(74)

where we used a new vacuum state
$\Psi^{(\textup{T})}_{0}=\hat{S}\Phi^{(\textup{T})}_{0}$. In this basis, three
well-defined quantum numbers are $k$, $M_{z}$, and $S_{e}$, i.e. the
Hamiltonian is block-diagonal in terms of these quantum numbers. The
projection of the total angular momentum reads

$\displaystyle M_{z}=\underbrace{n_{1}-m}_{\textup{from
$\mathbf{r}$}}+\underbrace{n_{2}-k}_{\textup{from
$\mathbf{R}$}}-\underbrace{(n_{h}-l)}_{\textup{from $\mathbf{r}_{h}$}}.$ (75)

The spin quantum number is $S_{e}=0$ ($S_{e}=1$) if $n_{1}-m$ is even (odd).
Moreover, hereinafter we will assume $k=0$ since the Hamiltonian
$H^{(\textup{T})}$ commutes with $\widetilde{B}^{\dagger}_{e}(\mathbf{R})$.
The matrix elements do not depend on $k$, so the energy also turns out to be
$k$-independent.

In order to obtain an explicit form of the new vacuum state, it is necessary
to rewrite the expression for the operator $\hat{S}$. Using the
Baker–Campbell–Hausdorff formula (see, e.g., Ref. Kirzhnits (1967)) for the
operator equality
$\mathrm{exp}[\Theta(\hat{a}+\hat{b})]=\mathrm{exp}(\Theta\hat{d})$, one can
present the operator $\hat{S}$ in the following form:

$\hat{S}=\mathrm{exp}\left[\tanh\Theta\hat{B}^{\dagger}_{h}(\mathbf{r}_{h})\hat{B}^{\dagger}_{e}(\mathbf{R})\right]\\\
\times\exp{-\ln{(\cosh\Theta)}\\!\left[\hat{B}^{\dagger}_{e}(\mathbf{R})\hat{B}_{e}(\mathbf{R})+\hat{B}^{\dagger}_{h}(\mathbf{r}_{h})\hat{B}_{h}(\mathbf{r}_{h})+1\right]}\\\
\times\mathrm{exp}\left[-\tanh\Theta\hat{B}_{e}(\mathbf{R})\hat{B}_{h}(\mathbf{r}_{h})\right].$
(76)

The new vacuum state then reads

$\displaystyle\Psi^{(\textup{T})}_{0}=\hat{S}\Phi^{(\textup{T})}_{0}=\dfrac{1}{\cosh\Theta}\mathrm{exp}\left[\tanh\Theta\hat{B}^{\dagger}_{h}(\mathbf{r}_{h})\hat{B}^{\dagger}_{e}(\mathbf{R})\right]\Phi^{(\textup{T})}_{0}.$
(77)

The state (65) is given by

$\displaystyle\Phi^{(\textup{T})}_{0}(\mathbf{r},\mathbf{R},\mathbf{r}_{h})=\dfrac{\mathrm{exp}(-\xi\xi^{*}-\xi_{R}\xi^{*}_{R}-\xi_{h}\xi^{*}_{h})}{\left(2\pi\right)^{3/2}\lambda^{3}}.$
(78)

Taking into account the explicit form of the operators
$B^{\dagger}_{h}(\mathbf{r}_{h})$ and $B^{\dagger}_{e}(\mathbf{R})$, one
obtains

$\Psi^{(\textup{T})}_{0}=\exp[\dfrac{1}{2\sqrt{2}}\left(\xi_{h}-\partialderivative{\xi^{*}_{h}}\right)\left(\xi^{*}_{R}-\partialderivative{\xi_{R}}\right)]\\\
\quad\quad\times\dfrac{\exp(-\xi\xi^{*}-\xi_{R}\xi^{*}_{R}-\xi_{h}\xi^{*}_{h})}{\sqrt{2}\left(2\pi\right)^{3/2}\lambda^{3}}.$
(79)

Let us label the function
$\sqrt{2}\Psi^{(\textup{T})}_{0}(\mathbf{r},\mathbf{R},\mathbf{r}_{h})$ as a
new normalized vacuum state
$\tilde{\Phi}^{(\textup{T})}_{0}(\mathbf{r},\mathbf{R},\mathbf{r}_{h})$, which
reads

$\displaystyle\tilde{\Phi}^{(\textup{T})}_{0}\\!(\mathbf{r},\\!\mathbf{R},\mathbf{r}_{h})\\!=\\!\dfrac{\exp{-\xi\xi^{*}\\!-\\!\xi_{R}\xi^{*}_{R}\\!-\\!\xi_{h}\xi^{*}_{h}\\!+\\!\sqrt{2}\xi_{h}\xi^{*}_{R}}}{\sqrt{2}\left(2\pi\right)^{3/2}\lambda^{3}}.\\!\\!$
(80)

Now one can generate an arbitrary wave function by means of the creation
operators in Eq. (74). Regarding the normalization of the wave functions, we
prefer to keep only factor $1/[\left(2\pi\right)^{3/2}\lambda^{3}]$ explicitly
in contrast to other normalization factors which will be taken into account at
the very end. Thus, we first employ the operator
$\hat{B}^{\dagger}_{e}(\mathbf{r})^{m}\hat{\widetilde{B}}_{h}^{\dagger}(\mathbf{r}_{h})^{l}$.
Note that

$\displaystyle\hat{\widetilde{B}}^{\dagger}_{h}(\mathbf{r}_{h})^{l}\tilde{\Phi}^{(\textup{T})}_{0}\\!(\mathbf{r},\\!\mathbf{R},\mathbf{r}_{h})\\!=\dfrac{\left(\xi_{h}\right)^{l}\tilde{\Phi}^{(\textup{T})}_{0}\\!(\mathbf{r},\\!\mathbf{R},\mathbf{r}_{h})\\!}{\sqrt{2}}.$
(81)

Then, we obtain

$\hat{B}^{\dagger}_{e}(\mathbf{r})^{m}\hat{\widetilde{B}}^{\dagger}_{h}(\mathbf{r}_{h})^{l}\tilde{\Phi}^{(\textup{T})}_{0}\\!(\mathbf{r},\\!\mathbf{R},\mathbf{r}_{h})\\\
=2^{(m-1)/2}\left(\xi^{*}\right)^{m}\left(\xi_{h}\right)^{l}\tilde{\Phi}^{(\textup{T})}_{0}\\!(\mathbf{r},\\!\mathbf{R},\mathbf{r}_{h}).$
(82)

It is convenient to introduce the following function:

$\displaystyle\psi_{ml}(\mathbf{r},\mathbf{R},\mathbf{r}_{h})=\left(\xi^{*}\right)^{m}\left(\xi_{h}\right)^{l}\tilde{\Phi}^{(\textup{T})}_{0}\\!(\mathbf{r},\\!\mathbf{R},\mathbf{r}_{h}).\\!$
(83)

According to Eq. (74), the general wave function up to the normalization
factor reads

$\displaystyle\psi_{n_{1}n_{2}n_{h}ml}(\mathbf{r},\mathbf{R},\mathbf{r}_{h})$
$\displaystyle=\left[\dfrac{1}{\sqrt{2}}\left(\xi-\partialderivative{\xi^{*}}\right)\right]^{n_{1}}\left[\dfrac{1}{\sqrt{2}}\left(\xi_{R}-\partialderivative{\xi^{*}_{R}}\right)\right]^{n_{2}}\left[\dfrac{1}{\sqrt{2}}\left(\xi^{*}_{h}-\partialderivative{\xi_{h}}\right)\right]^{n_{h}}\psi_{ml}(\mathbf{r},\mathbf{R},\mathbf{r}_{h}).$
(84)

Let us now calculate the normalization factor:

$\displaystyle
N_{n_{1}n_{2}n_{h}ml}=\langle\psi_{n_{1}n_{2}n_{h}ml}\text{\textbar}\psi_{n_{1}n_{2}n_{h}ml}\rangle.$
(85)

This expression involves three 2D integrals (over $\mathbf{r}$, $\mathbf{R}$,
and $\mathbf{r}_{h})$. Any wave function can be represented as the following
linear combination involving products of
$\tilde{\Phi}^{(\textup{T})}_{0}(\mathbf{r},\mathbf{R},\mathbf{r}_{h})$ and
monomials:

$\displaystyle\psi_{n_{1}n_{2}n_{h}ml}(\mathbf{r},\mathbf{R},\mathbf{r}_{h})=\sum
C^{n_{1}n_{2}n_{h}ml}_{a_{1}a_{2}b_{1}b_{2}c_{1}c_{2}}(\xi^{*})^{a_{1}}(\xi)^{a_{2}}$
$\displaystyle\quad\quad\times(\xi_{h})^{b_{1}}(\xi_{h}^{*})^{b_{2}}(\xi_{R})^{c_{1}}(\xi_{R}^{*})^{c_{2}}\tilde{\Phi}^{(\textup{T})}_{0}\\!(\mathbf{r},\\!\mathbf{R},\mathbf{r}_{h}).$
(86)

Therefore,

$\displaystyle N_{n_{1}n_{2}n_{h}ml}$ $\displaystyle=\sum\sum
C^{n_{1}n_{2}n_{h}ml*}_{a^{\prime}_{1}a^{\prime}_{2}b^{\prime}_{1}b^{\prime}_{2}c^{\prime}_{1}c^{\prime}_{2}}C^{n_{1}n_{2}n_{h}ml}_{a_{1}a_{2}b_{1}b_{2}c_{1}c_{2}}\int\\!\differential{\mathbf{r}}\int\\!\differential{\mathbf{R}}\int\\!\differential{\mathbf{r}_{h}}\big{|}\tilde{\Phi}^{(\textup{T})}_{0}\\!(\mathbf{r},\\!\mathbf{R},\mathbf{r}_{h})\big{|}^{2}$
$\displaystyle\qquad\qquad\times(\xi^{*})^{a_{1}+a_{2}^{\prime}}(\xi)^{a_{2}+a_{1}^{\prime}}(\xi_{h})^{b_{1}+b_{2}^{\prime}}(\xi_{h}^{*})^{b_{2}+b_{1}^{\prime}}(\xi_{R})^{c_{1}+c_{2}^{\prime}}(\xi_{R}^{*})^{c_{2}+c_{1}^{\prime}}$
$\displaystyle=\sum\sum
C^{n_{1}n_{2}n_{h}ml*}_{a^{\prime}_{1}a^{\prime}_{2}b^{\prime}_{1}b^{\prime}_{2}c^{\prime}_{1}c^{\prime}_{2}}C^{n_{1}n_{2}n_{h}ml}_{a_{1}a_{2}b_{1}b_{2}c_{1}c_{2}}I^{(\textup{T},\textup{N})}_{a_{1}+a_{2}^{\prime},a_{2}+a_{1}^{\prime},b_{1}+b_{2}^{\prime},b_{2}+b_{1}^{\prime},c_{1}+c_{2}^{\prime},c_{2}+c_{1}^{\prime}},$
(87)

where we define the master integral
$I^{(\textup{T},\textup{N})}_{p_{1}p_{2}q_{1}q_{2}r_{1}r_{2}}$ as

$\displaystyle I^{(\textup{T},\textup{N})}_{p_{1}p_{2}q_{1}q_{2}r_{1}r_{2}}$
$\displaystyle=\frac{1}{(2\pi)^{3}\lambda^{6}}\int\\!\differential{\mathbf{r}}\int\\!\differential{\mathbf{R}}\int\\!\differential{\mathbf{r}_{h}}(\xi^{*})^{p_{1}}(\xi)^{p_{2}}(\xi_{h})^{q_{1}}(\xi_{h}^{*})^{q_{2}}(\xi_{R})^{r_{1}}(\xi_{R}^{*})^{r_{2}}$
$\displaystyle\times\mathrm{exp}\Big{[}-2(\xi\xi^{*}+\xi_{R}\xi^{*}_{R}+\xi_{h}\xi^{*}_{h})+\sqrt{2}(\xi_{h}\xi^{*}_{R}+\xi_{h}^{*}\xi_{R})\Big{]}.$
(88)

The expansion coefficients
$C^{n_{1}n_{2}n_{h}ml}_{a_{1}a_{2}b_{1}b_{2}c_{1}c_{2}}$ for an arbitrary wave
function in terms of monomials
$(\xi^{*})^{a_{1}}(\xi)^{a_{2}}(\xi_{h})^{b_{1}}(\xi_{h}^{*})^{b_{2}}(\xi_{R})^{c_{1}}(\xi_{R}^{*})^{c_{2}}$
can be obtained automatically. We attach the corresponding subroutine in
Wolfram Mathematica format. Thus, we should evaluate only the master integral
(88). First, taking into account $\xi=r\mathrm{e}^{{i\mkern
1.0mu}\varphi}/(2\lambda)$, we integrate out the $\mathbf{r}$-dependent part:

$\displaystyle\int\differential{\mathbf{r}}\left(\xi^{*}\right)^{p_{1}}\left(\xi\right)^{p_{2}}\exp(-2\xi\xi^{*})=\int\limits^{\infty}_{0}r\differential{r}\int\limits^{2\pi}_{0}\differential{\varphi}\exp[{i\mkern
1.0mu}\varphi(p_{2}-p_{1})]\left(\dfrac{r}{2\lambda}\right)^{p_{1}+p_{2}}\exp(-\dfrac{r^{2}}{2\lambda^{2}})=2\pi\delta_{p_{1}p_{2}}\dfrac{p_{1}!\lambda^{2}}{2^{p_{1}}}.$

Equation (88) takes the form

$I^{(\textup{T},\textup{N})}_{p_{1}p_{2}q_{1}q_{2}r_{1}r_{2}}=\dfrac{2\pi\delta_{p_{1}p_{2}}p_{1}!}{(2\pi)^{3}2^{p_{1}}\lambda^{4}}\int\limits^{\infty}_{0}R\differential{R}\int\limits^{\infty}_{0}r_{h}\differential{r_{h}}\left(\dfrac{r_{h}}{2\lambda}\right)^{q_{1}+q_{2}}\left(\dfrac{R}{2\lambda}\right)^{r_{1}+r_{2}}\int\limits_{0}^{2\pi}\differential{\varphi_{h}}\int\limits_{0}^{2\pi}\differential{\varphi_{R}}\exp[{i\mkern
1.0mu}\varphi_{h}(q_{1}-q_{2})]\exp[{i\mkern
1.0mu}\varphi_{R}(r_{1}-r_{2})]\\\
\times\exp[-\dfrac{R^{2}}{2\lambda^{2}}-\dfrac{r_{h}^{2}}{2\lambda^{2}}+\dfrac{r_{h}R\cos{(\varphi_{R}-\varphi_{h})}}{\sqrt{2}\lambda^{2}}].$
(89)

Using the Jacobi-Anger relation $\mathrm{exp}({i\mkern
1.0mu}z\cos{\varphi})=\\!\\!\\!\\!\sum\limits_{n=-\infty}^{\infty}\\!\\!\\!\\!{i\mkern
1.0mu}^{n}J_{n}(z)\mathrm{exp}({i\mkern 1.0mu}n\varphi)$, we can integrate
over the angles as follows:

$\displaystyle\int\limits_{0}^{2\pi}\differential{\varphi_{R}}\int\limits_{0}^{2\pi}\differential{\varphi_{h}}\exp[{i\mkern
1.0mu}\varphi_{h}(q_{1}-q_{2})]\exp[{i\mkern
1.0mu}\varphi_{R}(r_{1}-r_{2})]\exp[{i\mkern
1.0mu}n(\varphi_{R}-\varphi_{h})]=(2\pi)^{2}\delta_{0,q_{1}-q_{2}-n}\delta_{0,r_{1}-r_{2}+n}.$
(90)

Let us assume that $q_{1}\geq q_{2}$. One obtains

$\displaystyle I^{(\textup{T},\textup{N})}_{p_{1}p_{2}q_{1}q_{2}r_{1}r_{2}}$
$\displaystyle=\frac{1}{\lambda^{4}}\,\delta_{p_{1}p_{2}}\delta_{q_{1}-q_{2},r_{2}-r_{1}}\,\dfrac{p_{1}!}{2^{p_{1}}}\,{i\mkern
1.0mu}^{q_{1}-q_{2}}\int\limits^{\infty}_{0}\differential{R}R\left(\dfrac{R}{2\lambda}\right)^{r_{1}+r_{2}}\exp(-\dfrac{R^{2}}{2\lambda^{2}})$
$\displaystyle\times\int\limits^{\infty}_{0}\differential{r_{h}}r_{h}\left(\dfrac{r_{h}}{2\lambda}\right)^{q_{1}+q_{2}}\exp(-\dfrac{r_{h}}{2\lambda^{2}})J_{q_{1}-q_{2}}\left(-{i\mkern
1.0mu}\dfrac{r_{h}R}{\lambda^{2}\sqrt{2}}\right)=\dfrac{\delta_{p_{1}p_{2}}\delta_{q_{1}-q_{2},r_{2}-r_{1}}p_{1}!}{2^{p_{1}+q_{1}+r_{1}+\frac{q_{1}-q_{2}}{2}-2}}\int\limits^{\infty}_{0}\differential{\left(\dfrac{R}{\sqrt{2}\lambda}\right)}\left(\dfrac{R}{\sqrt{2}\lambda}\right)$
$\displaystyle\times\int\limits^{\infty}_{0}\differential{\left(\dfrac{r_{h}}{\sqrt{2}\lambda}\right)}\left(\dfrac{r_{h}}{\sqrt{2}\lambda}\right)\sum\limits^{\infty}_{m=0}\dfrac{1}{2^{m}m!}\dfrac{1}{(q_{1}-q_{2}+m)!}\exp(-\dfrac{r^{2}_{h}}{2\lambda^{2}})\exp(-\dfrac{R^{2}}{2\lambda^{2}})\left(\dfrac{r_{h}}{\sqrt{2}\lambda}\right)^{2q_{1}+2m}$
$\displaystyle\times\left(\dfrac{R}{\sqrt{2}\lambda}\right)^{2r_{1}+2m+2q_{1}-2q_{2}}=\dfrac{\delta_{p_{1}p_{2}}\delta_{q_{1}-q_{2},r_{2}-r_{1}}p_{1}!}{2^{p_{1}+q_{1}+r_{1}+\frac{q_{1}-q_{2}}{2}}}\int\limits^{\infty}_{0}\differential{x}\int\limits^{\infty}_{0}\differential{y}\sum\limits^{\infty}_{m=0}\dfrac{y^{q_{1}+m}\exp(-y)x^{m+r_{1}+q_{1}-q_{2}}\exp(-x)}{2^{m}m!\left(q_{1}-q_{2}+m\right)!}=$
$\displaystyle=\dfrac{\delta_{p_{1}p_{2}}\delta_{q_{1}-q_{2},r_{2}-r_{1}}p_{1}!}{2^{p_{1}+q_{1}+r_{1}+\frac{q_{1}-q_{2}}{2}}}\sum\limits_{m=0}^{\infty}\dfrac{\left(q_{1}+m\right)!(r_{1}+m+q_{1}-q_{2})!}{m!2^{m}\left(q_{1}-q_{2}+m\right)!}=\dfrac{\delta_{p_{1}p_{2}}\delta_{q_{1}-q_{2},r_{2}-r_{1}}p_{1}!q_{1}!r_{2}!}{(q_{1}-q_{2})!2^{p_{1}+q_{1}+r_{1}+\frac{q_{1}-q_{2}}{2}}}$
$\displaystyle\times{}_{2}F_{1}\left(q_{1}+1,r_{2}+1,q_{1}-q_{2}+1;\dfrac{1}{2}\right).$
(91)

If $q_{1}<q_{2}$, then one can perform complex conjugation of the initial
expression for $I^{(\textup{T},\textup{N})}_{p_{1}p_{2}q_{1}q_{2}r_{1}r_{2}}$
and calculate it via the same steps as described above. The result can also be
obtained by simply substituting $q_{1}\leftrightarrow q_{2}$ and
$r_{1}\leftrightarrow r_{2}$. Then, for universality, we introduce
$q_{>}=\textup{max}(q_{1},q_{2})$, $r_{>}=\textup{max}(r_{1},r_{2})$, and
$r_{<}=\textup{min}(r_{1},r_{2})$. In terms of these quantities, the master
integral reads

$\displaystyle
I^{(\textup{T},\textup{N})}_{p_{1}p_{2}q_{1}q_{2}r_{1}r_{2}}=\\!\dfrac{\delta_{p_{1}p_{2}}\delta_{q_{1}-q_{2},r_{2}-r_{1}}p_{1}!q_{>}!r_{>}!}{|q_{1}-q_{2}|!2^{p_{1}+q_{>}+r_{<}+\frac{|q_{1}-q_{2}|}{2}}}\\!\
_{2}F_{1}\\!\\!\left(\\!q_{>}\\!+\\!1,r_{>}\\!+\\!1,|q_{1}\\!-\\!q_{2}|\\!+1;\\!\dfrac{1}{2}\\!\right).$
(92)

Having obtained the expression for the normalization factor
$N_{n_{1}n_{2}n_{h}ml}$, we introduce the normalized wave functions as

$\displaystyle\varphi^{(\textup{T})}_{n_{1}n_{2}n_{h}ml}(\mathbf{r},\mathbf{R},\mathbf{r}_{h})=\dfrac{\psi_{n_{1}n_{2}n_{h}ml}(\mathbf{r},\mathbf{R},\mathbf{r}_{h})}{\sqrt{N_{n_{1}n_{2}n_{h}ml}}}.$
(93)

With these wave functions, we can calculate the matrix elements of interest.
First, it immediately follows that the matrix elements of the operator
$\hat{H}_{0}$ read

$\langle\varphi^{(\textup{T})}_{n_{1}n_{2}n_{h}ml}\text{\textbar}\hat{H}_{0}\text{\textbar}\varphi^{(\textup{T})}_{n_{1}^{\prime}n_{2}^{\prime}n_{h}^{\prime}m^{\prime}l^{\prime}}\rangle\\\
=\delta_{n_{1}n_{1}^{\prime}}\delta_{n_{2}n_{2}^{\prime}}\delta_{n_{h}n_{h}^{\prime}}\delta_{mm^{\prime}}\delta_{ll^{\prime}}E^{(\textup{T},0)}_{n_{1}+n_{2},n_{h}},$
(94)

where

$E^{(\textup{T},0)}_{n_{e}n_{h}}=\hbar\omega_{e}\bigg{(}1+n_{e}\bigg{)}+\hbar\omega_{h}\bigg{(}\frac{1}{2}+n_{h}\bigg{)}.$
(95)

Next, one has to compute the matrix elements involving the interparticle-
interaction terms. Let us start with the contribution $\hat{V}_{eh}$. Taking
into account Eq. (C.2), we obtain

$\displaystyle
V^{(eh)}_{n_{1}n_{2}n_{h}ml;n_{1}^{\prime}n_{2}^{\prime}n_{h}^{\prime}m^{\prime}l^{\prime}}\\!\\!=\\!\\!\langle\varphi^{(\textup{T})}_{n_{1}n_{2}n_{h}ml}\text{\textbar}\hat{V}_{eh}\text{\textbar}\varphi^{(\textup{T})}_{n_{1}^{\prime}n_{2}^{\prime}n_{h}^{\prime}m^{\prime}l^{\prime}}\rangle\\!=\\!\\!\dfrac{\sum\sum
C^{n_{1}n_{2}n_{h}ml*}_{a^{\prime}_{1}a^{\prime}_{2}b^{\prime}_{1}b^{\prime}_{2}c^{\prime}_{1}c^{\prime}_{2}}C^{n_{1}n_{2}n_{h}ml}_{a_{1}a_{2}b_{1}b_{2}c_{1}c_{2}}I^{(\textup{T},eh)}_{a_{1}+a_{2}^{\prime},a_{2}+a_{1}^{\prime},b_{1}+b_{2}^{\prime},b_{2}+b_{1}^{\prime},c_{1}+c_{2}^{\prime},c_{2}+c_{1}^{\prime}}}{\sqrt{N_{n_{1}n_{2}n_{h}ml}}\sqrt{N_{n_{1}^{\prime}n_{2}^{\prime}n_{h}^{\prime}m^{\prime}l^{\prime}}}},$
(96)

where

$I^{(\textup{T},eh)}_{p_{1}p_{2}q_{1}q_{2}r_{1}r_{2}}=-\dfrac{e^{2}}{\lambda}\\!\int\\!\\!\differential{\mathbf{r}}\\!\\!\int\\!\\!\differential{\mathbf{R}}\\!\\!\int\\!\\!\differential{\mathbf{r}}_{h}(\xi^{*})^{p_{1}}(\xi)^{p_{2}}(\xi_{h})^{q_{1}}(\xi_{h}^{*})^{q_{2}}(\xi_{R})^{r_{1}}(\xi_{R}^{*})^{r_{2}}\\\
\times\dfrac{\exp{-2(\xi\xi^{*}+\xi_{R}\xi^{*}_{R}+\xi_{h}\xi^{*}_{h})+\sqrt{2}(\xi_{h}\xi^{*}_{R}+\xi_{h}^{*}\xi_{R})}}{\left(2\pi\right)^{3}\lambda^{6}}\left[V(|\mathbf{r}_{h}-\mathbf{r}_{1}|)+V(|\mathbf{r}_{h}-\mathbf{r}_{2}|)\right].$
(97)

It is convenient to introduce the variable
$\widetilde{\mathbf{r}}=\sqrt{2}\mathbf{r}_{h}-\mathbf{R}$. We perform further
calculations for an arbitrary spherically symmetric potential, while the
Coulomb potential will be considered only at the very end of our calculations.
In terms of $\widetilde{\mathbf{r}}$ , the sum in square brackets can be
written as
$V(|\widetilde{\mathbf{r}}-\mathbf{r}|/\sqrt{2})+V(|\widetilde{\mathbf{r}}+\mathbf{r}|/2)$.
Passing to the coordinates
$\\{\widetilde{\mathbf{r}},\mathbf{r},\mathbf{R}\\}$ and also making an
additional change of variables in the second term [reflection
$\mathbf{r}\to-\mathbf{r}$ leading to the factor $(-1)^{p_{1}+p_{2}}$], we
obtain the following expression:

$\displaystyle I^{(\textup{T},eh)}_{p_{1}p_{2}q_{1}q_{2}r_{1}r_{2}}$
$\displaystyle=-\dfrac{e^{2}[1+(-1)^{p_{1}+p_{2}}]}{2\lambda(2\pi)^{3}\lambda^{6}}\int\differential{\mathbf{r}}\int\differential{\mathbf{R}}\int\differential{\widetilde{\mathbf{r}}}(\xi)^{p_{2}}(\xi^{*})^{p_{1}}(\xi_{R})^{r_{1}}(\xi_{R}^{*})^{r_{2}}\left(\dfrac{\widetilde{\xi}+\xi_{R}}{\sqrt{2}}\right)^{q_{1}}\left(\dfrac{\widetilde{\xi}^{*}+\xi^{*}_{R}}{\sqrt{2}}\right)^{q_{2}}$
$\displaystyle\times
V\left(\dfrac{|\widetilde{\mathbf{r}}-\mathbf{r}|}{\sqrt{2}}\right)\exp(-\dfrac{r^{2}}{2\lambda^{2}}-\dfrac{\widetilde{r}^{2}}{4\lambda^{2}}-\dfrac{R^{2}}{4\lambda^{2}})=-\dfrac{e^{2}(1+(-1)^{p_{1}+p_{2}})}{2\lambda(2\pi)^{3}\lambda^{6}2^{(q_{1}+q_{2})/2}}\sum\limits_{k=0}^{q_{1}}\sum\limits_{n=0}^{q_{2}}C^{k}_{q_{1}}C^{n}_{q_{2}}\int\differential{\mathbf{r}}\int\differential{\widetilde{\mathbf{r}}}$
$\displaystyle\times
V\left(\dfrac{|\widetilde{\mathbf{r}}-\mathbf{r}|}{\sqrt{2}}\right)\exp(-\dfrac{r^{2}}{2\lambda^{2}}-\dfrac{\widetilde{r}^{2}}{4\lambda^{2}})(\xi)^{p_{2}}(\xi^{*})^{p_{1}}(\widetilde{\xi})^{q_{1}-k}(\widetilde{\xi}^{*})^{q_{2}-n}\int\differential{\mathbf{R}}(\xi_{R})^{r_{1}+k}(\xi^{*}_{R})^{r_{2}+n}\exp(-\dfrac{R^{2}}{4\lambda^{2}}),$
(98)

where $\widetilde{\xi}=\sqrt{2}\xi_{h}-\xi_{R}$. Let us separately consider
$\mathbf{R}$ and $\widetilde{\mathbf{r}}$–$\mathbf{r}$ parts. For the
$\mathbf{R}$ integration, we obtain the following expression:

$\displaystyle\int\limits^{\infty}_{0}\differential{R}R\exp(-\dfrac{R^{2}}{4\lambda^{2}})\left(\dfrac{R}{2\lambda}\right)^{r_{1}+r_{2}+k+n}\int\limits^{2\pi}_{0}\differential{\varphi_{R}}\\!\exp{{i\mkern
1.0mu}\varphi_{R}[r_{1}+k-(r_{2}+n)]}=4\pi\lambda^{2}\delta_{r_{1}+k,r_{2}+n}(k+r_{1})!.$
(99)

Now we consider the $\widetilde{\mathbf{r}}$–$\mathbf{r}$ integration in Eq.
(98). In order to factorize the expression, it suffices here to perform the
following change of variables: $\mathbf{r}=\bm{\rho}_{1}/3+\bm{\rho}_{2}/3$,
$\widetilde{\mathbf{r}}=\bm{\rho}_{1}/3-2\bm{\rho}_{2}/3$. Then, the
corresponding integral reads

$\displaystyle\int\differential{\mathbf{r}}\int\differential{\widetilde{\mathbf{r}}}V\left(\dfrac{|\widetilde{\mathbf{r}}-\mathbf{r}|}{\sqrt{2}}\right)\exp(-\dfrac{r^{2}}{2\lambda^{2}}-\dfrac{\widetilde{r}^{2}}{4\lambda^{2}})(\xi)^{p_{2}}(\xi^{*})^{p_{1}}(\widetilde{\xi})^{q_{1}-k}(\widetilde{\xi}^{*})^{q_{2}-n}=\dfrac{1}{3^{2+p_{2}+p_{1}+q_{1}+q_{2}-k-n}}\int\differential\bm{\rho}_{1}\\!\\!\int\differential\bm{\rho}_{2}V\\!\\!\left(\dfrac{|\bm{\rho}_{2}|}{\sqrt{2}}\right)$
$\displaystyle\times\exp(-\dfrac{\rho_{1}^{2}+2\rho_{2}^{2}}{12\lambda^{2}})\left(\xi_{\rho_{1}}+\xi_{\rho_{2}}\right)^{p_{2}}\left(\xi^{*}_{\rho_{1}}+\xi_{\rho_{2}}^{*}\right)^{p_{1}}\left(\xi_{\rho_{1}}-2\xi_{\rho_{2}}\right)^{q_{1}-k}\left(\xi^{*}_{\rho_{1}}-2\xi^{*}_{\rho_{2}}\right)^{q_{2}-n}=\dfrac{1}{3^{2+p_{2}+p_{1}+q_{1}+q_{2}-k-n}}$
$\displaystyle\times\sum\limits_{t_{2}=0}^{p_{2}}\sum\limits_{t_{1}=0}^{p_{1}}\sum\limits_{s_{1}=0}^{q_{1}-k}\sum\limits_{s_{2}=0}^{q_{2}-n}C^{t_{2}}_{p_{2}}C^{t_{1}}_{p_{1}}C^{s_{1}}_{q_{1}-k}C^{s_{2}}_{q_{2}-n}\int\limits_{0}^{\infty}\rho_{2}\differential{\rho_{2}}V\left(\dfrac{\rho_{2}}{\sqrt{2}}\right)\exp(-\dfrac{\rho_{2}^{2}}{6\lambda^{2}})(\xi_{\rho_{2}})^{p_{2}-t_{2}}(\xi^{*}_{\rho_{2}})^{p_{1}-t_{1}}(-2\xi_{\rho_{2}})^{q_{1}-l-s_{1}}$
$\displaystyle\times(-2\xi^{*}_{\rho_{2}})^{q_{2}-n-s_{2}}\int\rho_{1}\differential{\rho_{1}}(\xi_{\rho_{1}})^{t_{2}+s_{1}}(\xi_{\rho_{1}}^{*})^{t_{1}+s_{2}}\exp(-\dfrac{\rho^{2}_{1}}{12\lambda^{2}})=(*).$
(100)

The $\rho_{1}$\- and $\rho_{2}$-integrals are factorized, so they can be
easily evaluated separately. First,

$\displaystyle\int\limits_{0}^{\infty}\differential{\rho_{1}}\rho_{1}\\!\left(\dfrac{\rho_{1}}{2\lambda}\right)^{t_{2}+t_{1}+s_{1}+s_{2}}\\!\\!\\!\\!\\!\exp(-\dfrac{\rho_{1}^{2}}{12\lambda^{2}})\\!\\!\int\limits^{2\pi}_{0}\\!\differential{\varphi_{\rho_{1}}}\\!\exp[{i\mkern
1.0mu}\varphi_{\rho_{1}}(t_{2}\\!+\\!s_{1}\\!-\\!(t_{1}\\!+\\!s_{2}))]=12\pi\lambda^{2}3^{t_{1}+s_{2}}(t_{1}+s_{2})!\delta_{t_{2}+s_{1},t_{1}+s_{2}}.$
(101)

Second,

$\displaystyle\int\limits_{0}^{\infty}\differential{\rho_{2}}\rho_{2}V\left(\dfrac{\rho_{2}}{\sqrt{2}}\right)\left(\dfrac{\rho_{2}}{2\lambda}\right)^{p_{2}-t_{2}+p_{1}-t_{1}+q_{1}-k-s_{1}+q_{2}-n-s_{2}}\exp(-\dfrac{\rho_{2}^{2}}{6\lambda^{2}})\int\limits^{2\pi}_{0}\differential{\varphi_{\rho_{2}}}(-2)^{q_{1}-k-s_{1}+q_{2}-n-s_{2}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\times\exp[{i\mkern
1.0mu}\varphi_{\rho_{2}}(p_{2}-t_{1}+q_{1}-k-s_{1}-(p_{1}-t_{1}+q_{2}-n-s_{2}))]=$
$\displaystyle=2\pi\delta_{p_{2}+q_{1}-k,p_{1}+q_{2}-n}(-2)^{q_{1}-k-s_{1}+q_{2}-n-s_{2}}\int\limits_{0}^{\infty}\differential{\rho_{2}}V\\!\left(\dfrac{\rho_{2}}{\sqrt{2}}\right)\rho_{2}\exp(-\dfrac{\rho_{2}^{2}}{6\lambda^{2}})\left(\dfrac{\rho_{2}^{2}}{4\lambda^{2}}\right)^{p_{2}-t_{2}+q_{1}-k-s_{1}}=$
$\displaystyle=2\pi\delta_{p_{2}+q_{1}-k,p_{1}+q_{2}-n}(-2)^{q_{1}-k-s_{1}+q_{2}-n-s_{2}}\lambda^{2}\mathfrak{I}_{eh}(p_{2}-t_{2}+q_{1}-k-s_{1}),$
(102)

where we have introduced

$\displaystyle\mathfrak{I}_{eh}(\alpha)=\dfrac{1}{\lambda^{2}}\int\limits_{0}^{\infty}V\left(\dfrac{\rho_{2}}{\sqrt{2}}\right)\exp(-\dfrac{\rho_{2}^{2}}{6\lambda^{2}})\left(\dfrac{\rho_{2}^{2}}{4\lambda^{2}}\right)^{\alpha}\rho_{2}\differential{\rho_{2}}.$
(103)

In the case of Coulomb potential, this integral yields

$\displaystyle\mathfrak{I}^{\textup{C}}_{eh}(\alpha)=\dfrac{\sqrt{3\pi}}{\varepsilon}\left(\dfrac{3}{4}\right)^{\alpha}\left(2\alpha-1\right)!!.$
(104)

Collecting the contributions for the $\widetilde{\mathbf{r}}$–$\mathbf{r}$
part, for the expression (100), we finally obtain

$(*)=\dfrac{1}{3^{2+p_{2}+p_{1}+q_{1}+q_{2}-k-n}}\sum\limits_{t_{2}=0}^{p_{2}}\sum\limits_{t_{1}=0}^{p_{1}}\sum\limits_{s_{1}=0}^{q_{1}-k}\sum\limits_{s_{2}=0}^{q_{2}-n}C^{t_{2}}_{p_{2}}C^{t_{1}}_{p_{1}}C^{s_{1}}_{q_{1}-k}C^{s_{2}}_{q_{2}-n}12\pi\lambda^{2}3^{t_{1}+s_{2}}(t_{1}+s_{2})!\delta_{t_{2}+s_{1},t_{1}+s_{2}}\\\
\times
2\pi\delta_{p_{2}+q_{1}-k,p_{1}+q_{2}-n}(-2)^{q_{1}-k-s_{1}+q_{2}-n-s_{2}}\lambda^{2}\mathfrak{I}^{C}_{eh}(p_{2}-t_{2}+q_{1}-k-s_{1}).$
(105)

Thus, we find the final expression for
$I^{(\textup{T},eh)}_{p_{1}p_{2}q_{1}q_{2}r_{1}r_{2}}$. Continuing the
equation (98), we come to the following expression:

$\displaystyle I^{(\textup{T},eh)}_{p_{1}p_{2}q_{1}q_{2}r_{1}r_{2}}=$
$\displaystyle-\delta_{p_{2}+q_{1}+r_{1},p_{1}+q_{2}+r_{2}}\dfrac{2e^{2}}{3\lambda}\dfrac{(1+(-1)^{p_{1}+p_{2}})}{2^{(q_{1}+q_{2})/2}}\sum\limits_{k=0}^{q_{1}}\sum\limits_{n=0}^{q_{2}}C^{k}_{q_{1}}C^{n}_{q_{2}}\dfrac{\delta_{r_{1}+k,r_{2}+n}(k+r_{1})!}{3^{p_{1}+p_{2}+q_{1}+q_{2}-k-n}}$
$\displaystyle\times\sum\limits_{t_{2}=0}^{p_{2}}\sum\limits_{t_{1}=0}^{p_{1}}\sum\limits_{s_{1}=0}^{q_{1}-k}\sum\limits_{s_{2}=0}^{q_{2}-n}C^{t_{2}}_{p_{2}}C^{t_{1}}_{p_{1}}C^{s_{1}}_{q_{1}-k}C^{s_{2}}_{q_{2}-n}3^{t_{1}+s_{2}}(t_{1}+s_{2})!\delta_{t_{2}+s_{1},t_{1}+s_{2}}(-2)^{q_{1}-k-s_{1}+q_{2}-n-s_{2}}$
$\displaystyle{}\times\mathfrak{I}_{eh}(p_{2}-t_{2}+q_{1}-k-s_{1}).$ (106)

This is the most general expression allowing one to obtain the matrix element
(96) for any spherically symmetric potential. In the case of the Coulomb one,
it takes the following form:

$\displaystyle
I^{(\textup{T},eh,\textup{C})}_{p_{1}p_{2}q_{1}q_{2}r_{1}r_{2}}=-\delta_{p_{2}+q_{1}+r_{1},p_{1}+q_{2}+r_{2}}\dfrac{2\sqrt{2}E_{0}}{\sqrt{3}}\dfrac{(-1)^{q_{1}+q_{2}}[1+(-1)^{p_{1}+p_{2}}]}{3^{p_{1}+q_{2}}2^{p_{1}+p_{2}+(q_{1}+q_{2})/2}}\sum\limits_{k=0}^{q_{1}}\sum\limits_{n=0}^{q_{2}}C^{k}_{q_{1}}C^{n}_{q_{2}}\delta_{r_{1}+k,r_{2}+n}(k+r_{1})!$
(107)
$\displaystyle\times\sum\limits_{t_{2}=0}^{p_{2}}\sum\limits_{t_{1}=0}^{p_{1}}\sum\limits_{s_{1}=0}^{q_{1}-k}\sum\limits_{s_{2}=0}^{q_{2}-n}C^{t_{2}}_{p_{2}}C^{t_{1}}_{p_{1}}C^{s_{1}}_{q_{1}-k}C^{s_{2}}_{q_{2}-n}\delta_{t_{2}+s_{1},t_{1}+s_{2}}(-1)^{k+s_{1}+n+s_{2}}3^{t_{1}+s_{2}+n}2^{t_{1}+t_{2}}(t_{1}+s_{2})!\\!\left[2(p_{2}\\!-\\!t_{2}\\!+\\!q_{1}\\!-\\!k\\!-\\!s_{1}\\!)\\!-\\!1\\!\right]!!.$

The matrix elements corresponding to the electron-electron interaction is
determined via

$\displaystyle
V^{(ee)}_{n_{1}n_{2}n_{h}mln_{1}^{\prime}n_{2}^{\prime}n_{h}^{\prime}m^{\prime}l^{\prime}}\\!\\!=\\!\\!\langle\varphi^{(\textup{T})}_{n_{1}n_{2}n_{h}ml}\text{\textbar}\hat{V}_{ee}\text{\textbar}\varphi^{(\textup{T})}_{n_{1}^{\prime}n_{2}^{\prime}n_{h}^{\prime}m^{\prime}l^{\prime}}\rangle\\!=\\!\\!\dfrac{\sum\sum
C^{n_{1}n_{2}n_{h}ml*}_{a^{\prime}_{1}a^{\prime}_{2}b^{\prime}_{1}b^{\prime}_{2}c^{\prime}_{1}c^{\prime}_{2}}C^{n_{1}n_{2}n_{h}ml}_{a_{1}a_{2}b_{1}b_{2}c_{1}c_{2}}I^{(\textup{T},ee)}_{a_{1}+a_{2}^{\prime},a_{2}+a_{1}^{\prime},b_{1}+b_{2}^{\prime},b_{2}+b_{1}^{\prime},c_{1}+c_{2}^{\prime},c_{2}+c_{1}^{\prime}}}{\sqrt{N_{n_{1}n_{2}n_{h}ml}}\sqrt{N_{n_{1}^{\prime}n_{2}^{\prime}n_{h}^{\prime}m^{\prime}l^{\prime}}}},$
(108)

where

$I^{(\textup{T},ee)}_{p_{1}p_{2}q_{1}q_{2}r_{1}r_{2}}=-\dfrac{e^{2}}{\lambda}\\!\int\\!\\!d\mathbf{r}\\!\\!\int\\!\\!d\mathbf{R}\\!\\!\int\\!\\!d\mathbf{r}_{h}(\xi^{*})^{p_{1}}(\xi)^{p_{2}}(\xi_{h})^{q_{1}}(\xi_{h}^{*})^{q_{2}}(\xi_{R})^{r_{1}}(\xi_{R}^{*})^{r_{2}}\\\
\times\dfrac{\exp[-2(\xi\xi^{*}+\xi_{R}\xi^{*}_{R}+\xi_{h}\xi^{*}_{h})+\sqrt{2}(\xi_{h}\xi^{*}_{R}+\xi_{h}^{*}\xi_{R})]}{\left(2\pi\right)^{3}\lambda^{6}}V(|\mathbf{r}_{2}-\mathbf{r}_{1}|).$
(109)

One then obtains

$\displaystyle I^{(\textup{T},ee)}_{p_{1}p_{2}q_{1}q_{2}r_{1}r_{2}}$
$\displaystyle=\dfrac{1}{(2\pi)^{3}\lambda^{6}}\dfrac{e^{2}}{\lambda}\int\limits_{0}^{\infty}r\differential{r}\int\limits_{0}^{2\pi}\differential{\varphi}V(\sqrt{2}r)\left(\dfrac{r}{2\lambda}\right)^{p1+p2}\exp(-\dfrac{r^{2}}{2\lambda^{2}})\exp[{i\mkern
1.0mu}\varphi(p_{2}-p_{1})]\int\limits
r_{h}\differential{r_{h}}\left(\dfrac{r_{h}}{2\lambda}\right)^{q_{1}+q_{2}}$
$\displaystyle\times\exp(-\dfrac{r_{h}^{2}}{2\lambda^{2}})\int\limits_{0}^{\infty}\\!\differential{R}R\left(\dfrac{R}{2\lambda}\right)^{r_{1}+r_{2}}\exp(-\dfrac{R^{2}}{2\lambda^{2}})\int\limits_{0}^{2\pi}\differential{\varphi_{R}}\int\limits_{0}^{\infty}\differential{\varphi_{h}}\exp[{i\mkern
1.0mu}\varphi_{R}(r_{1}-r_{2})]\exp[{i\mkern 1.0mu}\varphi_{h}(q_{1}-q_{2})]$
$\displaystyle\times\exp[{i\mkern 1.0mu}u\cos{(\varphi_{R}-\varphi_{h})}],$
(110)

where $u=-{i\mkern 1.0mu}Rr_{h}/(\sqrt{2}\lambda^{2})$. The integration over
the angles yields

$\int\limits_{0}^{2\pi}\differential{\varphi_{r}}\exp[{i\mkern
1.0mu}\varphi_{r}(p_{2}-p_{1})]\int\limits_{0}^{2\pi}\differential{\varphi_{R}}\int\limits_{0}^{\infty}\differential{\varphi_{r_{h}}}\exp[{i\mkern
1.0mu}\varphi_{R}(r_{1}-r_{2})]\exp[{i\mkern
1.0mu}\varphi_{r_{h}}(q_{1}-q_{2})]\exp[{i\mkern
1.0mu}u\cos{(\varphi_{R}-\varphi_{r_{h}})}]\\\
=\delta_{p_{1}p_{2}}\delta_{q_{1}-q_{2},r_{1}-r_{2}}{i\mkern
1.0mu}^{r_{2}-r_{1}}(2\pi)^{3}J_{r_{2}-r_{1}}(u).$ (111)

Expanding the Bessel function in a power series and integrating over $R$ and
$r_{h}$, we obtain the final expression for an arbitrary spherically symmetric
potential:

$\displaystyle
I^{(\textup{T},ee)}_{p_{1}p_{2}q_{1}q_{2}r_{1}r_{2}}=\dfrac{e^{2}}{\lambda}\dfrac{1}{2^{q_{2}+(r_{2}+r_{1})/2}}\delta_{p_{1}p_{2}}\delta_{r_{2}-r_{1},q_{1}-q_{2}}\mathfrak{I}_{ee}(2p_{1})\sum\limits_{m=0}^{\infty}\dfrac{1}{2^{m}m!}\dfrac{(m+q_{1})!(m+r_{2})!}{(m+r_{2}-r_{1})!},$
(112)

where we define the $r$ part via

$\displaystyle\mathfrak{I}_{ee}(\alpha)=\dfrac{1}{\lambda^{2}}\int\limits_{0}^{\infty}r\differential{r}V(\sqrt{2}r)\left(\dfrac{r}{2\lambda}\right)^{\alpha}\exp(-\dfrac{r^{2}}{2\lambda^{2}}).$
(113)

In the case of the Coulomb potential, this integral is given by

$\displaystyle\mathfrak{I}_{ee}^{\textup{C}}(\alpha)=\dfrac{\sqrt{\pi}}{2\varepsilon}\dfrac{\left(\alpha-1\right)!!}{2^{\alpha}}.$
(114)

Thus, for the Coulomb potential, the master integral of the electron-electron
interaction reads

$\displaystyle
I^{(\textup{T},ee,\textup{C})}_{p_{1}p_{2}q_{1}q_{2}r_{1}r_{2}}=\dfrac{E_{0}}{\sqrt{2}}\dfrac{\left(2p_{1}-1\right)!!}{2^{q_{2}+(r_{2}+r_{1})/2+2p_{1}}}\delta_{p_{1}p_{2}}\delta_{r_{2}-r_{1},q_{1}-q_{2}}\sum\limits_{m=0}^{\infty}\dfrac{1}{2^{m}m!}\dfrac{(m+q_{1})!(m+r_{2})!}{(m+r_{2}-r_{1})!}.$
(115)

Now we have all of the necessary ingredients for the Hamiltonian matrix
elements (55).

##  Appendix D Exciton and a free electron $(\textup{X},e)$

When solving the three-particle problem discussed in Appendix C, one obtains
discrete energy levels corresponding to bound trion states and a continuous
spectrum. The lower boundary of the latter indicates the energy at which the
trion disintegrates and turns into a combination of an exciton and a free
electron. Since these two constituents do not interact with each other, the
free electron contributes to the total energy according to Eq. (20), while the
exciton contribution can be extracted from the matrix problem (30). Let the
free-electron LL be characterized by quantum number $n_{0}$. The Hamiltonian
is diagonal with respect to $n_{0}$ because the Coulomb interaction relates
solely to the exciton part. The matrix has the following form:

$H^{(\textup{X},e)}_{\mathbf{K},n_{0},n^{\prime}m^{\prime};nm}=E^{(\textup{X},0)}_{nm}\delta_{nn^{\prime}}\delta_{mm^{\prime}}\\\
+V_{\mathbf{K},n^{\prime}-n_{0},m^{\prime};n-n_{0},m}.$ (116)

Here the second term is present only for
$n_{0}\leqslant\mathrm{min}\,(n,\,n^{\prime})$. To obtain the lower energy
level of the continuous spectrum, one also has to minimize the resulting
energy by varying $|\mathbf{K}|$. If the Coulomb mixing among different LLs of
the three-particle problem is neglected, then for given $n_{e}\equiv
n_{1}+n_{2}$ and $n_{h}$ in Eq. (55) one should choose $n=n^{\prime}=n_{e}$
and $m=m^{\prime}=n_{h}$ in Eq. (116). In order to fully take into account the
Coulomb interaction, it is necessary to vary $n$, $n^{\prime}$, $m$,
$m^{\prime}$ in Eq. (116) from zero to sufficiently large numbers, so the
energy eigenvalues converge.

##  Appendix E Examples of calculating the continuum onset in the simplest
cases

Let us discuss here the quantitative features of the energy spectrum of the
third panel in Fig. 1. First, we observe the regions of a discrete spectrum
and continuous bands. The discrete levels correspond to the trion states,
while the lower boundary of each continuous region (blue, violet, orange,
pink) relates to the onset of the exciton continuum. For instance, in the case
$n_{h}=0$ and $n_{e}=1$, there are no discrete levels. To find the lower bound
of each continuous region, we consider the problem involving non-interacting
exciton and electron (see Appendix D). Let us illustrate this approach for
$(n_{e},n_{h})=(0,0)$ and $(0,1)$. In these simplest cases all calculations
are reduced to the analytical formula instead of diagonalization of matrix.
The lower border $E^{\textup{X}}_{00}$ for the region $(0,0)$ can be found by
means of the following relation:

$\displaystyle E^{\textup{X}}_{00}$
$\displaystyle=E^{(\textup{X},0)}_{00}+V_{\mathbf{K}\rightarrow\mathbf{K}_{\textup{min}},00;00}$
$\displaystyle=\frac{1}{2}\hbar\omega_{0}+\left.I^{(\textup{X},eh,\textup{C})}_{00}\right\rvert_{r_{0}\rightarrow
r_{\textup{min}}},$ (117)

where $E^{(\textup{X},0)}_{00}$ is the corresponding free LL [Eq. (40)],
$V_{\mathbf{K},00;00}$ is the Coulomb matrix element [Eq. (52)],
$\omega_{0}=\omega_{e}+\omega_{h}$, and $I^{(\textup{X},eh,\textup{C})}_{00}$
is defined in Eq. (54). As was mentioned above, since an exciton is a neutral
system, the translational symmetry gives rise to a conservation law for vector
$\mathbf{K}$. Then we introduce
$\mathbf{r}_{0}=-c(\mathbf{B}\times\mathbf{K})/(eB^{2})$ and vary the
magnitude of $\mathbf{r}_{0}$ because the energy is independent of its
direction (see Appendix B). In Eq. (117), the energy as a function of $r_{0}$
is minimal (see Fig. 5). For $(n_{e},n_{h})=(0,0)$, one finds
$r_{\textup{min}}=0$. In the case $(1,0)$, we have the following relation:

$\displaystyle E^{\textup{X}}_{01}$
$\displaystyle=E^{(\textup{X},0)}_{01}+V_{\mathbf{K}\rightarrow\mathbf{K}_{\textup{min}},01;01}$
$\displaystyle=\frac{1}{2}\big{(}\hbar\omega_{e}+3\hbar\omega_{h}\big{)}+\left.2I^{(\textup{X},eh,\textup{C})}_{11}\right\rvert_{r_{0}\rightarrow
r_{\textup{min}}},$ (118)

where factor $2$ in the last term comes from the coefficient $B^{01}_{10}$
introduced in Eq. (46). In Fig. 5 we display the functions
$I^{(\textup{X},eh,\textup{C})}_{00}$, and
$2I^{(\textup{X},eh,\textup{C})}_{11}$ versus the dimensionless parameter
$r_{0}/\lambda$ in units of $E_{0}=\sqrt{\pi/2}(e^{2}/\varepsilon\lambda)$,
which represents a characteristic energy scale of all of the matrix elements.
Since the Coulomb contributions to the exciton energy can be continuously
varied by changing the parameter $r_{0}$, the exciton spectrum contains
continuous bands, which are clearly seen in the third panel of Fig. 1. For
GaAs at $B=30~{}\textup{T}$, we obtain
$E^{\textup{X}}_{00}=28.665~{}\textup{meV}$ and
$E^{\textup{X}}_{01}=48.207~{}\textup{meV}$.

Figure 5:  Functions
$\tilde{I}^{(\textup{X},eh,\textup{C})}_{00}=I^{(\textup{X},eh,\textup{C})}_{00}$
and
$\tilde{I}^{(\textup{X},eh,\textup{C})}_{11}=2I^{(\textup{X},eh,\textup{C})}_{11}$
governing the Coulomb matrix elements and defined in Eq. (54) versus
$r_{0}=|\mathbf{r}_{0}|$ in units of $\lambda$. These quantities entering Eqs.
(117) and (118) determine the exciton energy. The vector $\mathbf{r}_{0}$ is
defined via $\mathbf{r}_{0}=-c(\mathbf{B}\times\mathbf{K})/(eB^{2})$, where
$\mathbf{K}$ is a conserved vector arising due to the translational symmetry
of the system. Such a continuous dependence on the parameter $r_{0}$ leads to
a band structure of the exciton spectrum.

##  Appendix F Rules for filling the trion Hamiltonian

Figure 6: Detailed spectrum which corresponds to the second (purple) part of
the third vertical panel in Fig. 1 for GaAs and $B=30~{}\textup{T}$. The
points are grouped in accordance with specific values of $Mz$ (changes along
the $x$ axis). The gray points are associated with triplets ($S_{e}=1$), while
purple points correspond to singlets $S_{e}=0$. Recall that the indices
associated with the free LLs are fixed and have the following values:
$n_{e}=0$, $n_{h}=1$. For each integer $M_{z}$ we have $50$ points for each
energy level.  Figure 7: Detailed spectrum which corresponds to the lower part
of the fourth vertical panel in Fig. 1 for GaAs and $B=30~{}\textup{T}$. The
points are grouped in accordance with specific values of $Mz$ (changes along
the $x$ axis). The gray points are associated with triplets ($S_{e}=1$), while
black points correspond to singlets $S_{e}=0$. The indices associated with the
free LLs change as $n_{e}=0\dots 4$, $n_{h}=0\dots 4$. For each integer
$M_{z}$ we have four points for each energy level.

In order to obtain the spectrum of a trion, it is necessary to construct the
Hamiltonian correctly. Let us look at the filling of the Hamiltonian using the
example of the Fig. 1. First, we analyze the third panel, where the matrix
(55) is diagonalized neglecting the matrix elements between different LLs.
Here, we employ the method proposed in Refs. Dzyubenko and Sivachenko (2000);
Dzyubenko (2000, 2001). In practice, the matrix is filled in such a way that
we do not get negative values of the indices, and also keep in mind that the
parity and the projection of the total orbital angular momentum are preserved.
Thus, if the arbitrary matrix element is given by
$H^{(\textup{T})}_{n_{1}n_{2}n_{h}ml;n_{1}^{\prime}n_{2}^{\prime}n_{h}^{\prime}m^{\prime}l^{\prime}}$,
then the indices should satisfy the following conditions: $n_{2}=n_{e}-n_{1}$,
$n_{2}^{\prime}=n_{e}^{\prime}-n_{1}^{\prime}$,
$m=\tilde{m}+(1+(-1)^{n_{1}+S_{e}})/2$,
$m^{\prime}=\tilde{m}^{\prime}+(1+(-1)^{n_{1}^{\prime}+S_{e}})/2$,
$l=\tilde{m}+(1+(-1)^{n_{1}+S_{e}})/2+M_{z}-n_{e}+n_{h}$,
$l^{\prime}=\tilde{m}^{\prime}+(1+(-1)^{n_{1}^{\prime}+S_{e}})/2+M_{z}-n_{e}^{\prime}+n_{h}^{\prime}$,
where $n_{e}$ and $n_{e}^{\prime}$ change in the range
$[n_{e,\textup{min}},n_{e,\textup{max}}]$ with step 1, $n_{1}$ and
$n_{1}^{\prime}$ change in the range $[0,n_{e,\textup{max}}]$, $n_{h}$ and
$n_{h}^{\prime}$ change in the range $[n_{h,\textup{min}},n_{h,\textup{max}}]$
with step 1, $\tilde{m}$ and $\tilde{m}^{\prime}$ change in the range $[0,M]$
with step 2. Moreover, the following additional requirements should be
verified: $n_{e}-n_{1}\geq 0$, $n_{e}^{\prime}-n_{1}^{\prime}\geq 0$,
$\tilde{m}+(1+(-1)^{n_{1}+S_{e}})/2+M_{z}-n_{e}+n_{h}\geq 0$,
$\tilde{m}^{\prime}+(1+(-1)^{n_{1}^{\prime}+S_{e}})/2+M_{z}-n_{e}^{\prime}+n_{h}^{\prime}\geq
0$. The parameter $M$ is responsible for the convergence of the results
affecting the size of the model Hilbert space.

Thus, we should choose the following parameters: the value of the projection
of the angular momentum $M_{z}$, the spin quantum number $S_{e}=0$ or
$S_{e}=1$, $M$ (for the calculations at third panel in Fig. 1 it is up to
$90$), and the range of the parameters associated with the free LLs
($n_{h,\textup{min}}$, $n_{h,\textup{max}}$, $n_{e,\textup{min}}$,
$n_{e,\textup{max}}$). For instance, in order to obtain the blue region of the
spectra, we fix
$n_{h,\textup{min}}=n_{h,\textup{max}}=n_{e,\textup{min}}=n_{e,\textup{max}}=0$,
$M=90$ and collect the eigenvalues for different $M_{z}\in\\{-7,\dots,0\\}$
and $S_{e}=0$, $1$. The same steps we do for other colored parts presented in
the third vertical panel in Fig. 1, the difference is only in the choice of
$n_{h,\textup{min}}$, $n_{h,\textup{max}}$, $n_{e,\textup{min}}$,
$n_{e,\textup{max}}$.

If the Coulomb mixing is taken into account, the corresponding spectrum is
presented in the fourth vertical panel in Fig. 1. Here, we chose
$n_{h,\textup{min}}=n_{e,\textup{min}}=0$, and
$n_{h,\textup{max}}=n_{e,\textup{max}}=4$, while $M$ is varied from $0$ to
$12$. The smaller values of $M$ are employed due to the larger dimension of
the model Hilbert space and corresponding computational complexity. Note that
despite such a reduction in the dimension, this turns out to be sufficient for
the convergence of the lower discrete levels.

To present a more detailed spectrum structure, we provided two examples of
spectra where different parts corresponding to specific values of $M_{z}$ and
$S_{e}$ are distinguished. One of them (Fig. 6) relates to the purple part of
the third vertical panel in Fig. 1, while the second example (Fig. 7)
corresponds to the lower part of the fourth vertical panel in the same figure,
where we take into account the Coulomb-induced mixing of the different LLs.
The corresponding details can be found in Appendix G.

##  Appendix G Detailing of spectrum structure

Here we would like to illuminate the technicalities of our spectrum
calculations. Of course, in practical trion computations, the regions that we
associate with the exciton continuum are not continuous since the basis set
employed is always finite, if huge. However, in order to correctly determine
the continuous zones of the spectrum and establish the true discrete levels
that we associate with the trion states, it is sufficient to (i) follow the
behavior of the levels as the basis increases with increasing parameter $M$,
(ii) carry out independent calculations of the lower exciton level in order to
verify the lower border of the exciton continuum. Let us consider the first
example in Fig. 6. In order to construct it, we resolved the details of the
purple sector of the third vertical panel in Fig. 1 with respect to the values
of the projection of the total angular momentum $M_{z}$. The points are not
merged into a single area as was done previously (in Fig. 1 we assigned a
horizontal line of the same width to each point).

Figure 8: Dependence of the lowest trion level (gray circles) and lowest
exciton level (gray stars) on the maximal value of the LLs taken into account
when constructing the model Hilbert space in the case of GaAs at
$B=10~{}\textup{T}$. Here we assume
$n_{e}^{\textup{max}}=n_{h}^{\textup{max}}=n^{\textup{max}}$.

In Fig. 6, we indicate the corresponding values of the projection $M_{z}$
below each group of the points at the horizontal axis. Blurring is due to the
fact that when moving along the $x$ axis in the vicinity of each value of
$M_{z}$, we change the value of the parameter $M$ demonstrating (i) the
emergence of new levels, (ii) how the existing levels change their positions.
Moreover, levels that have different spin number $S_{e}$ are marked by
different colors (gray and purple). The spectrum shows two qualitatively
different types of behavior. Above a certain level, with increasing parameter
$M$, we observe a significant change in levels and the appearance of more and
more new ones. However, in the region below the mentioned border, the levels
practically do not change from the very beginning and stand apart from each
other. Furthermore, new levels do not appear there. As one might guess, the
upper part is associated with the exciton continuum, while the lower part
yields the trion levels. Let us also note that the upper region tends to fill
the entire part of the spectrum up to the corresponding free LL, which is also
marked by the upper horizontal line. As was already mentioned, to confirm this
qualitative separation, we calculated the lower exciton level within the
analysis of noninteracting exciton and electron (see Appendix D). In Fig. 6
this level is labeled by $E^{\textup{X}}_{01}$. In Fig. 7 we show another
example of detailing the projection of the total angular momentum $M_{z}$ and
$S_{e}$ in the case of a more complicated consideration when mixing of various
free LLs is taken into account. In order to obtain a clear structure of the
discrete levels (trion ones), we limit ourselves only by the lower part of the
spectrum presented in Fig. 1 within the fourth vertical panel.

Figure 9: Behavior of the lowest exciton level with growth of Hilbert space
(gray dots) for GaAs at $B=10~{}\textup{T}$. By the increasing we mean the
number of LLs $n$ taken into account. These points are fitted using different
extrapolations, built on the basis of the data of different length. The
horizontal (dashed-dotted) lines indicate the values of corresponding
$a_{\textup{X}}(n_{\textup{max}})$, for $n_{\textup{max}}=3,5,7,9$.

It is evident that with the inclusion of higher values of the free LLs $n_{e}$
and $n_{h}$, the dimension of the model Hilbert space grows significantly. For
this reason, we did not achieve such values of $M$ as in the previous case.
However, as it turned out, this was not necessary. Moreover, it is easier to
identify the discrete levels since they demonstrate stability with respect to
increasing $M$ from the very beginning. An example of such stability is
demonstrated for the two lower trion levels. Moreover, for the lowest one, we
display an inset in Fig. 7 which clearly shows the converged value of the
energy level.

##  Appendix H Convergence of energy levels and extrapolation procedure

Our predictions should converge with respect to the number of the LLs
incorporated within the model Hilbert space. It turns out that even within the
setup that takes into account all levels up to the fifth ones
($n_{e}=0,\dots,5$ and $n_{h}=0,\dots,5$), we cannot conclude that the
convergence has been achieved. In order to overcome this problem, we utilize
an extrapolation procedure described below. In Fig. 8, we demonstrate the
energy of one of the lowest trion levels and the lowest exciton level as
functions of the number of LLs included within the model Hilbert space in the
case of GaAs at $B=10~{}\textup{T}$. We depict first the energies calculated
in the zero approximation ($n^{\textup{max}}_{e}=n^{\textup{max}}_{h}=0$, blue
dots), as was done in Ref. Dzyubenko and Sivachenko (2000), and then we
present the results computed via our general technique: for the trion problem
we reached $n^{\textup{max}}_{e}=n^{\textup{max}}_{h}=5$, while for the
exciton problem we took into acoount the states up to
$n^{\textup{max}}_{e}=n^{\textup{max}}_{h}=9$. One observes that the
convergence is very slow and, based on the existing number of points, it is
difficult to precisely determine the values of the true trion and exciton
levels. However, using these points, we performed an extrapolation procedure
allowing us to obtain more accurate predictions. In Fig. 8, they are depicted
by the horizontal lines. We note that although extrapolating the data can
never yield absolutely accurate values, it is not possible to extract the
energy levels without extrapolation since the convergence is unlikely to be
achieved by a brute-force approach, i.e, by simply increasing the dimension of
the model Hilbert space via including additional LLs, due to the slow
convergence rate and computational limitations.

Thus, to overcome this problem, we extrapolate the energy values in order to
determine more accurate estimates. This procedure is obviously not
unambiguous, so the predictions will always possess, in fact, some
uncertainties. Here we will describe the approach we used for Fig. 12. Below,
we consider one point where $B=10~{}\textup{T}$.

Figure 10: Gray dots show the behavior of the free constant
$a_{\textup{X}}(n)$ entering extrapolation formula (119) with the growth of
Hilbert space (increasing the number of taken into account LLs). The dotted
green line presents the highest extrapolation $f_{a_{\textup{X}}}(9,n)$
obtained on the basis of all gray dots, while the red dashed line is the
asymptotic value of the corresponding free constant which is interpreted by us
as the final estimate for $E^{X}_{\infty}$.  Figure 11: Gray dots demonstrate
the behavior of one of the discrete trion levels with growth of Hilbert space
(increasing of $n_{\textup{max}}$). These points are fitted using different
extrapolations, built on the basis of the data of different length. The
horizontal (dashed-dotted) lines indicate the values of corresponding
$a_{\textup{T}}(n_{\textup{max}})$, for $n_{\textup{max}}=2,3,4$. The value of
$a_{\textup{T}}(5)$ is identified separately by solid red line and considered
by us as the final estimate of $E^{\textup{T,}1}_{\infty}$ at fixed $B$.
Figure 12: The dependence of the energy spectrum of the system as function of
the magnetic field B. Solid red lines correspond to the two discrete trion
states, that we get in our calculation the dashed red line corresponds to the
lower boundary of the energy continuum solid blue line corresponds to the
trion bound state reported in Ref. Dzyubenko (2000), and dashed blue line to
corresponding energetic continuum boundary. Solid black line shows the
position of the lowest Landau level corresponding to non-interacting
particles. The increase of magnetic field does not lead to the merging of red
and blue lines, which means, that the role of the Coulomb mixing between
different Landau levels does not decrease with magnetic field.

The basis size is governed by the number $n_{\textup{max}}$, which is a
maximal value of the LL quantum numbers $n_{e}$ and $n_{h}$. Let us first
analyze the behavior of the lowest exciton level. The change in its position
as a function of $n_{\textup{max}}$ is shown in Fig. 9 (see gray dots). An
important issue here is the choice of the form of the approximation function.
Our analysis of the data obtained suggests that the natural approximation
reads

$\displaystyle\\!\\!f_{\textup{X}}(n_{\textup{max}},n)\\!=\\!a_{\textup{X}}(n_{\textup{max}})+\dfrac{b_{\textup{X}}(n_{\textup{max}})}{n^{k_{\textup{X}}(n_{\textup{max}})}+c_{\textup{X}}(n_{\textup{max}})}.$
(119)

It is clear that the fitting constants should alter with the appearance of new
data points regarding $n_{\textup{max}}$. If these constants converge
relatively well, we could conclude that the approximation formula is chosen
correctly. Here we are only interested in the parameter $a_{\textup{X}}$ and
its asymptotic value $a_{\textup{X}}(n\rightarrow\infty)$ which gives us the
final estimate for the lowest exciton level.

In Fig. 10, we present values of $a_{\textup{X}}(n_{\textup{max}})$ involved
in the different extrapolation formulas $f_{\textup{X}}(n_{\textup{max}},n)$.
Since the dimension of the model Hilbert space increases, the value of the
constant itself also changes but the shifts are observed in a convergent
manner. We construct a similar extrapolation for
$a_{\textup{X}}(n_{\textup{max}})$ which we plot in the figure. The
corresponding extrapolation function is labeled by $f_{a_{\textup{X}}}$. In
the figure we also depict the red line which gives the asymptotic value
$a_{\textup{X}}(n\rightarrow\infty)$ and is considered as a final estimate for
$E^{\textup{X}}_{\infty}$. The same steps were applied for extracting
$E^{\textup{X}}_{\infty}$ for other values of the magnetic field $B$ presented
in Fig. 12.

Finally, let us analyze the behavior of one of the lowest trion level. The
dependence of its position as a function of $n_{\textup{max}}$ is demonstrated
in Fig. 11. In Table 1, we also provide the numbers we will extrapolate.

Table 1: Values of one of the lowest trion levels calculated at $B=10~{}\textup{T}$ for GaAs. The numbers correspond to different values of the LLs quantum number $n$. The corresponding values of the conserved quantities $S_{e}$ and $M_{z}$ are also shown. $n$ | $E^{\textup{T,}1}_{nn}\textup{(meV)}$ | $M_{z}$ | $S_{e}$ | $n$ | $E^{\textup{T,}1}_{nn}\textup{(meV)}$ | $M_{z}$ | $S_{e}$
---|---|---|---|---|---|---|---
$0$ | $1.5172$ | $-1$ | $1$ | $3$ | $-4.8842$ | $0$ | $0$
$1$ | $-1.0890$ | $0$ | $0$ | $4$ | $-5.9807$ | $0$ | $0$
$2$ | $-3.3553$ | $0$ | $0$ | $5$ | $-6.8054$ | $0$ | $0$

The form of the approximation function is chosen the same way as in the
exciton case:

$\displaystyle\\!\\!f_{\textup{T},1}(n_{\textup{max}},n)\\!=\\!a_{\textup{T}}(n_{\textup{max}})+\dfrac{b_{\textup{T}}(n_{\textup{max}})}{n^{k_{\textup{T}}(n_{\textup{max}})}+c_{\textup{T}}(n_{\textup{max}})}.$
(120)

As in the case of the lowest exciton level, it is necessary to monitor the
convergence of the extrapolation constants when $n_{\textup{max}}$ increases.
If they converge relatively well, we can conclude that the approximation
formula is chosen correctly. We are only interested in the parameter
$a_{\textup{T}}$. For the approximations that are available to us, we extract
the following values: $a_{\textup{T}}(2)=-7.9902~{}\textup{meV}$,
$a_{\textup{T}}(3)=-9.8591~{}\textup{meV}$,
$a_{\textup{T}}(4)=-10.400~{}\textup{meV}$,
$a_{\textup{T}}(5)=-10.783~{}\textup{meV}$. Unfortunately, in contrast to the
exciton problem, such a number of points is not enough to construct an
extrapolation for the behavior of the constant itself; here we decided to
limit ourselves to only its upper estimate and consider
$a_{\textup{T}}(5)=-10.783~{}\textup{meV}$ as the final estimate for the
lowest trion level $E^{\textup{T},1}_{\infty}$ for $B=10~{}\textup{T}$. The
same procedure is repeated for all points in Fig. 12

Finally, let us discuss the intersection between the two trion levels revealed
in Fig. 12. First, we note that they may indeed intersect since they have
different values of the spin quantum number $S_{e}$. Second, we underline that
the exact treatment of the three-particle problem is very complicated from the
computational viewpoint, so the quantitative predictions are, in fact, not
sufficiently accurate to properly resolve the relative position of the two
levels. Although one can evidently observe both of them, it might be difficult
to identify the corresponding energy difference for a given value of the
magnetic field strength and, therefore, to accurately localize the
intersection point. Providing a high-resolution picture represents a
formidable task and an important challenge for future studies.

##  Appendix I Mathematica subroutine

As Supplemental Material, we attach Mathematica files with the help of which
one can calculate any matrix elements and the spectra presented in the paper.
The archive contains two files: func_spec.m which includes basic functions and
a Mathematica notebook Trion_exciton_spectrum_in_megnetic_field.nb which
contains examples of calculations. To speed up the computations, we calculated
an extensive database of the matrix elements for exciton and trion problems,
which are located in directories data_ex and data_tr.

## References

  * Miura (2008) N. Miura,  _Physics of Semiconductors in High Magnetic Fields_ (Oxford University Press Inc., 2008).
  * Arora (2021) A. Arora, “Magneto-optics of layered two-dimensional semiconductors and heterostructures: Progress and prospects,” J. Appl. Phys 129 (2021).
  * Landau (1930) L. Landau, “Diamagnetismus der metalle,” Z. Phys. 64, 629 (1930).
  * Akimoto and Hasegawa (1967) O. Akimoto and H. Hasegawa, “Interband optical transitions in extremely anisotropic semiconductors. II. coexistence of exciton and the landau levels,” J. Phys. Soc. Jpn 22, 181 (1967).
  * Edelstein  _et al._ (1989) W. Edelstein, H. N. Spector, and R. Marasas, “Two-dimensional excitons in magnetic fields,” Phys. Rev. B 39, 7697 (1989).
  * MacDonald and Ritchie (1986) A. H. MacDonald and D. S. Ritchie, “Hydrogenic energy levels in two dimensions at arbitrary magnetic fields,” Phys. Rev. B 33, 8336 (1986).
  * Potemski  _et al._ (1991) M. Potemski, L. Via, G. E. W. Bauer, J. C. Maan, K. Ploog, and G. Weimann, “Magnetoexcitons in narrow gaas/${\mathrm{ga}}_{1\mathrm{-}\mathit{x}}$${\mathrm{al}}_{\mathit{x}}$as quantum wells,” Phys. Rev. B 43, 14707 (1991).
  * Yang and Sham (1987) S.-R. E. Yang and L. J. Sham, “Theory of magnetoexcitons in quantum wells,” Phys. Rev. Lett. 58, 2598 (1987).
  * Bodnar  _et al._ (2017) S. Yu. Bodnar, P. S. Grigoryev, D. K. Loginov, V. G. Davydov, Yu. P. Efimov, S. A. Eliseev, V. A. Lovtcius, E. V. Ubyivovk, V. Yu. Mikhailovskii, and I. V. Ignatiev, “Exciton mass increase in a gaas/algaas quantum well in a transverse magnetic field,” Phys. Rev. B 95, 195311 (2017).
  * Gor’kov and Dzyaloshinskii (1968) L. P. Gor’kov and I. E. Dzyaloshinskii, “Contribution to the theory of the mott exciton in a strong magnetic field,” JETP 26, 449 (1968).
  * Avron and Herbst (1978) J. E. Avron and I. W. Herbst, “Separation of center of mass in homogeneous magnetic fields,” Ann. Phys. (N.Y.) 114, 114 (1978).
  * Lerner and Lozovik (1980) I. V. Lerner and Yu. E. Lozovik, “Mott exciton in a quasi-two-dimensional semiconductor in a strong magnetic field,” JETP 51, 588 (1980), [Zh. Eksp. Teor. Fiz. 78, 1167 (1980)].
  * He  _et al._ (2016) Y.-M. He, O. Iff, N. Lundt, V. Baumann, M. Davanco, K. Srinivasan, S. Höfling, and C. Schneider, “Identifying multiexcitons in $\mathrm{Mo}{\mathrm{s}}_{2}$ monolayers at room temperature,” Nat. Comm. 7, 13409 (2016).
  * Stevens  _et al._ (2018) C. E. Stevens, J. Paul, T. Cox, P. K. Sahoo, Humberto R. Gutiérrez, V. Turkowski, D. Semenov, S. A. McGill, M. D. Kapetanakis, I. E. Perakis, D. J. Hilton, and D. Karaiskaj, “Biexcitons in monolayer transition metal dichalcogenides tuned by magnetic fields,” Nat. Comm. 9, 3720 (2018).
  * Torche and Bester (2021) A. Torche and G. Bester, “Biexcitons fine structure and non-equilibrium effects in transition metal dichalcogenides monolayers from first principles,” Comm. Phys. 4, 67 (2021).
  * Finkelstein and Bar-Joseph (1995) G. Finkelstein and I. Bar-Joseph, “Charged excitons in GaAs quantum wells,” Il Nuovo Cimento D 17, 1239 (1995).
  * Shields  _et al._ (1995) A. J. Shields, M. Pepper, M. Y. Simmons, and D. A. Ritchie, “Spin-triplet negatively charged excitons in gaas quantum wells,” Phys. Rev. B 52, 7841 (1995).
  * Sanvitto  _et al._ (2002) D. Sanvitto, D. M. Whittaker, A. J. Shields, M. Y. Simmons, D. A. Ritchie, and M. Pepper, “Origin of the oscillator strength of the triplet state of a trion in a magnetic field,” Phys. Rev. Lett. 89, 246805 (2002).
  * Smoleński  _et al._ (2019) T. Smoleński, O. Cotlet, A. Popert, P. Back, Y. Shimazaki, P. Knüppel, N. Dietler, T. Taniguchi, K. Watanabe, M. Kroner, and A. Imamoglu, “Interaction-induced shubnikov–de haas oscillations in optical conductivity of monolayer ${\mathrm{mose}}_{2}$,” Phys. Rev. Lett. 123, 097403 (2019).
  * Munteanu  _et al._ (2000) F. M. Munteanu, D. G. Rickel, C. H. Perry, Y. Kim, J. A. Simmons, and J. L. Reno, “Bright and dark triplet states of the negatively charged magnetoexcitons revealed in photoluminescence and time-resolved measurements,” Phys. Rev. B 62, 16835 (2000).
  * Andronikov  _et al._ (2005) D. Andronikov, V. Kochereshko, A. Platonov, T. Barrick, S. A. Crooker, and G. Karczewski, “Singlet and triplet trion states in high magnetic fields: Photoluminescence and reflectivity spectra of modulation-doped $\mathrm{Cd}\mathrm{Te}{\mathrm{cd}}_{0.7}{\mathrm{mg}}_{0.3}\mathrm{Te}$ quantum wells,” Phys. Rev. B 72, 165339 (2005).
  * Klein  _et al._ (2022) J. Klein, M. Florian, A. Hötger, A. Steinhoff, A. Delhomme, T. Taniguchi, K. Watanabe, F. Jahnke, A. W. Holleitner, M. Potemski, C. Faugeras, A. V. Stier, and J. J. Finley, “Trions in ${\mathrm{mos}}_{2}$ are quantum superpositions of intra- and intervalley spin states,” Phys. Rev. B 105, L041302 (2022).
  * Van Tuan  _et al._ (2022) D. Van Tuan, S.F. Shi, X. Xu, S. A. Crooker, and H. Dery, “Six-body and eight-body exciton states in monolayer ${\mathrm{wse}}_{2}$,” Phys. Rev. Lett. 129, 076801 (2022).
  * Wang  _et al._ (2020) T. Wang, Z. Li, Z. Lu, Y. Li, S. Miao, Z. Lian, Y. Meng, M. Blei, T. Taniguchi, K. Watanabe, S. Tongay, W. Yao, D. Smirnov, C. Zhang, and S.-F. Shi, “Observation of quantized exciton energies in monolayer ${\mathrm{wse}}_{2}$ under a strong magnetic field,” Phys. Rev. X 10, 021024 (2020).
  * Li  _et al._ (2020) Z. Li, T. Wang, S. Miao, Y. Li, Z. Lu, C. Jin, Z. Lian, Y. Meng, M. Blei, T. Taniguchi, K. Watanabe, S. Tongay, W. Yao, D. Smirnov, C. Zhang, and S.-F. Shi, “Phonon-exciton interactions in WSe2 under a quantizing magnetic field,” Nat. Commun. 11 (2020).
  * Stébé  _et al._ (1987) B. Stébé, E. Feddi, and G. Munschy, “Excitonic trions in a low magnetic field,” Phys. Rev. B 35, 4331 (1987).
  * Chapman  _et al._ (1998) J. R. Chapman, N. F. Johnson, and V. N. Nicopoulos, “Photoluminescence signatures of negatively charged magnetoexcitons,” Phys. Rev. B 57, 1762 (1998).
  * Finkelstein  _et al._ (1995) G. Finkelstein, H. Shtrikman, and I. Bar-Joseph, “Optical spectroscopy of a two-dimensional electron gas near the metal-insulator transition,” Phys. Rev. Lett. 74, 976 (1995).
  * Wojs and Hawrylak (1995) A. Wojs and P. Hawrylak, “Negatively charged magnetoexcitons in quantum dots,” Phys. Rev. B 51, 10880 (1995).
  * Palacios  _et al._ (1996) J. J. Palacios, D. Yoshioka, and A. H. MacDonald, “Long-lived charged multiple-exciton complexes in strong magnetic fields,” Phys. Rev. B 54, R2296 (1996).
  * Whittaker and Shields (1997) D. M. Whittaker and A. J. Shields, “Theory of ${X}^{-}$ at high magnetic fields,” Phys. Rev. B 56, 15185 (1997).
  * Stébé  _et al._ (1998) B. Stébé, E. Feddi, A. Ainane, and F. Dujardin, “Optical and magneto-optical absorption of negatively charged excitons in three- and two-dimensional semiconductors,” Phys. Rev. B 58, 9926 (1998).
  * Kozhevnikov  _et al._ (1997) M. Kozhevnikov, E. Cohen, A. Ron, H. Shtrikman, and L. N. Pfeiffer, “Comparative cyclotron resonance and luminescence study of photoexcited particles in mixed type-i–type-ii gaas/alas multiple quantum wells,” Phys. Rev. B 56, 2044 (1997).
  * Hayne  _et al._ (1999) M. Hayne, C. L. Jones, R. Bogaerts, C. Riva, A. Usher, F. M. Peeters, F. Herlach, V. V. Moshchalkov, and M. Henini, “Photoluminescence of negatively charged excitons in high magnetic fields,” Phys. Rev. B 59, 2927 (1999).
  * Riva  _et al._ (2002) C. Riva, F. M. Peeters, and K. Varga, “Theory of trions in quantum wells,” Phys. E: Low-Dimens. Syst. Nanostructures. 12, 543 (2002).
  * Riva  _et al._ (2001) C. Riva, F. M. Peeters, and K. Varga, “Magnetic field dependence of the energy of negatively charged excitons in semiconductor quantum wells,” Phys. Rev. B 63, 115302 (2001).
  * Kuznetsov  _et al._ (2018) V. A. Kuznetsov, L. V. Kulik, M. D. Velikanov, A. S. Zhuravlev, A. V. Gorbunov, S. Schmult, and I. V. Kukushkin, “Three-particle electron-hole complexes in two-dimensional electron systems,” Phys. Rev. B 98, 205303 (2018).
  * Filikhin  _et al._ (2018) I. Filikhin, R. Ya. Kezerashvili, Sh. M. Tsiklauri, and B. Vlahovic, “Trions in bulk and monolayer materials: Faddeev equations and hyperspherical harmonics,” Nanotechnology 29, 124002 (2018).
  * Chang and Chang (2022) Y.-W. Chang and Y.-C. Chang, “Theory of magnetic-field effect on trions in two-dimensional materials,” J. Chem. Phys. 157 (2022).
  * Chang and Chang (2021) Y.-W. Chang and Y.-C. Chang, “Variationally optimized orbital approach to trions in two-dimensional materials,” J. Chem. Phys. 155 (2021).
  * Stébé and Moradi (2000) B. Stébé and A. Moradi, “Charged excitons in a low magnetic field in $\mathrm{GaAs}/{\mathrm{ga}}_{1-x}{\mathrm{al}}_{x}\mathrm{As}$ and ${\mathrm{c}\mathrm{d}\mathrm{t}\mathrm{e}/\mathrm{c}\mathrm{d}}_{1-x}{\mathrm{zn}}_{x}\mathrm{Te}$ semiconductor quantum wells,” Phys. Rev. B 61, 2888 (2000).
  * Kaur  _et al._ (2019) H. Kaur, S. Sharma, P. Aggarwal, and R. K. Hazra, “Binding energies and chemical potential of neutral and charged exciton-complexes of transition metal dichalcogenides/anisotropic 2-d quantum dots in magnetic field by exact multi-pole expansion of coulomb correlations,” Phys. E: Low-Dimens. Syst. Nanostructures. 108, 347 (2019).
  * Schmelcher and Cederbaum (1991) P. Schmelcher and L. S. Cederbaum, “Interaction of the landau orbitals of atomic ions in a magnetic field with electronic motion,” Phys. Rev. A 43, 287 (1991).
  * Dzyubenko and Sivachenko (2000) A. B. Dzyubenko and A. Yu. Sivachenko, “Charged magnetoexcitons in two-dimensions: Magnetic translations and families of dark states,” Phys. Rev. Lett. 84, 4429 (2000).
  * Dzyubenko (2000) A. B. Dzyubenko, “Two-dimensional charged electron–hole complexes in magnetic fields: keeping magnetic translations preserved,” Solid State Commun. 113, 683 (2000).
  * Dzyubenko (2001) A. B. Dzyubenko, “Charged hydrogenic problem in a magnetic field: Noncommutative translations, unitary transformations, and coherent states,” Phys. Rev. B 65, 035318 (2001).
  * Landau and Lifshitz (1977) L. D. Landau and E. M. Lifshitz,  _Quantum Mechanics: Non-Relativistic Theory_ (Pergamon Press, Oxford, 1977).
  * Kirzhnits (1967) D. A. Kirzhnits,  _Field Theoretical Methods in Many-Body Systems_ (Pergamon Press, Oxford, 1967) p. 372.