text:
# Variational bounds on the ground-state energy of three electrons and one hole in two-dimension

## I Introduction

The research on trions in 2D systems has shown many interesting results in recent years. In a broader sense, the trion problem, not limited to 2D, is linked to the existence and the stability of negative hydrogen ion H- which has not only important implication for atomic physics and astrophysics, but also vital application to other fields such as physical chemistry, accelerator engineering, and controlled fusion.

In a trion, or charged exciton, two electrons are bound to one hole despite the strong Coulomb repulsion between electrons. (Likewise two holes can be bound to one electron as well but the discussion that follows is the same if “hole” is switched with “electron.”) It is quite a marvel that the two electrons would not fall apart from the strong Coulomb repulsion between them. More surprising is that the ground-state energy of such a complex is lower than that of an exciton, at least in the limit of heavy hole (i.e., the mass of the hole is much larger than that of electrons, which effectively results into the H- problem).

First let us try to understand in terms of classical mechanics why trion forms despite the strong Coulomb repulsion.

Suppose that two electrons are on two sides of one hole as shown in Fig. 1. The electrostatic forces felt by electron 1 (e1) are from the hole and electron 2 (e2). Because the forces are governed by the inverse square Coulomb’s law, the repulsion exerted on e1 by e2 is only $`1/4`$ that of the attraction given by the hole. Similar argument applies to e2 and therefore these two electrons tend to come together to form a bound state in spite of the strong repulsion among like charges. In order for electrons not to fall into the hole we can assume in a classical fashion that the electrons rotate around the hole in such a way that the inward electrostatic force is balanced by the centrifugal force. In quantum mechanics such tendency of electron’s coming together will be stopped by the penalty in kinetic energy which increases as $`a^2`$, where $`a`$ is roughly the average distance between the hole and one of the electrons. Note that the gain in Coulomb energy varies at a slower pace approximately as $`a^1`$ so that it would be hard for the electrons to come near the hole without limit, if the possible electron-hole recombination is not taken into account. (The trion will not be stable if the recombination is spontaneous. We will discuss this problem further in Sec. VI.)

Therefore the bound state of trion is not quite a surprise after all. It basically results from the inverse square Coulomb’s law and a favorable geometric configuration of two electrons and one hole.

## II Three electrons and one hole

Now we want to take one step further to see if a similar physical picture generalizes to the case of three electrons and one hole. In the heavy hole limit, this problem is in effect that for the doubly charged negative hydrogen ion H. There is a much celebrated theorem by E. H. Lieb saying that H is not stable against dissociation into H- and e- in vacuum. Here vacuum means that there is not any other charged particle nearby at finite distance except for the particles forming the H ion. We need a few more words on this theorem. In Lieb’s original papers, a bound system is defined as an eigenstate of the Hamiltonian. This is a strong condition, which leaves room for metastable states for H. The existence of possible metastable states for H is currently a highly debated issue. Lieb’s theorem either does not rule out any possible H state whose total energy is negative. Here the energy zero is taken when particles are infinitely separated and thus with arbitrarily small kinetic and potential energy.

In this paper we study the bound on the ground-state energy of three electrons and one hole at close distance in 2D. Instead of having a heavy positive charge, we assume that the masses of the four are all the same. Our result shows strictly that the ground-state energy is negative when the four particles are relatively localized one to another.

Before starting quantum mechanical calculation, we want to show in a classical picture why such a state of negative energy might be possible. In Fig. 2, three electrons are at the vertices of an equilateral triangle, with one hole in the center. According to the simple knowledge of electrostatics, the force on the central hole is 0 while those on the electrons are $`\frac{(\sqrt{3}1)e^2}{4\sqrt{3}\pi ϵa^2}`$ inwards to the center, where $`a`$ is the distance between the hole and one electron. Not to make themselves fall into the hole the electrons may rotate around the hole in such a way that the inward electrostatic force is balanced by the centrifugal force, just as in a trion. Intuitively a four-body system with a net charge of $`2e`$ may not be stable because of the strong Coulomb repulsion. However under the special geometry as shown, the peculiarity of Coulomb’s law (an inverse square law) actually favors a bound state of 3 electrons and 1 hole (3E+1H), at least within the framework of classical electrodynamics. Thus these electrons tend to move to the hole, giving rise to an effective binding among them.

Note that the total charge of such a complex is $`2e`$ which just reminds us of the Cooper pair that brings in superconductivity, because the complex of four fermions (three electrons and one hole) is now effectively a charged boson.

The stability of such complex against electron-hole recombination or annihilation actually depends on the energy band structure of the host conductors. While electrons and holes recombine spontaneously in semiconductors, they coexist peacefully inside semimetals without destroying one another. The stability is governed by semi-metallic energy band structures. We will deliberate on this later in Sec. VI.

For the purpose of convenience, let us introduce a new term, tetron, to refer to the 3E+1H complex, in a sense similar to the term trion. We will use tetron and 3E+1H interchangeably thereafter.

Our strategy for this paper is as follows. We will first mention relevant 2D quantum mechanical results for the Coulomb potential in Sec. III. Then as a warm-up exercise we try to find out the upper bound on the ground-state energy of a 2D trion in the heavy hole approximation, which is given in Sec. IV. Based on techniques developed for trions, we try a few fully antisymmetrized wavefunctions in Sec. V to show that states of negative energy indeed happen to tetron.

## III Quantum mechanics of an electron in 2D Coulomb potential

This problem has been worked out numerous times but we just repeat some relevant results here for quick reference.

The Hamiltonian for an electron in a 2D Coulomb potential reads

$$H=\frac{\mathrm{}^2}{2m}^2\frac{e^2}{4\pi ϵr},$$
(1)
where $`ϵ`$ is the permittivity, $`e`$ the elementary charge and $`m`$ the mass of the electron.

We try to solve the problem in a polar coordinate, where the gradient operator is given by

$$=\widehat{𝐫}\frac{}{r}+\widehat{𝜽}\frac{1}{r}\frac{}{\theta },$$
(2)
with the corresponding Laplacian

$$^2=\frac{1}{r}\frac{}{r}(r\frac{}{r})+\frac{1}{r^2}\frac{^2}{\theta ^2}.$$
(3)
The Schrödinger equation reads

$$\{\frac{\mathrm{}^2}{2m}^2\frac{e^2}{4\pi ϵr}\}\mathrm{\Psi }=E\mathrm{\Psi }.$$
(4)
We want to simplify Eq. (4) by making it dimensionless. Note that $`\mathrm{}^2/2m`$ has a dimension of (energy)(length)2 while $`e^2/4\pi ϵ`$ is in (energy)(length). We select the units (or the energy and length scales) so that both $`\mathrm{}^2/2m`$ and $`e^2/4\pi ϵ`$ are 1. Under such a convention length will be measured in the unit $`2\pi ϵ\mathrm{}^2/me^2`$ while energy is measured in $`me^4/8\pi ^2ϵ^2\mathrm{}^2`$. This unit system is summarized in Table I.

Now the Schrödinger equation takes a simple form,

$$\{^2\frac{1}{r}\}\mathrm{\Psi }=E\mathrm{\Psi }.$$
(5)
Eq. (5) has the following solutions as well known,

$$\mathrm{\Psi }_{nm}(r,\theta )=\frac{w_{nm}(r)}{\sqrt{r}}\frac{\mathrm{exp}(\pm im\theta )}{\sqrt{2\pi }},m=0,1,2,3,\mathrm{},$$
(6)
where $`w_{nm}(r)`$ is given by

$`w_{nm}(r)`$ $`=`$ $`{\displaystyle \frac{(2Kr)^{m+1/2}\mathrm{exp}(Kr)}{(2m)!}}[{\displaystyle \frac{K(n+m1/2)!}{n(nm1/2)!}}]^{1/2}`$ (8)
$`\times _1F_1(n+m+{\displaystyle \frac{1}{2}};2m+1;2Kr).`$
Here $`{}_{1}{}^{}F_{1}^{}`$ is the confluent hypergeometric function, and

$$K=\frac{1}{2n},n=\frac{1}{2},\frac{3}{2},\mathrm{},m=0,1,\mathrm{},n\frac{1}{2}.$$
(9)
The energy $`E_{nm}`$ of eigenstate $`\mathrm{\Psi }_{nm}`$ is given by

$$E_{nm}=\frac{1}{4n^2}.$$
(10)
For the purpose of convenience, we list the first four eigenstates as follows,

$`\mathrm{\Psi }_{\frac{1}{2},0}(r,\theta )`$ $`=`$ $`\sqrt{{\displaystyle \frac{2}{\pi }}}\mathrm{exp}(r),`$ (11)
$`\mathrm{\Psi }_{\frac{3}{2},0}(r,\theta )`$ $`=`$ $`\sqrt{{\displaystyle \frac{2}{27\pi }}}[1{\displaystyle \frac{2}{3}}r]\mathrm{exp}({\displaystyle \frac{r}{3}}),`$ (12)
$`\mathrm{\Psi }_{\frac{3}{2},\pm 1}(r,\theta )`$ $`=`$ $`\sqrt{{\displaystyle \frac{4}{243\pi }}}r\mathrm{exp}({\displaystyle \frac{r}{3}})\mathrm{exp}(\pm i\theta ).`$ (13)
## IV Upper bound on ground-state energy of a 2D trion

We take this problem as a warm-up exercise. This is a very illustrative example showing how to balance between the Coulomb and kinetic energy. It also gives us some intuition in choosing good trial wavefunctions and on how to deal with two-body Coulomb integrals. The basic line of reasoning follows that of H.A. Bethe, E.A. Hylleraas and others. Although there have been recent researches in the same direction, our methods in dealing with the integrals are somehow different and hopefully easier to understand.

For the purpose of simplicity we take the mass of the hole as infinity. The Schrödinger for this problem is now clearly

$$\{_1^2_2^2\frac{1}{r_1}\frac{1}{r_2}+\frac{1}{r_{12}}\}\mathrm{\Psi }=E\mathrm{\Psi }.$$
(15)
where $`_1^2`$ and $`_2^2`$ are respectively short-hands for

$$\frac{1}{r_1}\frac{}{r_1}(r_1\frac{}{r_1})+\frac{1}{r_1^2}\frac{^2}{\theta _1^2},$$
(16)
and

$$\frac{1}{r_2}\frac{}{r_2}(r_2\frac{}{r_2})+\frac{1}{r_2^2}\frac{^2}{\theta _2^2}.$$
(17)
Here subscripts refer to electrons such that $`r_1`$, $`\theta _1`$ and $`r_2`$, $`\theta _2`$ are the coordinates of e1 and e2 respectively.

We choose the following trial wavefunction

$$\mathrm{\Psi }=\mathrm{exp}[\lambda (r_1+r_2)][1+k(r_1r_2)^2],$$
(18)
where $`\lambda `$ and $`k`$ are parameters to be determined later by variational principle. The first term in this wavefunction suggests that e1 and e2 are both in $`1s`$ states. Here $`1s`$ state are not exactly the $`\mathrm{\Psi }_{\frac{1}{2},0}`$ given in Eqs. (LABEL:3psi) but rather that with a screening constant $`\lambda `$. The physical origin of $`\lambda `$ is clear: the Coulomb field of the hole felt by one electron is reduced (screened) by the presence of the other. For this reason we also require that $`\lambda `$ be positive and smaller than 1. The existence of the second term is to reduce the Coulomb interaction between e1 and e2. Note that this term is small when $`r_1`$ and $`r_2`$ are close which is another way of saying that electrons 1 and 2 tend to avoid each other. So why not just keep the second term only? The reason is that while this term is good for reducing the Coulomb energy it causes some penalty in kinetic energy when compared with the first one. Thus the viable solution is to combine the merits of the two terms, which results in the trial wavefunction (18).

To facilitate our calculation we want better management on the trial wavefunction. Note that functions like $`\mathrm{exp}[\lambda (r_1+r_2)]`$ are not normalized and would be a headache for the algebra and calculus to come. The solution is to generalize the eigenstates of Eqs. (LABEL:3psi) to accommodate the screening effect. The first four solutions now read

$`|1s`$ $`=`$ $`\sqrt{{\displaystyle \frac{2}{\pi }}}\lambda \mathrm{exp}(\lambda r),`$ (19)
$`|2s`$ $`=`$ $`\sqrt{{\displaystyle \frac{2}{27\pi }}}[\lambda {\displaystyle \frac{2}{3}}\lambda ^2r]\mathrm{exp}({\displaystyle \frac{\lambda r}{3}}),`$ (20)
$`|2p^\pm `$ $`=`$ $`\sqrt{{\displaystyle \frac{4}{243\pi }}}\lambda ^2r\mathrm{exp}({\displaystyle \frac{\lambda r}{3}})\mathrm{exp}(\pm i\theta ).`$ (21)
where we have used notations $`1s`$, $`2s`$, and $`2p^\pm `$ in place of the cumbersome $`\mathrm{\Psi }`$. The meaning of these symbols should be self-evident given a background knowledge of hydrogen eigenstates in three-dimension (3D). It is not hard to verify that these wavefunctions are now normalized.

The trial wavefunction is then instead written as

$$\mathrm{\Psi }=[1+k(r_1r_2)^2]|1s_11s_2,$$
(22)
where $`|1s_11s_2=|1s_1|1s_2`$ stands for

$$\sqrt{\frac{2}{\pi }}\lambda \mathrm{exp}(\lambda r_1)\sqrt{\frac{2}{\pi }}\lambda \mathrm{exp}(\lambda r_2),$$
(23)
and $`|1s_1`$ is the short-hand for

$$\sqrt{\frac{2}{\pi }}\lambda \mathrm{exp}(\lambda r_1).$$
(24)
Similarly we define $`|1s_2`$.

This trial wave function is symmetric in the spatial coordinates and in order to make the total wavefunction antisymmetric we only need to multiply $`\mathrm{\Psi }`$ by a spin singlet state $`𝒮_{12}^a=\frac{1}{\sqrt{2}}[|_1_2|_1_2]`$ where the subscripts are electron indices with the superscript $`a`$ indicating antisymmetry. Note that the spin wavefunction in this case will not affect the calculation of either $`\mathrm{\Psi }𝒮_{12}^a|\mathrm{\Psi }𝒮_{12}^a=\mathrm{\Psi }|\mathrm{\Psi }`$ or $`\mathrm{\Psi }𝒮_{12}^a|H|\mathrm{\Psi }𝒮_{12}^a=\mathrm{\Psi }|H|\mathrm{\Psi }`$ thus we can drop it for now. We also tabulate some often used integrals in Table II to facilitate our calculation.

The Ritz variational principle states

$$E(\mathrm{\Psi })=\frac{\mathrm{\Psi }|H|\mathrm{\Psi }}{\mathrm{\Psi }|\mathrm{\Psi }}=\mathrm{min}.,$$
(25)
where

$$\mathrm{\Psi }|H|\mathrm{\Psi }=\mathrm{\Psi }|_1^2_2^2\frac{1}{r_1}\frac{1}{r_2}+\frac{1}{r_{12}}|\mathrm{\Psi }.$$
(26)
First we need to calculate $`\mathrm{\Psi }|\mathrm{\Psi }`$,

$`\mathrm{\Psi }|\mathrm{\Psi }`$ $`=`$ $`1s_11s_2|\{1+2k(r_1^2+r_2^22r_1r_2)`$ (28)
$`+k^2[r_1^4+r_2^4+6r_1^2r_2^24r_1^3r_24r_1r_2^3]\}|1s_11s_2`$
$`=`$ $`1+2k\lambda ^2+4.5k^2\lambda ^4.`$ (29)
The first term on the right-hand-side (RHS) of Eq. (26) is

$$\mathrm{\Psi }|_1^2|\mathrm{\Psi }=_1\mathrm{\Psi }|_1\mathrm{\Psi },$$
(30)
where $`|_1\mathrm{\Psi }`$ is to be understood as $`_1|\mathrm{\Psi }`$, while $`_1\mathrm{\Psi }|`$, as a notation, is the Hermitian conjugate of $`|_1\mathrm{\Psi }`$. This identity can be proved in rectangular Cartesian coordinate using integration by parts when $`\mathrm{\Psi }`$ is zero at infinity which our trial wavefunction satisfies.

Because

$`_1|\mathrm{\Psi }`$ $`=`$ $`\{\widehat{𝐫}_1{\displaystyle \frac{}{r_1}}+\widehat{𝜽}_1{\displaystyle \frac{1}{r_1}}{\displaystyle \frac{}{\theta _1}}\}|1s_11s_2`$ (31)
$`=`$ $`\{\lambda [1+k(r_1^2+r_2^22r_1r_2)]+2k(r_1r_2)\}\widehat{𝐫}_1|1s_11s_2`$ (32)
$`=`$ $`\lambda \mathrm{\Psi }\widehat{𝐫}_1+2k(r_1r_2)\widehat{𝐫}_1|1s_11s_2,`$ (33)
we have

$`_1\mathrm{\Psi }|_1\mathrm{\Psi }`$ $`=`$ $`\lambda ^2\mathrm{\Psi }|\mathrm{\Psi }+4k^2(3\lambda ^22\lambda ^2)`$ (34)
$`=`$ $`\lambda ^2+2k+8.5k^2\lambda ^2.`$ (35)
Since $`\mathrm{\Psi }`$ is symmetric in $`r_1`$ and $`r_2`$ it should not be a surprise that the second term on RHS of Eq. 26, $`\mathrm{\Psi }|_2^2|\mathrm{\Psi }`$ is equal to $`\mathrm{\Psi }|_1^2|\mathrm{\Psi }`$.

The third term in $`\mathrm{\Psi }|H|\mathrm{\Psi }`$ \[i.e., Eq. (26)\] is

$`\mathrm{\Psi }|{\displaystyle \frac{1}{r_1}}|\mathrm{\Psi }`$ $`=`$ $`1s_11s_2|{\displaystyle \frac{1}{r_1}}\{1+2k(r_1^2+r_2^22r_1r_2)`$ (37)
$`+k^2[r_1^4+r_2^4+6r_1^2r_2^24r_1^3r_24r_1r_2^3]\}|1s_11s_2`$
$`=`$ $`(2\lambda +4k\lambda ^1+9k^2\lambda ^3)`$ (38)
to which the fourth term $`\mathrm{\Psi }|\frac{1}{r_2}|\mathrm{\Psi }`$ is equal because of the symmetry in $`r_1`$ and $`r_2`$.

The last term of $`\mathrm{\Psi }|H|\mathrm{\Psi }`$ is a two-body Coulomb integral involving elliptic integrals, and

$`\mathrm{\Psi }|{\displaystyle \frac{1}{r_{12}}}|\mathrm{\Psi }`$ $`=`$ $`1s_11s_2|{\displaystyle \frac{1}{r_{12}}}\{1+2k(r_1^2+r_2^22r_1r_2)`$ (40)
$`+k^2[r_1^4+r_2^4+6r_1^2r_2^24r_1^3r_24r_1r_2^3]\}|1s_11s_2`$
$`=`$ $`1s_11s_2|{\displaystyle \frac{1}{r_{12}}}\{1+2k(2r_1^22r_1r_2)`$ (42)
$`+k^2[2r_1^4+6r_1^2r_2^28r_1^3r_2]\}|1s_11s_2,`$
where we have used identities

$`1s_11s_2|{\displaystyle \frac{r_2^2}{r_{12}}}|1s_11s_2`$ $`=`$ $`1s_11s_2|{\displaystyle \frac{r_1^2}{r_{12}}}|1s_11s_2,`$ (43)
$`1s_11s_2|{\displaystyle \frac{r_2^4}{r_{12}}}|1s_11s_2`$ $`=`$ $`1s_11s_2|{\displaystyle \frac{r_1^4}{r_{12}}}|1s_11s_2,`$ (44)
$`1s_11s_2|{\displaystyle \frac{r_1r_2^3}{r_{12}}}|1s_11s_2`$ $`=`$ $`1s_11s_2|{\displaystyle \frac{r_1^3r_2}{r_{12}}}|1s_11s_2,`$ (45)
which are easy to understand because of the symmetry of these integrals in $`r_1`$ and $`r_2`$.

As an example, we explicitly show in the following how to work out the first term on RHS of Eq. (42).

$`1s_11s_2|{\displaystyle \frac{1}{r_{12}}}|1s_11s_2`$ (47)
$`=({\displaystyle \frac{2}{\pi }})^2\lambda ^4{\displaystyle _0^{\mathrm{}}}r_1𝑑r_1{\displaystyle _0^{\mathrm{}}}r_2𝑑r_2\mathrm{exp}[2\lambda (r_1+r_2)]{\displaystyle _0^{2\pi }}𝑑\theta _1{\displaystyle _0^{2\pi }}𝑑\theta _2{\displaystyle \frac{1}{\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}\theta _{12}}}},`$ (48)
where $`\theta _{12}=\theta _1\theta _2`$.

We want to change the integration variables in the integral

$$_0^{2\pi }𝑑\theta _1_0^{2\pi }𝑑\theta _2\frac{1}{\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}\theta _{12}}}$$
(49)
from $`\theta _1`$ and $`\theta _2`$ to $`\theta _{12}`$ and $`\theta _2`$. To understand how this is done, we refer to Fig. 3. Shown in Fig. 3 is a plane formed by $`\theta _1`$ and $`\theta _2`$. The integration limits of $`\theta _1`$ and $`\theta _2`$ in Eq. 49 or the integration domain corresponds to the square ABCD.

Because the integrand in integral (49)

$$\frac{1}{\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}\theta _{12}}}$$
(50)
is periodic in $`\theta _1`$ with period $`2\pi `$, the integration over triangle ABC is the same as that over triangle DCE. Note that DCE is the translate of ABC by $`2\pi `$ in the $`\theta _1`$ direction. Thus we can now change the domain of integration from the square ABCD to the parallelogram ACED. Since the segments AC and DE respectively satisfy the relation $`\theta _1\theta _2=\theta _{12}=0`$ and $`\theta _1\theta _2=\theta _{12}=2\pi `$, if we transform integral variables to $`\theta _{12}`$ and $`\theta _2`$ along with the change in integration domain, we will have

$`{\displaystyle _0^{2\pi }}𝑑\theta _1{\displaystyle _0^{2\pi }}𝑑\theta _2{\displaystyle \frac{1}{\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}\theta _{12}}}}`$ (51)
$`=`$ $`{\displaystyle _0^{2\pi }}𝑑\theta _2{\displaystyle _0^{2\pi }}𝑑\theta _{12}{\displaystyle \frac{1}{\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}\theta _{12}}}}|{\displaystyle \frac{(\theta _1,\theta _2)}{(\theta _2,\theta _{12})}}|,`$ (52)
where $`|\frac{(\theta _1,\theta _2)}{(\theta _2,\theta _{12})}|`$ is the absolute of Jacobian determinant

$$\frac{(\theta _1,\theta _2)}{(\theta _2,\theta _{12})}=\left|\begin{array}{cc}\frac{\theta _1}{\theta _2}& \frac{\theta _1}{\theta _{12}}\\ \frac{\theta _2}{\theta _2}& \frac{\theta _2}{\theta _{12}}\end{array}\right|=\left|\begin{array}{cc}1& 1\\ 1& 0\end{array}\right|=1,$$
(53)
which is identically $`1`$ in its absolute value.

Thus integral (49) now satisfies

$`{\displaystyle _0^{2\pi }}𝑑\theta _2{\displaystyle _0^{2\pi }}𝑑\theta _{12}{\displaystyle \frac{1}{\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}\theta _{12}}}}`$ (54)
$`=`$ $`2\pi {\displaystyle _0^{2\pi }}𝑑\theta _{12}{\displaystyle \frac{1}{\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}\theta _{12}}}}.`$ (55)
Inserting Eq. (55) into Eq. (48), we get

$`1s_11s_2|{\displaystyle \frac{1}{r_{12}}}|1s_11s_2`$ (56)
$`={\displaystyle \frac{8}{\pi }}\lambda ^4{\displaystyle _0^{\mathrm{}}}𝑑r_1{\displaystyle _0^{\mathrm{}}}𝑑r_2\mathrm{exp}[2\lambda (r_1+r_2)]r_1r_2{\displaystyle _0^{2\pi }}𝑑\theta _{12}{\displaystyle \frac{1}{\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}\theta _{12}}}}.`$ (57)
Eq. (57) can be further simplified by the following change of variables,

$`r_1=r,`$ (58)
$`r_2=ur,`$ (59)
with corresponding Jacobian

$$\frac{(r_1,r_2)}{(r,u)}=\left|\begin{array}{cc}\frac{r_1}{r}& \frac{r_1}{u}\\ \frac{r_2}{r}& \frac{r_2}{u}\end{array}\right|=\left|\begin{array}{cc}1& 0\\ u& r\end{array}\right|=r.$$
(60)
This transformation changes RHS of Eq. (57) to

$$\frac{8}{\pi }\lambda ^4_0^{\mathrm{}}u𝑑u_0^{2\pi }𝑑\theta _{12}\frac{1}{\sqrt{1+u^22u\mathrm{cos}\theta _{12}}}_0^{\mathrm{}}𝑑rr^2\mathrm{exp}[2\lambda (1+u)r],$$
(61)
which can be first integrated with respect to $`r`$ to give

$$\frac{8}{\pi }\lambda ^42!_0^{\mathrm{}}[\frac{2}{\lambda (1+u)}]^3u𝑑u_0^{2\pi }𝑑\theta _{12}\frac{1}{\sqrt{1+u^22u\mathrm{cos}\theta _{12}}},$$
(62)
where

$$_0^{2\pi }𝑑\theta _{12}\frac{u}{\sqrt{1+u^22u\mathrm{cos}\theta _{12}}}=\frac{4𝑲(\frac{4u}{12u+u^2})}{\sqrt{12u+u^2}}$$
(63)
after integration over $`\theta _{12}`$, with $`𝑲`$ the complete elliptic integral of the first kind. Integral (62) is therefore reduced to

$$\frac{16}{\pi }\lambda _0^{\mathrm{}}[\frac{2}{1+u}]^3\frac{4u𝑲(\frac{4u}{12u+u^2})}{\sqrt{12u+u^2}}𝑑u$$
(64)
which can be integrated numerically to give the final answer

$`1s_11s_2|{\displaystyle \frac{1}{r_{12}}}|1s_11s_2=1.1780972\lambda .`$ (65)
Similarly we work out other terms in Eq. (42)

$`2k1s_11s_2|{\displaystyle \frac{2r_1^22r_1r_2}{r_{12}}}|1s_11s_2=1.0308351k\lambda ^1,`$ (66)
and

$`k^21s_11s_2|{\displaystyle \frac{2r_1^4+6r_1^2r_2^28r_1^3r_2}{r_{12}}}|1s_11s_2`$ (67)
$`=1.5186410k^2\lambda ^3.`$ (68)
Substituting Eqs. (65),(66) and (68) into (42) and combining terms from (35), (38), we find the expectation value of the Hamiltonian

$`\mathrm{\Psi }|H|\mathrm{\Psi }`$ (69)
$`=2\lambda ^2+4k+17k^2\lambda ^2`$ (70)
$`2.8219028\lambda 6.9691649k\lambda ^116.4813590k^2\lambda ^3,`$ (71)
where the first three terms are from kinetic energy and the last three are Coulomb ones.

We want to minimize the energy

$`E`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Psi }|H|\mathrm{\Psi }}{\mathrm{\Psi }|\mathrm{\Psi }}}`$ (72)
$`=`$ $`{\displaystyle \frac{2\lambda ^2+4k+17k^2\lambda ^22.8219028\lambda 6.9691649k\lambda ^116.4813590k^2\lambda ^3}{1+2k\lambda ^2+4.5k^2\lambda ^4}},`$ (73)
with respect to $`\lambda `$ and $`k`$, which has numerical results as follows,

$$E_{min}=1.08147;k=0.111270;\lambda =0.693001.$$
(74)
This value, as an upper bound of the ground-state energy for our specific trion, is lower than that of an exciton with heavy hole, which is $`1`$. The trion is therefore more stable than exciton and it forms spontaneously given one exciton and an extra electron.

Our calculation shows that the trion binding energy in two-dimension $`>0.08147`$ 2D Rydberg, larger than that of 3D one which is $`0.055`$ 3D Rydberg. This is not a surprise because it is well known that lower dimensionality favors the formation of bound states. As a result, pure excitons in 2D have a binding energy 4 times as big as that in 3D.

We are so far confined to heavy hole limit. It is interesting to see quite some reports showing trion states more stable than those of exciton regardless of the electron-hole mass ratio although the author of this paper have not obtained such results using a slightly different but more rigorous Hamiltonian. Not reaching a low energy after a few trials in variational methods does not rule out its actual existence and this problem will not be further probed here because it is off our main topic for this paper. Such issue aside, we have been familiar with the techniques needed for further development and let us now investigate the possible negative-energy state of tetron or 3E+1H.

## V The variational treatment for tetron

### A The Hamiltonian

We assume that all particles have the same mass $`m`$. The full Hamiltonian is

$`H`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}(\underset{¯}{}_0^2+\underset{¯}{}_1^2+\underset{¯}{}_2^2+\underset{¯}{}_3^2)`$ (75)
$`+`$ $`{\displaystyle \frac{e^2}{4\pi ϵ}}({\displaystyle \frac{1}{\underset{¯}{r}_{12}}}+{\displaystyle \frac{1}{\underset{¯}{r}_{23}}}+{\displaystyle \frac{1}{\underset{¯}{r}_{31}}}{\displaystyle \frac{1}{\underset{¯}{r}_{10}}}{\displaystyle \frac{1}{\underset{¯}{r}_{20}}}{\displaystyle \frac{1}{\underset{¯}{r}_{30}}}),`$ (76)
where subscripts 1, 2, and 3 are indices to the three electrons while 0 designates the hole. $`\underset{¯}{r}_{12}`$ is the distance between particles 1 and 2. The meanings of $`\underset{¯}{r}_{23}`$, $`\underset{¯}{r}_{31}`$, $`\underset{¯}{r}_{10}`$, $`\underset{¯}{r}_{20}`$, and $`\underset{¯}{r}_{30}`$ are similarly defined. Using the underlined $`r`$ instead of pure $`r`$ is simply for notational purpose; we will in a moment transform the coordinates into a more convenient form. Under the current notation $`\underset{¯}{}_0^2`$ is given by

$$\underset{¯}{}_0^2=\frac{1}{\underset{¯}{r}_0}\frac{}{\underset{¯}{r}_0}(\underset{¯}{r}_0\frac{}{\underset{¯}{r}_0})+\frac{1}{\underset{¯}{r}_0^2}\frac{^2}{\underset{¯}{\theta }_0^2}$$
(77)
and $`\underset{¯}{}_1^2`$, $`\underset{¯}{}_2^2`$, and $`\underset{¯}{}_3^2`$ are understood similarly.

Performing the dimensionless unit transform same as made before in Sec. III , the Hamiltonian is now

$`H`$ $`=`$ $`(\underset{¯}{}_0^2+\underset{¯}{}_1^2+\underset{¯}{}_2^2+\underset{¯}{}_3^2)`$ (79)
$`+({\displaystyle \frac{1}{\underset{¯}{r}_{12}}}+{\displaystyle \frac{1}{\underset{¯}{r}_{23}}}+{\displaystyle \frac{1}{\underset{¯}{r}_{31}}}{\displaystyle \frac{1}{\underset{¯}{r}_{10}}}{\displaystyle \frac{1}{\underset{¯}{r}_{20}}}{\displaystyle \frac{1}{\underset{¯}{r}_{30}}}).`$
We introduce a new coordinate system

$`𝐫_0`$ $`=`$ $`{\displaystyle \frac{\underset{¯}{𝐫}_0+\underset{¯}{𝐫}_1+\underset{¯}{𝐫}_2+\underset{¯}{𝐫}_3}{4}},`$ (80)
$`𝐫_1`$ $`=`$ $`{\displaystyle \frac{\underset{¯}{𝐫}_0+\underset{¯}{𝐫}_1\underset{¯}{𝐫}_2\underset{¯}{𝐫}_3}{2}},`$ (81)
$`𝐫_2`$ $`=`$ $`{\displaystyle \frac{\underset{¯}{𝐫}_0\underset{¯}{𝐫}_1+\underset{¯}{𝐫}_2\underset{¯}{𝐫}_3}{2}},`$ (82)
$`𝐫_3`$ $`=`$ $`{\displaystyle \frac{\underset{¯}{𝐫}_0\underset{¯}{𝐫}_1\underset{¯}{𝐫}_2+\underset{¯}{𝐫}_3}{2}},`$ (83)
which are treated as vector transforms.

These transforms result into the following relations,

$`𝐫_0𝐫_1`$ $`=`$ $`\underset{¯}{𝐫}_2+\underset{¯}{𝐫}_3,`$ (84)
$`𝐫_0𝐫_2`$ $`=`$ $`\underset{¯}{𝐫}_3+\underset{¯}{𝐫}_1,`$ (85)
$`𝐫_0𝐫_3`$ $`=`$ $`\underset{¯}{𝐫}_1+\underset{¯}{𝐫}_2,`$ (86)
and

$`𝐫_1𝐫_2`$ $`=`$ $`\underset{¯}{𝐫}_1\underset{¯}{𝐫}_2,`$ (87)
$`𝐫_2𝐫_3`$ $`=`$ $`\underset{¯}{𝐫}_2\underset{¯}{𝐫}_3,`$ (88)
$`𝐫_3𝐫_1`$ $`=`$ $`\underset{¯}{𝐫}_3\underset{¯}{𝐫}_1.`$ (89)
With the transformed coordinate system in place, the Hamiltonian now takes a new form,

$`H`$ $`=`$ $`({\displaystyle \frac{1}{4}}_0^2+_1^2+_2^2+_3^2)`$ (90)
$`+`$ $`({\displaystyle \frac{1}{r_{12}}}+{\displaystyle \frac{1}{r_{23}}}+{\displaystyle \frac{1}{r_{31}}}{\displaystyle \frac{1}{r_{1+2}}}{\displaystyle \frac{1}{r_{2+3}}}{\displaystyle \frac{1}{r_{3+1}}}).`$ (91)
where for example

$$_0^2=\frac{1}{r_0}\frac{}{r_0}(r_0\frac{}{r_0})+\frac{1}{r_0^2}\frac{^2}{\theta _0^2},$$
(92)
and

$$r_{12}=\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}\theta _{12}},$$
(93)
while $`r_{1+2}`$ is the short-hand for $`|𝐫_1+𝐫_2|`$ or $`\sqrt{r_1^2+r_2^2+2r_1r_2\mathrm{cos}\theta _{12}}`$.

$`𝐫_0`$ now refers to the center-of-mass (c.m.). After discarding the irrelevant first term on RHS of Eq. (91), we reduce the internal Hamiltonian in the c.m. frame to

$`H`$ $`=`$ $`(_1^2+_2^2+_3^2)`$ (95)
$`+({\displaystyle \frac{1}{r_{12}}}+{\displaystyle \frac{1}{r_{23}}}+{\displaystyle \frac{1}{r_{31}}}{\displaystyle \frac{1}{r_{1+2}}}{\displaystyle \frac{1}{r_{2+3}}}{\displaystyle \frac{1}{r_{3+1}}}),`$
which is effectively a three-body problem.

### B Spin-singlet trial wavefunction of $`1s`$, $`2p^+`$, and $`2p^{}`$ states

It is easy to see that if the trial wavefunction for the three-body problem is antisymmetrized, the corresponding wavefunction for the our original four-body problem will be antisymmetric as well under the permutation of the original electron coordinates. Now we want to see how to make a three-body antisymmetrized wavefunctions.

The trial wavefunctions will be initially made only from the first four eigenstates for the purpose of simplicity.

$`|1s`$ $`=`$ $`\sqrt{{\displaystyle \frac{2}{\pi }}}\lambda \mathrm{exp}(\lambda r),`$ (96)
$`|2s`$ $`=`$ $`\sqrt{{\displaystyle \frac{2}{27\pi }}}[\lambda {\displaystyle \frac{2}{3}}\lambda ^2r]\mathrm{exp}({\displaystyle \frac{\lambda r}{3}}),`$ (97)
$`|2p^\pm `$ $`=`$ $`\sqrt{{\displaystyle \frac{4}{243\pi }}}\lambda ^2r\mathrm{exp}({\displaystyle \frac{\lambda r}{3}})\mathrm{exp}(\pm i\theta ).`$ (98)
which are copied from the previous section for quick reference here. As a notation convenience which will be extensively used later, we also define

$$|2p=\sqrt{\frac{4}{243\pi }}\lambda ^2r\mathrm{exp}(\frac{\lambda r}{3}).$$
(100)
Note that $`2p|2p=1`$ although $`|2p`$ is not necessarily orthogonal to other eigenstates.

We now try to make an antisymmetrized wavefunction $`\mathrm{\Psi }`$ out of Eqs. (LABEL:62). $`\mathrm{\Psi }`$ should satisfy our main object, $`\mathrm{\Psi }|H|\mathrm{\Psi }/\mathrm{\Psi }|\mathrm{\Psi }<0`$.

The nice symmetry property of the Hamiltonian (95) helps us rule out $`\mathrm{\Psi }`$ only made of s states. To see this, we need to understand that the kinetic terms in Hamiltonian are positive definite, or $`\mathrm{\Psi }|(_1^2+_2^2+_3^2)|\mathrm{\Psi }/\mathrm{\Psi }|\mathrm{\Psi }>0`$. Thus the only possibility left for a negative energy state is by a negative Coulomb term $`\mathrm{\Psi }|(\frac{1}{r_{12}}+\frac{1}{r_{23}}+\frac{1}{r_{31}}\frac{1}{r_{1+2}}\frac{1}{r_{2+3}}\frac{1}{r_{3+1}})|\mathrm{\Psi }/\mathrm{\Psi }|\mathrm{\Psi }`$. We only need to consider terms like $`\mathrm{\Psi }|(\frac{1}{r_{12}}\frac{1}{r_{1+2}})|\mathrm{\Psi }/\mathrm{\Psi }|\mathrm{\Psi }`$ because the total Coulomb energy is just three times as much because of the symmetry of $`\mathrm{\Psi }`$ in $`𝐫_1,𝐫_2,𝐫_3`$. Now we show and exploit a further symmetry hidden in $`r_{12}`$ and $`r_{1+2}`$.

Note that

$`r_{12}`$ $`=`$ $`\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}\theta _{12}},`$ (101)
$`r_{1+2}`$ $`=`$ $`\sqrt{r_1^2+r_2^2+2r_1r_2\mathrm{cos}\theta _{12}}`$ (102)
$`=`$ $`\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}(\pi +\theta _{12})}.`$ (103)
Both $`r_{12}`$ and $`r_{1+2}`$ are periodic in $`\theta _1`$ with a period $`2\pi `$. By the same argument we made in Sec. IV for two-body Coulomb integrals, we can integrate with respect to $`\theta _{12}`$ and $`\theta _2`$, both from $`0`$ to $`2\pi `$. However $`r_{1+2}`$ is identical to $`r_{12}`$ except for a phase shift of $`\pi `$ in $`\theta _{12}`$. Since $`\mathrm{cos}\theta _{12}`$ has a period of $`2\pi `$ in $`\theta _{12}`$, $`\mathrm{\Psi }|\frac{1}{r_{12}}|\mathrm{\Psi }/\mathrm{\Psi }|\mathrm{\Psi }`$ would be equal to $`\mathrm{\Psi }|\frac{1}{r_{1+2}}|\mathrm{\Psi }/\mathrm{\Psi }|\mathrm{\Psi }`$ if $`\mathrm{\Psi }`$ does not have any angular dependency such as those given by pure $`s`$ states. The total Coulomb energy would be 0 and we have no way whatsoever to obtain a negative energy. We will give a more general theorem later and exploit similar properties to greatly simplify our calculation.

So one of the states chosen in $`\mathrm{\Psi }`$ has to be $`|2p^+`$ or $`|2p^{}`$. From our trial-and-error experience, we find that it is generally better to have two $`2p`$ states and one $`s`$ state.

As the first trial we want to test wavefunctions made out of $`|1s`$, $`|2p^+`$ and $`|2p^{}`$.

By using variational approach, we also want to investigate the spin structure of such 3E+1H complex. Since there is no spin-orbit coupling term in the Hamiltonian, we need only to specify the total spin of three electrons for any trial wavefunction. Basically there are two possible options for three electrons: either they are all lined up in the same direction or one of them is against the other two. The hole is also a fermion, whose spin tends to be opposite to at least one of the electrons. Thus the total spin of this complex is either 0 or $`\mathrm{}`$, corresponding to spin singlet or triplet states. We would consider this explicitly in our calculation.

First let us start with three orthogonal states $`|1s`$, $`|2p^+`$ and $`|2p^{}`$. The arrows designate spin states. The complex is now in a spin singlet state because the mixed spin wavefunction $`|=\frac{1}{\sqrt{2}}|+\frac{1}{\sqrt{2}}|`$ should combine with that of the hole to give a singlet for one electron-hole pair. Writing $`|`$ also means that we do not care about the exact spin state; it can be either up or down.

The corresponding Slater determinant is

$`\sqrt{{\displaystyle \frac{1}{6}}}\left|\begin{array}{ccc}|1s_1_1& |2p_1^+_1& |2p_1^{}_1\\ |1s_2_2& |2p_2^+_2& |2p_2^{}_2\\ |1s_3_3& |2p_3^+_3& |2p_3^{}_3\end{array}\right|`$ (107)
$`=\sqrt{{\displaystyle \frac{1}{6}}}\left|\begin{array}{cc}|2p_2^+_2& |2p_2^{}_2\\ |2p_3^+_3& |2p_3^{}_3\end{array}\right||1s_1_1+\sqrt{{\displaystyle \frac{1}{6}}}\left|\begin{array}{cc}|2p_1^+_1& |2p_1^{}_1\\ |2p_2^+_2& |2p_2^{}_2\end{array}\right||1s_3_3\sqrt{{\displaystyle \frac{1}{6}}}\left|\begin{array}{cc}|2p_1^+_1& |2p_1^{}_1\\ |2p_3^+_3& |2p_3^{}_3\end{array}\right||1s_2_2`$ (114)
$`=\sqrt{{\displaystyle \frac{1}{6}}}\mathbf{\{}|2p_22p_3\{\mathrm{exp}[i(\theta _2\theta _3)]|_2_3\mathrm{exp}[i(\theta _3\theta _2)]|_3_2\}|1s_1_1`$ (115)
$`+|2p_12p_2\{\mathrm{exp}[i(\theta _1\theta _2)]|_1_2\mathrm{exp}[i(\theta _2\theta _1)]|_2_1\}|1s_3_3`$ (116)
$`|2p_12p_3\{\mathrm{exp}[i(\theta _1\theta _3)]|_1_3\mathrm{exp}[i(\theta _3\theta _1)]|_3_1\}|1s_2_2\mathbf{\}},`$ (117)
where for example $`_1`$ indicates a mixed spin state for electron 1 while $`2p_2^+`$ says that electron 2 is in a $`2p^+`$ state and so on so forth. Subscripts 1, 2, and 3 as usual are electron indices.

There will be 6 terms in this antisymmetrized trial wavefunction, which would render follow-on manipulation difficult. Here we do a small trick by adding to the above the following Slater determinant, which is similar to Eq. (117) except for the switch between $`2p^+`$ and $`2p^{}`$ states.

$`\sqrt{{\displaystyle \frac{1}{6}}}\left|\begin{array}{ccc}|1s_1_1& |2p_1^{}_1& |2p_1^+_1\\ |1s_2_2& |2p_2^{}_2& |2p_2^+_2\\ |1s_3_3& |2p_3^{}_3& |2p_3^+_3\end{array}\right|`$ (121)
$`=\sqrt{{\displaystyle \frac{1}{6}}}\left|\begin{array}{cc}|2p_2^{}_2& |2p_2^+_2\\ |2p_3^{}_3& |2p_3^+_3\end{array}\right||1s_1_1+\sqrt{{\displaystyle \frac{1}{6}}}\left|\begin{array}{cc}|2p_1^{}_1& |2p_1^+_1\\ |2p_2^{}_2& |2p_2^+_2\end{array}\right||1s_3_3\sqrt{{\displaystyle \frac{1}{6}}}\left|\begin{array}{cc}|2p_1^{}_1& |2p_1^+_1\\ |2p_3^{}_3& |2p_3^+_3\end{array}\right||1s_2_2`$ (128)
$`=\sqrt{{\displaystyle \frac{1}{6}}}\mathbf{\{}|2p_22p_3\{\mathrm{exp}[i(\theta _3\theta _2)]|_2_3\mathrm{exp}[i(\theta _2\theta _3)]|_3_2\}|1s_1_1`$ (129)
$`+|2p_12p_2\{\mathrm{exp}[i(\theta _2\theta _1)]|_1_2\mathrm{exp}[i(\theta _1\theta _2)]|_2_1\}|1s_3_3`$ (130)
$`|2p_12p_3\{\mathrm{exp}[i(\theta _3\theta _1)]|_1_3\mathrm{exp}[i(\theta _1\theta _3)]|_3_1\}|1s_2_2\mathbf{\}}.`$ (131)
The result is the following after normalization,

$`\mathrm{\Psi }={\displaystyle \frac{2}{\sqrt{6}}}\{`$ $`|2p_12p_21s_3\mathrm{cos}\theta _{12}|𝒮_{12,3}^a`$ (134)
$`+|2p_22p_31s_1\mathrm{cos}\theta _{23}|𝒮_{23,1}^a`$
$`+|2p_32p_11s_2\mathrm{cos}\theta _{31}|𝒮_{31,2}^a\},`$
where the spin wavefunction $`|𝒮_{12,3}^a=\frac{\sqrt{2}}{2}(|_1_2|_2_1)|_3`$ which is antisymmetric in the first two indices 1 and 2, and note that

$$𝒮_{12,3}^a|𝒮_{23,1}^a=𝒮_{12,3}^a|𝒮_{31,2}^a=𝒮_{23,1}^a|𝒮_{31,2}^a=\frac{1}{2}.$$
(135)
$`\mathrm{\Psi }`$ is already normalized and thus $`\mathrm{\Psi }|\mathrm{\Psi }=1`$. The kinetic energy of this trial wavefunction $`\mathrm{\Psi }|_1^2_2^2_3^2|\mathrm{\Psi }`$ is easy to work out

$$\mathrm{\Psi }|_1^2_2^2_3^2|\mathrm{\Psi }=\frac{11}{9}\lambda ^2.$$
(136)
Actually this can be figured out without explicit calculation. We have three particles in $`1s`$, $`2p^+`$ and $`2p^{}`$ states whose kinetic energy are respectively $`\lambda ^2`$, $`\lambda ^2/9`$ and $`\lambda ^2/9`$, thus the total kinetic energy should be the sum, which gives $`11\lambda ^2/9`$.

As for the Coulomb energy, we need only to calculate this term

$$\mathrm{\Psi }|\frac{1}{r_{12}}\frac{1}{r_{1+2}}|\mathrm{\Psi },$$
(137)
because

$$\mathrm{\Psi }|\frac{1}{r_{12}}\frac{1}{r_{1+2}}|\mathrm{\Psi }=\mathrm{\Psi }|\frac{1}{r_{23}}\frac{1}{r_{2+3}}|\mathrm{\Psi }=\mathrm{\Psi }|\frac{1}{r_{31}}\frac{1}{r_{3+1}}|\mathrm{\Psi },$$
(138)
considering symmetry in electron permutation.

We then write explicitly

$`\mathrm{\Psi }|{\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_{1+2}}}|\mathrm{\Psi }`$ (139)
$`={\displaystyle \frac{4}{6}}\{2p_12p_21s_3|\mathrm{cos}^2\theta _{12}({\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_{1+2}}})|2p_12p_21s_3`$ (140)
$`+2p_22p_31s_1|\mathrm{cos}^2\theta _{23}({\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_{1+2}}})|2p_22p_31s_1`$ (141)
$`+2p_32p_11s_2|\mathrm{cos}^2\theta _{31}({\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_{1+2}}})|2p_32p_11s_2`$ (142)
$`2p_22p_31s_1|\mathrm{cos}\theta _{12}\mathrm{cos}\theta _{23}({\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_{1+2}}})|2p_12p_21s_3`$ (143)
$`2p_32p_11s_2|\mathrm{cos}\theta _{12}\mathrm{cos}\theta _{31}({\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_{1+2}}})|2p_12p_21s_3`$ (144)
$`2p_32p_11s_2|\mathrm{cos}\theta _{23}\mathrm{cos}\theta _{31}({\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_{1+2}}})|2p_32p_21s_1\}.`$ (145)
For obvious identities from symmetry

$`2p_22p_31s_1|\mathrm{cos}^2\theta _{23}({\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_{1+2}}})|2p_22p_31s_1`$ (146)
$`=2p_32p_11s_2|\mathrm{cos}^2\theta _{31}({\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_{1+2}}})|2p_32p_11s_2,`$ (147)
$`2p_22p_31s_1|\mathrm{cos}\theta _{12}\mathrm{cos}\theta _{23}({\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_{1+2}}})|2p_12p_21s_3`$ (148)
$`=2p_32p_11s_2|\mathrm{cos}\theta _{12}\mathrm{cos}\theta _{31}({\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_{1+2}}})|2p_12p_21s_3,`$ (149)
Eq. (145) is simplified to

$`\mathrm{\Psi }|{\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_{1+2}}}|\mathrm{\Psi }`$ (150)
$`={\displaystyle \frac{4}{6}}\{2p_12p_21s_3|\mathrm{cos}^2\theta _{12}({\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_{1+2}}})|2p_12p_21s_3`$ (151)
$`+22p_22p_31s_1|\mathrm{cos}^2\theta _{23}({\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_{1+2}}})|2p_22p_31s_1`$ (152)
$`22p_22p_31s_1|\mathrm{cos}\theta _{12}\mathrm{cos}\theta _{23}({\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_{1+2}}})|2p_12p_21s_3`$ (153)
$`2p_32p_11s_2|\mathrm{cos}\theta _{23}\mathrm{cos}\theta _{31}({\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_{1+2}}})|2p_22p_31s_1\}.`$ (154)
The first term on RHS of Eq. (154) will be zero. To understand this, we prove the following theorem.

Theorem: Suppose that $`f_1(r_1,r_2,r_3)`$ and $`f_2(r_1,r_2,r_3)`$ are functions only of $`r_1`$, $`r_2`$, and $`r_3`$ (i.e., do not depend on $`\theta _1`$, $`\theta _2`$, and $`\theta _3`$). Also assume that $`f_1|\frac{1}{r_{12}}|f_2`$ exists where $`f_1`$ and $`f_2`$ are short-hands for $`f_1(r_1,r_2,r_3)`$ and $`f_2(r_1,r_2,r_3)`$. Then we have

$`f_1|\mathrm{cos}^n\theta _{12}({\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_{1+2}}})|f_2`$ (155)
$`=\{\begin{array}{cc}& 0\text{, if }n\text{ is even;}\hfill \\ & f_1|2\mathrm{cos}^n\theta _{12}(\frac{1}{r_{12}})|f_2\text{, if n is odd.}\hfill \end{array}`$ (156)
Proof:

$`f_1|\mathrm{cos}^n\theta _{12}({\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_{1+2}}})|f_2`$ (157)
$`={\displaystyle _0^{\mathrm{}}}r_1𝑑r_1{\displaystyle _0^{\mathrm{}}}r_2𝑑r_2{\displaystyle _0^{\mathrm{}}}r_3𝑑r_3{\displaystyle _0^{2\pi }}𝑑\theta _1{\displaystyle _0^{2\pi }}𝑑\theta _2{\displaystyle _0^{2\pi }}𝑑\theta _3f_1f_2\mathrm{cos}^n\theta _{12}({\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_{1+2}}})`$ (158)
$`={\displaystyle _0^{\mathrm{}}}r_1𝑑r_1{\displaystyle _0^{\mathrm{}}}r_2𝑑r_2{\displaystyle _0^{\mathrm{}}}r_3𝑑r_3{\displaystyle _0^{2\pi }}𝑑\theta _3{\displaystyle _0^{2\pi }}𝑑\theta _2{\displaystyle _0^{2\pi }}𝑑\theta _{12}f_1f_2\mathrm{cos}^n\theta _{12}({\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_{1+2}}})`$ (159)
where in the second step we have changed variables from $`\theta _1`$ to $`\theta _{12}`$ for the same reason that we discussed before in Sec. IV. We consider first the integration with respect to $`\theta _{12}`$ in Eq. (159), or the following expression

$`{\displaystyle _0^{2\pi }}𝑑\theta _{12}\mathrm{cos}^n\theta _{12}({\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_{1+2}}})`$ (160)
$`={\displaystyle _0^{2\pi }}𝑑\theta _{12}\mathrm{cos}^n\theta _{12}({\displaystyle \frac{1}{\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}\theta _{12}}}}{\displaystyle \frac{1}{\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}(\pi +\theta _{12})}}}).`$ (161)
The last term on RHS of Eq. (161) is

$$_0^{2\pi }𝑑\theta _{12}\mathrm{cos}^n\theta _{12}\frac{1}{\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}(\pi +\theta _{12})}}.$$
(162)
Changing variable from $`\theta _{12}`$ to an arbitrary $`x=\pi +\theta _{12}`$ in this integral, we have

$`{\displaystyle _0^{2\pi }}𝑑\theta _{12}\mathrm{cos}^n\theta _{12}{\displaystyle \frac{1}{\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}(\pi +\theta _{12})}}}`$ (163)
$`={\displaystyle _\pi ^{3\pi }}𝑑x\mathrm{cos}^n(x\pi ){\displaystyle \frac{1}{\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}x}}}`$ (164)
$`={\displaystyle _\pi ^{3\pi }}𝑑x(1)^n\mathrm{cos}^nx{\displaystyle \frac{1}{\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}x}}}.`$ (165)
Because the integrand on RHS of Eq. (165) is a function of $`x`$ with a period $`2\pi `$, integrating from $`0`$ to $`\pi `$ is equal to integrating from $`2\pi `$ to $`3\pi `$, or

$`{\displaystyle _{2\pi }^{3\pi }}𝑑x(1)^n\mathrm{cos}^nx{\displaystyle \frac{1}{\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}x}}}`$ (166)
$`={\displaystyle _0^\pi }𝑑x(1)^n\mathrm{cos}^nx{\displaystyle \frac{1}{\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}x}}},`$ (167)
or alternatively

$`{\displaystyle _\pi ^{3\pi }}𝑑x(1)^n\mathrm{cos}^nx{\displaystyle \frac{1}{\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}x}}}`$ (168)
$`={\displaystyle _0^{2\pi }}𝑑x(1)^n\mathrm{cos}^nx{\displaystyle \frac{1}{\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}x}}}.`$ (169)
After Eq. (169) is substituted in Eq. (165),

$`{\displaystyle _0^{2\pi }}𝑑\theta _{12}\mathrm{cos}^n\theta _{12}{\displaystyle \frac{1}{\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}(\pi +\theta _{12})}}}`$ (170)
$`={\displaystyle _0^{2\pi }}𝑑x(1)^n\mathrm{cos}^nx{\displaystyle \frac{1}{\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}x}}}`$ (171)
$`={\displaystyle _0^{2\pi }}𝑑\theta _{12}(1)^n\mathrm{cos}^n\theta _{12}{\displaystyle \frac{1}{\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}\theta _{12}}}},`$ (172)
where in the second step variable $`x`$ is relabelled as $`\theta _{12}`$.

Inserting Eq. (172) into (161), we have

$`{\displaystyle _0^{2\pi }}𝑑\theta _{12}\mathrm{cos}^n\theta _{12}({\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_{1+2}}})`$ (173)
$`={\displaystyle _0^{2\pi }}𝑑\theta _{12}\mathrm{cos}^n\theta _{12}[{\displaystyle \frac{1}{\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}\theta _{12}}}}{\displaystyle \frac{(1)^n}{\sqrt{r_1^2+r_2^22r_1r_2\mathrm{cos}\theta _{12}}}}]`$ (174)
$`=\{\begin{array}{cc}& 0\text{, if }n\text{ is even;}\hfill \\ & _0^{2\pi }𝑑\theta _{12}\frac{2\mathrm{cos}^n\theta _{12}}{r_{12}}\text{, if }n\text{ is odd.}\hfill \end{array}`$ (175)
Then we substitute Eq. (175) back to Eq. (159) to give the theorem. This concludes our proof.

The second term on RHS of Eq. (154) or $`22p_22p_31s_1|\mathrm{cos}\theta _{23}^2(\frac{1}{r_{12}}\frac{1}{r_{1+2}})|2p_22p_31s_1`$ has the factor

$`\mathrm{cos}^2\theta _{23}`$ $`=`$ $`[\mathrm{cos}\theta _2\mathrm{cos}\theta _3+\mathrm{sin}\theta _2\mathrm{sin}\theta _3]^2`$ (176)
$`=`$ $`\mathrm{cos}^2\theta _2\mathrm{cos}^2\theta _3+\mathrm{sin}^2\theta _2\mathrm{sin}^2\theta _3`$ (178)
$`+2\mathrm{cos}\theta _2\mathrm{cos}\theta _3\mathrm{sin}\theta _2\mathrm{sin}\theta _3,`$
where last term on the right gives 0 after integration with respect to $`\theta _3`$, while the first two terms, combined with $`\frac{1}{r_{12}}\frac{1}{r_{1+2}}`$ are 0 after integration with respect to $`\theta _{12}`$ from our theorem, because there is no odd power of $`\mathrm{cos}\theta _{12}`$ in the integrand.

The third term on RHS of Eq. (154) is 0 after integration with respect to $`\theta _3`$ and we have only the last term left.

Because

$`\mathrm{cos}`$ $`\theta _{23}\mathrm{cos}\theta _{31}`$ (179)
$`=`$ $`[\mathrm{cos}\theta _2\mathrm{cos}\theta _3+\mathrm{sin}\theta _2\mathrm{sin}\theta _3][\mathrm{cos}\theta _3\mathrm{cos}\theta _1+\mathrm{sin}\theta _3\mathrm{sin}\theta _1]`$ (180)
$`=`$ $`\mathrm{cos}\theta _2\mathrm{cos}\theta _1\mathrm{cos}^2\theta _3+\mathrm{sin}\theta _2\mathrm{sin}\theta _1\mathrm{sin}^2\theta _3`$ (182)
$`+\mathrm{cos}\theta _3\mathrm{cos}\theta _1\mathrm{sin}\theta _2\mathrm{sin}\theta _3+\mathrm{cos}\theta _2\mathrm{cos}\theta _3\mathrm{sin}\theta _3\mathrm{sin}\theta _1,`$
which after integration with respect to $`\theta _3`$ gives

$`\pi (\mathrm{cos}\theta _2\mathrm{cos}\theta _1+\mathrm{sin}\theta _2\mathrm{sin}\theta _1)`$ (183)
$`=\pi \mathrm{cos}\theta _{12}.`$ (184)
Thus

$`2p_32p_11s_2|\mathrm{cos}\theta _{23}\mathrm{cos}\theta _{31}({\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_{1+2}}})|2p_32p_21s_1`$ (185)
$`=2p_11s_2|{\displaystyle \frac{\mathrm{cos}\theta _{12}}{r_{12}}}|2p_21s_1`$ (186)
which is worked out in the same fashion as integral (48) in section IV to give

$$2p_11s_2|\frac{\mathrm{cos}\theta _{12}}{r_{12}}|2p_21s_1=0.0172573\lambda .$$
(187)
The total Coulomb energy is thus $`0.0345146\lambda `$ and we have

$`E`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Psi }|H|\mathrm{\Psi }}{\mathrm{\Psi }|\mathrm{\Psi }}}`$ (188)
$`=`$ $`{\displaystyle \frac{11\lambda ^2}{9}}0.0345146\lambda `$ (189)
$`=`$ $`{\displaystyle \frac{11}{9}}(\lambda 0.0141196)^20.00019936,`$ (190)
from which the upper bound for ground-state energy is therefore $`0.00019936`$.

This confirms our conjecture that 3E+1H indeed forms a quasi-bound state. Here quasi-bound merely means that the 4 particles are relatively localized or near one to another and the total energy of this system is lower than when they are infinitely apart. Note that the 4 particle as a whole could still move freely with an effective mass of $`4m`$.

We proceed on to try other wavefunctions and want to improve our result.

### C Spin-triplet trial wavefunction of $`1s`$, $`2p^+`$, and $`2p^{}`$ states

The next logic choice is still $`1s`$, $`2p^+`$, and $`2p^{}`$ states but with a spin-triplet configuration,

$`\mathrm{\Psi }={\displaystyle \frac{2}{\sqrt{6}}}\{`$ $`|2p_12p_21s_3\mathrm{sin}\theta _{12}|𝒮_{1,23}^s`$ (193)
$`+|2p_22p_31s_1\mathrm{sin}\theta _{23}|𝒮_{23,1}^s`$
$`+|2p_32p_11s_2\mathrm{sin}\theta _{31}|𝒮_{31,2}^s\},`$
where the normalized spin wavefunction $`|𝒮_{12,3}^s`$ is symmetrized in indices 1 and 2, and describes a spin triplet state for electrons 1 and 2. $`|𝒮_{123}^s`$ may be written $`|_1_2_3=\frac{\sqrt{2}}{2}|_1_2(|_3+|_3)`$ but can also take other forms, for example $`|_1_2_3`$ or $`|\frac{1}{\sqrt{2}}(_1_2+_2_1)_3`$. Three forms of $`|𝒮_{12,3}^s`$ all satisfy $`𝒮_{12,3}^s|𝒮_{12,3}^s`$ =1 and $`𝒮_{12,3}^s|𝒮_{23,1}^s=1/2`$ and so on.

This trial wavefunction gives exactly the same result as the spin singlet trial wavefunction,

$`E`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Psi }|H|\mathrm{\Psi }}{\mathrm{\Psi }|\mathrm{\Psi }}}`$ (194)
$`=`$ $`{\displaystyle \frac{11\lambda ^2}{9}}0.0345146\lambda `$ (195)
$`=`$ $`{\displaystyle \frac{11}{9}}(\lambda 0.0141196)^20.00019936,`$ (196)
which gives $`E_{min}=0.00019936`$. We are not yet able to tell whether the spin singlet or the triplet state is more stable.

### D Spin-singlet trial wavefunction of modified $`1s`$, $`2p^+`$, and $`2p^{}`$ states

We know in such a 3E+1H complex, the strong repulsive Coulomb force is a huge penalty we have pay attention to, thus it would be nice if our trial wavefunction is small when two electron coordinates are close. Take the spin-singlet wavefunction for example,

$`\mathrm{\Psi }={\displaystyle \frac{2}{\sqrt{6}}}\{`$ $`|2p_12p_21s_3\mathrm{cos}\theta _{12}|𝒮_{12,3}^a`$ (199)
$`+|2p_22p_31s_1\mathrm{cos}\theta _{23}|𝒮_{23,1}^a`$
$`+|2p_32p_11s_2\mathrm{cos}\theta _{31}|𝒮_{31,2}^a\},`$
if we can modify $`\mathrm{\Psi }`$ so that it is small when two electrons are close, the total energy then might be lower. So we choose the following,

$`\mathrm{\Psi }={\displaystyle \frac{2}{\sqrt{6}}}\{`$ $`|2p_12p_21s_3(1\mathrm{cos}\theta _{12})|𝒮_{12,3}^a`$ (202)
$`+|2p_22p_31s_1(1\mathrm{cos}\theta _{23})|𝒮_{23,1}^a`$
$`+|2p_32p_11s_2(1\mathrm{cos}\theta _{31})|𝒮_{31,2}^a\},`$
where the first term on RHS will be zero when $`\theta _{12}=0`$, or $`\theta _1=\theta _2`$, or in other words, when electrons 1 and 2 are near.

Note that after this modification, the three particles are no longer in pure $`1s`$, $`2p^+`$, and $`2p^{}`$ states, but we may still use these words for notational purpose, and call them modified $`1s`$, $`2p^+`$, and $`2p^{}`$ states.

Then

$`\mathrm{\Psi }|\mathrm{\Psi }=`$ $`22p_12p_21s_3|12\mathrm{cos}\theta _{12}+\mathrm{cos}^2\theta _{12}|2p_12p_21s_3`$ (204)
$`22p_12p_21s_3|1\mathrm{cos}\theta _{12}\mathrm{cos}\theta _{23}+\mathrm{cos}\theta _{12}\mathrm{cos}\theta _{23}|2p_22p_31s_1`$
$`=`$ $`32\{2p_11s_3|2p_31s_1+2p_12p_21s_3|\mathrm{cos}\theta _{12}\mathrm{cos}\theta _{23}|2p_22p_31s_1\},`$ (205)
where we have used identities

$`2p_12p_21s_3|12\mathrm{cos}\theta _{12}+\mathrm{cos}^2\theta _{12}|2p_12p_21s_3`$ (206)
$`=2p_22p_31s_1|12\mathrm{cos}\theta _{23}+\mathrm{cos}^2\theta _{23}|2p_22p_31s_1`$ (207)
$`=2p_32p_11s_2|12\mathrm{cos}\theta _{31}+\mathrm{cos}^2\theta _{31}|2p_32p_11s_2,`$ (208)
and

$`2p_12p_21s_3|1\mathrm{cos}\theta _{12}\mathrm{cos}\theta _{23}+\mathrm{cos}\theta _{12}\mathrm{cos}\theta _{23}|2p_22p_31s_1`$ (209)
$`=2p_22p_31s_1|1\mathrm{cos}\theta _{23}\mathrm{cos}\theta _{31}+\mathrm{cos}\theta _{23}\mathrm{cos}\theta _{31}|2p_32p_11s_2,`$ (210)
and the like.

The exact value of $`\mathrm{\Psi }|\mathrm{\Psi }`$ is found to be $`90/32`$.

Next we want to find out the kinetic energy

$`\mathrm{\Psi }|_1^2_2^2_3^2|\mathrm{\Psi }`$ (211)
$`=3\mathrm{\Psi }|_1^2|\mathrm{\Psi }`$ (212)
$`=3_1\mathrm{\Psi }|_1\mathrm{\Psi },`$ (213)
where

$$_1=\widehat{𝐫}_1\frac{}{r_1}+\widehat{𝜽}_1\frac{1}{r_1}\frac{}{\theta _1}.$$
(214)
Because

$`_1|\mathrm{\Psi }={\displaystyle \frac{2}{\sqrt{6}}}\{`$ $`(`$ $`{\displaystyle \frac{1}{r_1}}{\displaystyle \frac{\lambda }{3}})|2p_12p_21s_3(1\mathrm{cos}\theta _{12})𝒮^a_{12,3}`$ (219)
$`\lambda |2p_22p_31s_1(1\mathrm{cos}\theta _{23})𝒮_{23,1}^a`$
$`+({\displaystyle \frac{1}{r_1}}{\displaystyle \frac{\lambda }{3}})|2p_32p_11s_2(1\mathrm{cos}\theta _{31})𝒮_{31,2}^a\}\widehat{𝐫}_1`$
$`+{\displaystyle \frac{2}{\sqrt{6}}}{\displaystyle \frac{\widehat{𝜽}_1}{r_1}}\{|2p_12p_21s_3\mathrm{sin}\theta _{12}𝒮_{12,3}^a`$
$`|2p_32p_11s_2\mathrm{sin}\theta _{31}𝒮_{31,2}^a\},`$
we have

$`_1\mathrm{\Psi }|_1\mathrm{\Psi }`$ (220)
$`={\displaystyle \frac{4}{6}}\{22p_12p_2|({\displaystyle \frac{1}{r_1^2}}{\displaystyle \frac{2\lambda }{3r_1}}+{\displaystyle \frac{\lambda ^2}{9}})(12\mathrm{cos}\theta _{12}+\mathrm{cos}^2\theta _{12})|2p_12p_2`$ (221)
$`+2p_22p_3|\lambda ^2(1+\mathrm{cos}^2\theta _{23}2\mathrm{cos}\theta _{23})|2p_22p_3`$ (222)
$`22p_11s_3|({\displaystyle \frac{\lambda ^2}{3}}{\displaystyle \frac{\lambda }{r_1}})|1s_12p_3`$ (223)
$`2p_12p_21s_3|({\displaystyle \frac{1}{r_1^2}}{\displaystyle \frac{2\lambda }{3r_1}}+{\displaystyle \frac{\lambda ^2}{9}})|2p_12p_31s_2`$ (224)
$`+22p_12p_2|{\displaystyle \frac{\mathrm{sin}^2\theta _{12}}{r_1^2}}|2p_12p_2\},`$ (225)
where we have used symmetry identities similar to Eqs. (208) and (210)to merge some terms. Those terms which go to 0 after integration of one variable, say $`\theta _2`$, are dropped as well.

After somewhat tedious but straightforward calculation, the total kinetic energy is $`\frac{1405}{432}\lambda ^2`$

The Coulomb energy of the trial wavefunction is calculated as usual though with added complexity. We just give the result here, $`0.566523\lambda `$.

Thus

$`E`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Psi }|H|\mathrm{\Psi }}{\mathrm{\Psi }|\mathrm{\Psi }}}`$ (226)
$`=`$ $`1.15638\lambda ^20.201430\lambda `$ (227)
$`=`$ $`1.15638(\lambda 0.0870951)^20.00877178,`$ (228)
which gives $`E_{min}=0.00877178`$, a much improved result over previous cases.

Also recall that we have spin-singlet state in this case. In contrast, the spin-triplet trial wavefunction of Eq. (199) cannot be so modified as Eq. (202) because the first term on RHS of Eq. (199) is already small when electrons 1 and 2 are close. This suggests that antiparallel spins be energetically more favorable than parallel ones. We will from now on focus our attention to spin-singlet state only.

### E Spin-singlet trial wavefunction of modified $`2s`$, $`2p^+`$, and $`2p^{}`$ states

The next trial wavefunction and the final one in this paper is made of modified $`2s`$, $`2p^+`$, and $`2p^{}`$ states,

$`\mathrm{\Psi }={\displaystyle \frac{2}{\sqrt{6}}}\{`$ $`|2p_12p_22s_3(1\mathrm{cos}\theta _{12})|𝒮_{12,3}^a`$ (231)
$`+|2p_22p_32s_1(1\mathrm{cos}\theta _{23})|𝒮_{23,1}^a`$
$`+|2p_32p_12s_2(1\mathrm{cos}\theta _{31})|𝒮_{31,2}^a\}.`$
The treating of such $`\mathrm{\Psi }`$ is similar to the previous case, and we just list the final results as follows,

$`E`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Psi }|H|\mathrm{\Psi }}{\mathrm{\Psi }|\mathrm{\Psi }}}`$ (232)
$`=`$ $`0.632997\lambda ^20.168735\lambda `$ (233)
$`=`$ $`0.632997(\lambda 0.133283)^20.0112447.`$ (234)
Then the upper bound for ground-state energy is now improved to $`0.0112447`$.

## VI Stability against electron-hole recombination

The stability of a trion or a tetron against electron-hole recombination or annihilation actually is not a concern if the electrons and holes are from semi-metallic band-edge. Whether electrons ($`e`$) and holes ($`h`$) should recombine or not depends on the band structure. Here we just give a very simple example: semimetals. Electrons and holes coexist peacefully without annihilating one another, in semimetals  such as Bi, As, Sb. The reason is as follows and shown in Fig. 4.

In typical semiconductors, the conduction and valence bands are separated by a positive energy gap as shown in Fig. 4(a). By jumping from the conduction to the valence band which is at a lower energy (i.e., recombining with a hole), the thermally excited electron lowers its energy. This process thus occurs spontaneously to the $`e`$-$`h`$ pair and makes it unstable. (The same process can be interpreted as a hole jumping up and note that in contrast to electron, the higher is the hole in the band diagram, the lower is its energy.) For semimetals with different band configurations, as typified in Fig. 4(b) and (c), the situation is however totally different. If the electron were to jump to the hole state and then recombine, as indicated by the arrows, it would have to take extra energy to reach there. This process is clearly unfavorable in energy and such a mechanism secures the stability of electron-hole liquid against recombination in semimetals. Thus we expect the proposed tetrons to find their best playing ground in solids with semi-metallic band structures.

## VII Stability against dissociation

As mentioned before a tetron is not stable against the dissociation of one of its electrons into infinity in vacuum. To understand this, let us perform a thought experiment. The theorem by Lieb states that the ground-state energy of a 3E+1H complex must be higher than that of a trion plus an electron which is assumed to be infinitely far away, having arbitrarily small kinetic energy. Thus after we place a quasi-bound state of 3E+1H in vacuum, the complex tends to lower its energy by moving one of the electrons away to minimize the energy. This corresponds to the dissociation (detachment) of one electron as shown in Fig. 5.

However if the tetron is not in vacuum, for example, if many other electrons and holes are present nearby, then there is hardly any gain in energy for the whole system even if an electron dissociates from one tetron; this electron would simply increase Coulomb energy with other electrons anyway even if it decreases the energy of the original tetron as shown in Fig. 6. Such cases may happen to condensed matter systems where charge carrier densities are governed by electronic structure, or some types of external controls such as dopings. The electrons cannot leave the holes for infinity because statistically there must be some electrons close to the holes. When the densities of electrons and holes are low, the situation is more like a vacuum and therefore tetrons are not stable. With an increasing charge density, however, the story could be increasingly different. Then if the densities exceed a certain threshold, the whole system would lower its energy by taking advantage of the favorable geometric configuration offered by tetrons, provided the electrons outnumber holes by a ratio near or higher than 3. Then the tetron would be more stable in this high density electron-hole liquid. Actually Nature may have already suggested a similar mechanism for us: doubly charged negative ions until recently have not been seen in gas phase, or in isolation. In contrast, they have long been known to exist in condensed matter systems such as liquids and solutions.

If the densities are to further increase by too much, the energy of the tetron state will no longer be negative after a second threshold. The reason is simple, the gain in Coulomb energy increases as $`a^1`$ while the penalty paid for the kinetic energy scales as $`a^2`$, where $`a`$ is the average distance between particles. Then we would see a system close to ordinary metal because the quasi-bound state is no longer possible.

Here we just give the above qualitative argument to show how tetron states might be materialized. Further studies are needed in a more quantitative and rigorous way.

## VIII Possible connection to high-$`T_c`$ pairing mechanism?

The tetron state acts in a whole as a boson because it starts with 4 fermions. This is significant because the net charge is $`2e`$ ($`2e`$ if we start with 3h+1e which is a more realistic scenario for most copper-oxide superconductors), and we are now left with a Cooper pair in effect. Moreover our current variational result implies a singlet state for the spin part of wavefunction, which is exactly what has been found in the Cooper pairs of high-$`T_c`$ superconductors. This entices us into the possible connection between the tetron state and high-$`T_c`$ superconductivity, although this model itself could be oversimplified without considering the crystal structure of the host material.

Charged fermion pairing models based on exchange of excitons have been proposed before by some notable theorists, often in a paradigm of perturbative approaches. What we have developed here could be a non-perturbative way of treating similar problems.

High-$`T_c`$ superconductivity found in copper oxides or cuprates are likely from the Bose-Einstein condensation of pre-formed pairs. In this regard, the tetron complex reasonably fits into such a picture, with its small size comparable to the coherence length of cuprate Cooper pairs. Also high-$`T_c`$ is intimately related to low dimensionality. Recall that the ground-state energy of a 2D exciton is four times as low as that in 3D, while $`T_c`$ decreases in cuprates when the superconductors are increasingly doped from a 2D conductor to 3D one.

There are mounting evidences showing the coexistence of electrons and holes in cuprate superconductors from key experiments such as angle-resolved photoemission spectroscopy and transport properties. Finally the qualitative discussion given in Sec. VII actually quite fits to the doping-dependent phase diagram of cuprate superconductors from underdoping to overdoping. Although the case at this moment is far from certain, the author sincerely hopes that this paper may generate some interests along this line of reasoning.

## IX Conclusion

In conclusion, various trial functions for 3E+1H complex (tetron) were worked out by rigorous variational method to give an upper limit on the ground-state energy. The calculation suggests quasi-bound states for such a complex. The possible connection of such states to a pairing mechanism of high-$`T_c`$ is discussed.