Three-Body Coulomb Problem

R. Combescot

► To cite this version:

R. Combescot. Three-Body Coulomb Problem. Physical Review X, 2017, 7 (4), pp.041035. 10.1103/PhysRevX.7.041035. hal-01651622

HAL Id: hal-01651622

https://hal.sorbonne-universite.fr/hal-01651622

Submitted on 29 Nov 2017

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Three-Body Coulomb Problem

R. Combecot

Laboratoire de Physique Statistique, Ecole Normale Supérieure, PSL Research University,
UPMC Paris 06, Université Paris Diderot, CNRS, 24 rue Lhomond, 75005 Paris, France
(Received 23 May 2017; revised manuscript received 27 July 2017; published 15 November 2017)

We present a general approach for the solution of the three-body problem for a general interaction and apply it to the case of the Coulomb interaction. This approach is exact, simple, and fast. It makes use of integral equations derived from the consideration of the scattering properties of the system. In particular, this makes full use of the solution of the two-body problem, the interaction appearing only through the corresponding known *T* matrix. In the case of the Coulomb potential, we make use of a very convenient expression for the *T* matrix obtained by Schwinger. As a check, we apply this approach to the well-known problem of the helium atom ground state and obtain a perfect numerical agreement with the known result for the ground-state energy. The wave function is directly obtained from the corresponding solution. We expect our method to be, in particular, quite useful for the trion problem in semiconductors.

DOI: 10.1103/PhysRevX.7.041035

Subject Areas: Atomic and Molecular Physics,
Condensed Matter Physics,
Nuclear Physics

I. INTRODUCTION

Few-body systems and problems [1,2] are ubiquitous in
almost all fields of physics; they arise, for example, in par-
ticle physics, nuclear physics, atomic physics, condensed-
matter physics, and so on. Since two-body problems are
easily solved analytically or numerically, the first level of
nontrivial problems arises with three-body problems.
Among the first examples in quantum mechanics has been
the helium atom, more specifically its ground-state energy,
where one deals with the Coulomb interaction. This was
first addressed by Hylleraas [3] by variational methods and
pushed recently to extraordinary precision [4]. Another
quite similar case is the H⁻ ion [5], which is remarkable for
its very weakly bound ground state and is of astrophysical
interest [6]. Yet another example is found in semiconductor
physics, where the trion, i.e., a bound state of an exciton
and an electron (or a hole) [7], is observed through its
absorption or emission spectrum [8,9]. This is again a case
where the interaction is essentially the Coulomb interac-
tion. Three-body systems arise also because they may have
their own intrinsic interest, such as the well-known Efimov
trimers [2,10], with their remarkable scaling properties,
which have been the subject of much recent activity in
nuclear physics and in cold-atom physics.

The recent surge of activity in ultracold gases [11,12],
following the achievement of Bose-Einstein condensation

in these systems, has led to a renewed interest in few-body
physics [1,2]. Indeed the situation is much simplified in
these cases because these systems are dilute and the relevant
atomic energies are very low. As a result, in most cases the
interaction can be considered essentially as a contact
interaction, and in the scattering amplitude the contributions
other than *s* wave can be safely ignored. All the possible
complexities of the interaction potential disappear and the
interaction is fully characterized by the scattering length.
There is no dependence of the scattering amplitude on wave
vectors; it depends only on energy. This makes the three-
body Schrödinger equation much simpler to solve since one
has to deal with free atoms except for a boundary condition
when two atoms are at the same position. This does not mean
that the three-body problem is so easy to solve, but this has
allowed the theoretical work to be pushed quite far, in
particular, in the case where the three bodies are in a
harmonic trap [13–16]. Similarly, with this contact inter-
action, the scattering properties are much easier to find and,
for example, the dimer-atom scattering length is obtained by
solving a one-dimensional integral equation, as initiated a
long time ago by Skorniakov and Ter-Martirosian [17] for
the neutron-deuteron problem. These problems are very
convenient to formulate in a diagrammatic formalism [18]
and to generalize to four-body problems [19], leading again
to fairly simple integral equations.

This diagrammatic approach has a further very attractive
interest. Indeed it makes full use of the solution of the
two-body problem. Actually, the interaction potential never
appears explicitly in the equations; it comes in only through
the two-body propagator corresponding to the solution of
the two-body problem. This looks to be a very reasonable
and attractive approach to the solution of the three-body

Published by the American Physical Society under the terms of
the Creative Commons Attribution 4.0 International license.
Further distribution of this work must maintain attribution to
the author(s) and the published article's title, journal citation,
and DOI.
problem: it makes much more sense to use the already known solution of the two-body problem rather than start again from the beginning, as if the two-body problem had not been solved. This feature is so attractive that it is worthwhile to explore whether it can be extended with the same advantages to the case of a general interaction potential, getting rid of the simplified contact interaction suited to cold gases. Actually, this spirit is very close to another approach to the many-body problem, the “composite boson” formalism [20], where the eigenstates of the two-body problem are taken as a new basis in which the whole many-body problem is rewritten. In this approach one again makes full use of the solution of the two-body problem. In particular, the trion problem has already been addressed within this approach [21].

It is the purpose of the present paper to explore this generalization. We find that this extension can indeed be done with minimal increase in complexity. As a result we find a new method to solve the three-body problem which is at the same time exact, simple, and fast. In practice, when we come to explicit use, we consider the specific case of the 3D Coulomb potential, which is appropriate to the case of the helium ground state we consider explicitly, and also to the case of the trion, which we have mainly in mind. This Coulomb potential case turns out to be particularly convenient since there is a simple analytic expression found by Schwinger [22] for the T matrix, which sums up the solution of the two-body Coulomb problem. However, there is no real problem to extend our method to any interaction potential, and also to any dimension D. One has merely to obtain, analytically or numerically, the corresponding T matrix for the two-body problem. An interpolation method can then be used, for example, to store the result for practical use in the numerical calculation. Beyond providing an efficient way to solve any three-body problem, we hope that this approach can be extended to the four-body problem along the same lines. But exploration of this path is naturally left for future work. Such an extension would naturally be extremely useful for many problems, in particular, to appropriately treat exciton-exciton interaction in semiconductors, which is of importance for the Bose-Einstein condensation of excitons [23].

The paper is organized as follows. In Sec. II, as an introduction, we consider the case of cold gases, where the interaction is short-ranged, and review the calculation of the atom-dimer scattering length, which contains the backbone of our procedure. Then in Sec. III, we generalize the approach to a general interaction potential, leading to an integral equation for a three-body scattering amplitude whose poles give the bound-state energies and the eigenfunctions of the three-body problem. In Sec. IV, it is shown explicitly how the wave function of a bound state is obtained from the solution of the integral equation. We then specialize to the Coulomb potential and review in Sec. V the derivation of the corresponding T matrix by Schwinger.

Finally, we make use of our results in Sec. VI to obtain the helium atom ground-state energy, which is found in perfect agreement with known results. We also give our results for the ground-state wave function. Section VII is a summary and conclusion. In many cases the technical details are given in appendixes in order to simplify and keep clear the flow of the paper.

To summarize, the present paper is devoted to presenting our approach and to checking it on a very well-known case, the helium atom ground state. Its application to other interesting cases, in particular, the energy of the trion, is left for future work. For convenience and to be definite, we prefer to adopt for our presentation the semiconductor vocabulary specific to the case of the trion, since we have it in mind, rather than keep a general, vague, and unspecific wording. Hence, our three particles are one hole and two electrons, which have in most of the paper opposite spins ↑ and ↓. The translation to any other physical situations of interest is obvious.

## II. SHORT-RANGE INTERACTION

Let us first recall what happens when we replace by a short-range interaction the Coulomb interaction between the hole, with mass $m_h$, and the electrons, with mass $m_e$, which we have mostly in mind. In addition, we omit in this section the interaction between the two electrons, and we assume the simplest situation where these electrons are identical (they have the same spin). This case is useful since this is the simplest one in our class of problems. This is basically the problem handled a long time ago by Skorniakov and Ter-Martirosian [17] to obtain the deuteron-neutron scattering length $a_3$. This is also the situation found in cold gases. Here, we treat it by making use of the diagrammatic method [18,19].

This case of the short-range interaction is quite simple because, at the low energy and wave vectors we are interested in, the dimer (or exciton) propagator does not depend on the entering or outgoing wave vectors but only on the total energy $\Omega$ and momentum $\mathbf{P}$ of the dimer. The situation is even simpler since only the energy $\Omega_r$ of the relative motion is actually entering. It is given by $\Omega_r = \Omega - \mathbf{P}^2/2M$, where $M = m_e + m_h$ is the total mass of the particles making up the dimer (or exciton). Its explicit expression is, within a factor, the scattering amplitude for the relative motion and it depends only on the scattering length $a$, which in this way is the only parameter necessary to fully characterize the interaction. Specifically, it is given diagrammatically by Fig. 1 and its explicit expression is

$$T_2(P) = \frac{2\pi}{\mu} \frac{1}{a^{-1} - \sqrt{2\mu(\mathbf{P}^2/2M - \Omega - i0_+)}}, \quad (1)$$

where $P = \{\Omega, \mathbf{P}\}$ is the energy-momentum four-vector and $\mu = m_e m_h / (m_e + m_h)$ is the reduced mass. We set $\hbar = 1$ as we do everywhere in the paper. The above
<figure><img src="image_2.png" /><figcaption>FIG. 1. Exciton or dimer propagator. The full line is an electron propagator and the dashed line a hole propagator. The wavy line corresponds to an interaction.</figcaption></figure>

expression has a single pole for $\Omega_r = -1/(2\mu a^2) \equiv -E_0$ corresponding to the single bound state of the two particles. This happens only when $a > 0$, which we assume in the following; otherwise there is no bound state and accordingly no dimer.

We now want to write an integral equation for the scattering amplitude of an electron on the exciton. Qualitatively this quantity is analogous to the $T_2$ considered just above, except that the hole is replaced by the exciton. In the case of short-range interaction the situation with respect to the variables coming in the exciton-electron vertex $T_3$ is quite simple. Let us call $P$ the total energy-momentum of the electron and the exciton, and $p$ the corresponding value for the entering electron. Hence, the energy-momentum of the exciton is $P - p$, and we know from Eq. (1) that this quantity is enough to fully characterize the entering exciton through its propagator. Similarly, if the energy-momentum of the outgoing electron is $\bar{p}$, the outgoing exciton has a energy-momentum $P - \bar{p}$. Accordingly, $T_3$ depends only on $P$, $p$, and $\bar{p}$. Note that, just as we have done above for $T_2(P)$ (see Fig. 1), we do not include in the expression for $T_3(p, \bar{p}; P)$ the entering and outgoing propagators for the electron and the exciton, since they would anyway be factored out in the equation we are looking for.

In order to obtain our equation, we note that the simplest process arising in the electron-exciton scattering is merely the exchange of the incoming electron with the one making up the exciton. More precisely, since all the interactions between the electron and the hole forming the entering exciton are already taken into account in the entering exciton propagator, the only possible interaction that can arise in $T_3$ is between the hole of the exciton and the incoming electron. However, this interaction may be followed by another one between the same particles, and another one, and so on. Summing up all these possible interactions gives an exciton propagator. The simplest case arises if this exciton propagator coincides with the outgoing exciton propagator, as described by the first term in the right-hand side of Fig. 2. We see that it corresponds indeed

to an electron exchange. However, any other process may also occur with this exciton propagator and the electron before the final state. But the sum of all these processes is precisely $T_3$ by definition. This is described by the second term in the right-hand side of Fig. 2. This leads to the integral equation represented diagrammatically in Fig. 2. It reads algebraically as

$$ T_3(p, \bar{p}; P) = -g_h(P - p - \bar{p}) - \sum_q g_h(P - p - q)g_e(q) \times T_2(P - q)T_3(q, \bar{p}; P), \quad (2) $$

where, according to Feynmann diagram rules, $\sum_q \equiv i \int dq d\omega_q / (2\pi)^4$.

Here, $g_h(p) \equiv g(\{\omega_p, \mathbf{p}\}) = 1/(\omega_p - \mathbf{p}^2/2m_h + i0_+)$ is the hole Green's function, while $g_e(q) = 1/(\omega_q - \mathbf{q}^2/2m_e + i0_+)$ is the electron Green's function. Finally, the minus signs in the right-hand side of Eq. (2) come from the fact that we are exchanging the two electrons and that this permutation of these two identical fermions implies a sign change.

We can integrate on the frequency $\omega_q$ in Eq. (2) by closing the integration contour in the lower complex half-plane for the variable $\omega_q$. Indeed from their definition, $g_h(P - p - q)$ and $T_2(P - q)$ are analytical functions of $\omega_q$ in this domain. Moreover, as given by Eq. (2) itself, $T_3(p, \bar{p}; P)$ is an analytical function of $\omega_p$ for $\mathrm{Im}\omega_p < 0$, so $T_3(q, \bar{p}; P)$ is also analytical in the lower complex half-plane for the variable $\omega_q$. Accordingly, the only singularity in this region is the simple pole coming from $g_e(q)$. Hence, by residue integration, only the on-the-shell value (i.e., evaluated for $\omega_q = \mathbf{q}^2/2m_e$) of the integrand comes in. This leads us to consider the simpler problem of finding $T_3(\{\mathbf{p}^2/2m_e, \mathbf{p}\}, \bar{p}; P)$ by restricting $p$ to also be taken on the shell.

Moreover, for the problems of physical interest, such as finding the scattering length or the ground-state energy, we do not need to consider general values for $\bar{p}$ and $P$.

<figure><img src="image_3.png" /><figcaption>FIG. 2. Diagrammatic representation of integral equation Eq. (2) for the exciton-electron scattering vertex $T_3$. Full line, electron propagator; dashed line, hole propagator; shaded line, exciton propagator.</figcaption></figure>
The scattering length corresponds to a situation where all the momenta in $T_3(p, \bar{p}; P)$ go to zero while the total energy is the exciton ground-state energy, since the energy of the scattering electron goes to zero. Similarly, we find the energy $-E$ (with $E > 0$) of the bound states of the exciton-electron system, and in particular the ground-state energy, by looking for resonances of the exciton-electron scattering amplitude when all the momenta are zero. Hence, we can restrict ourselves to the specific case $\bar{p} = 0$ and $\mathbf{P} = \mathbf{0}$, so that $P = \{-E, \mathbf{0}\}$. On the other hand, we cannot set from the start $\mathbf{p} = \mathbf{0}$, since in this case we could not write an integral equation. We have to consider $\mathbf{p} \neq \mathbf{0}$, and after having found the solution possibly let $\mathbf{p}$ go to zero.

As a result, setting $T_3(\{\mathbf{p}^2/2m_e, \mathbf{p}\}, 0; \{-E, \mathbf{0}\}) = t(\mathbf{p})$ (we do not write explicitly the dependence on $E$) we obtain from Eq. (2) the simpler integral equation:

$$t(\mathbf{p}) = \frac{2\mu}{2\mu E + \mathbf{p}^2} + \frac{2\mu}{(2\pi)^3} \int d\mathbf{q} \frac{T_2(\{-E + \mathbf{q}^2/2m_e), \mathbf{q}\})}{2\mu E + \mathbf{p}^2 + \mathbf{q}^2 + 2\mu \mathbf{p} \cdot \mathbf{q} / m_h} t(\mathbf{q}). \quad (3)$$

Actually, since $t(\mathbf{p})$ depends only on the single variable $|\mathbf{p}|$, as it is obvious by rotational invariance, the angular integrations are easily performed explicitly in the right-hand side, and we are left with a single variable integration. The integral equation is very easily solved numerically. For example, for the scattering length considered by Skorniakov and Ter-Martirosian [17], all the fermion masses are equal $m_e = m_h = m$, and one has to set $E = E_0 = 1/ma^2$. The scattering length $a_3$ is then related to the solution by $a_3 = 8t(0)/(3ma)$. One finds $a_3 = 1.18a$.

### III. GENERAL INTERACTION POTENTIAL

We consider now a general interaction potential $V(\mathbf{r})$, with $\mathbf{r} = \mathbf{r}_e - \mathbf{r}_h$, and only later on specialize to the Coulomb interaction. Such a general case immediately brings formal complications, which remain in the Coulomb case. They appear as soon as we consider the exciton propagator $T_2$. Indeed, the entering and outgoing wave vectors are relevant variables, whereas for the short-range interaction they are always small enough compared to the cutoff wave vector to be taken equal to zero. This is clear when we notice that the two-body propagator is directly related to the one-body propagator corresponding to the relative motion of the electron and the hole. Let us call $\mathbf{k}_e$ and $\omega_e$ the entering wave vector and energy of the electron, with similarly $\mathbf{k}_h$ and $\omega_h$ for the entering hole, together with $\mathbf{k}'_e$ and $\omega'_e$, and similarly $\mathbf{k}'_h$ and $\omega'_h$ for the outgoing electron and hole, respectively. From momentum and energy conservation we have $\mathbf{P} = \mathbf{k}_e + \mathbf{k}_h = \mathbf{k}'_e + \mathbf{k}'_h$ and $\Omega = \omega_e + \omega_h = \omega'_e + \omega'_h$. For the relative motion the relevant

energy is again $\Omega_r = \Omega - \mathbf{P}^2/2M$, while the entering momentum is $\mathbf{k} = (m_h\mathbf{k}_e - m_e\mathbf{k}_h)/M = \mathbf{k}_e - (m_e/M)\mathbf{P}$ and the outgoing momentum is $\mathbf{k}' = (m_h\mathbf{k}'_e - m_e\mathbf{k}'_h)/M = \mathbf{k}'_e - (m_e/M)\mathbf{P}$. For our purpose the relative motion problem is solved as soon as we have the corresponding Green's function $G(\omega, \mathbf{k}, \mathbf{k}')$ defined by

$$G(\omega, \mathbf{k}, \mathbf{k}') = \langle \mathbf{k} | \frac{1}{\omega - H} | \mathbf{k}' \rangle, \quad (4)$$

where $H = \mathbf{p}^2/2\mu + V(\mathbf{r})$ is the Hamiltonian corresponding to the relative motion. This Green's function is basically obtained by solving the Schrödinger equation for $H$.

Formally, this Green's function satisfies the equation

$$G = g + gVG, \quad (5)$$

where $g(\omega, \mathbf{p}) = 1/(\omega - \mathbf{p}^2/2\mu + i0_+)$ is the free-particle Green's function for the relative motion. Since our $T_2$ is precisely the $T$ matrix for this relative motion, which satisfies

$$T_2 = V + VgT_2, \quad (6)$$

it is easily obtained from $G$, since it is readily checked that if $G$ satisfies Eq. (5), then $T_2 = g^{-1}Gg^{-1} - g^{-1}$ satisfies Eq. (6). This gives explicitly

$$T_2(\omega, \mathbf{k}, \mathbf{k}') = \left(\omega - \frac{\mathbf{k}^2}{2\mu}\right) [G(\omega, \mathbf{k}, \mathbf{k}') - g(\omega, \mathbf{k})\delta_{\mathbf{k},\mathbf{k}'}] \times \left(\omega - \frac{\mathbf{k}'^2}{2\mu}\right). \quad (7)$$

Hence we see that, while for short-range interaction $T_2$ depends only on the four-vector $P = \{\Omega, \mathbf{P}\}$, we now have to also take into account the dependence on the incoming $\mathbf{k}$ and outgoing $\mathbf{k}'$ wave vectors for the relative motion. Accordingly, we denote it $T_2(P; \mathbf{k}, \mathbf{k}')$, the total energy $\Omega$ entering only, as above, through the combination $\Omega_r = \Omega - \mathbf{P}^2/2M$.

We note that, for a general interaction potential $V(\mathbf{r})$, which might even be known only numerically, Eq. (6) for $T_2$ reads explicitly

$$T_2(\omega, \mathbf{k}, \mathbf{k}') = V(\mathbf{k} - \mathbf{k}') + \int \frac{d\mathbf{k}''}{(2\pi)^3} \frac{V(\mathbf{k} - \mathbf{k}'')T_2(\omega, \mathbf{k}'', \mathbf{k}')}{\omega - \mathbf{k}''^2/2\mu + i0_+}, \quad (8)$$

where $V(\mathbf{q})$ is the Fourier transform of $V(\mathbf{r})$. This integral equation can be solved numerically once for all for the relevant range of variables, and the result for $T_2(\omega, \mathbf{k}, \mathbf{k}')$ can be stored for further use in the equations we obtain below.
The complication of now having for $T_2$ a wave vector dependence implies a corresponding complication for $T_3$. Previously, for short-range interaction we could specify only the total four-vector of the exciton. Now we also have to specify the wave vectors of the electron and the hole making up the exciton after its breakup since they will enter in the $T_2$ describing all the further processes. In practice, it is not more complicated to start by specifying the four-vector corresponding to these two particles.

To be more specific, we now specialize to the case relevant for our problem of finding the ground-state energy of our three-body system. Since for the actual cases we consider it will clearly be found for the two electrons having opposite spins, we restrict ourselves to this case. Accordingly, in contrast with the preceding section where the two electrons were indistinguishable, they will now be distinguishable particles and, in particular, we will have no exchange processes. This leads us to introduce two scattering vertices for the electron and the exciton, instead of a single one as in the preceding section. We call $T_{3↑}$ the vertex corresponding to an electron $↑$ scattering on an exciton made of a hole and a $↓$ spin electron. Similarly, we introduce $T_{3↓}$ corresponding to the situation where the electron spins are exchanged. We recall that, in $T_{3↑}$, since all the hole-$↓$ electron interactions have been taken into account in the incoming exciton propagator, the only possible interaction of the hole after the exciton breakup is with the $↑$ electron. However, there is another possibility, namely, that the two electrons scatter. This possible process did not enter in the preceding section for cold gases because, at very low energy, the scattering is dominantly $s$ wave, which is forbidden for two identical fermions. In contrast, we naturally want to take into account here the Coulomb repulsion between electrons.

Accordingly, we have to consider the $T$ matrix corresponding to the sum of the repeated scattering between the two electrons, which is analogous to $T_2(P; k, k')$ except that the interaction is now repulsive and there is naturally no bound state. We denote this matrix by $T_2^e(P; k, k')$. Correspondingly, we have to consider a $T_3$ vertex where the

entering electrons have just interacted repeatedly so that they can no longer interact and the first interaction to be considered is between the hole and one of these electrons. We denote this vertex $T_{3h}$. Ultimately, the three vertices $T_{3↑}$, $T_{3↓}$, and $T_{3h}$ are a way to describe all the possible scattering between the hole and the two electrons. On the other hand, we have to consider for these three vertices only the case where the outgoing particles are the $↑$ electron and the exciton (just as the entering particles in $T_{3↑}$). As above, we denote by $P$ the total energy-momentum four-vector, and by $\bar{p}$ the energy-momentum of the outgoing $↑$ electron. Actually, just as above, we finally take $\bar{p} = 0$ and $P = \{-E, 0\}$. So, for simplicity, we do not indicate these variables in the $T_3$ vertices. On the other hand, as we mention above, we have to indicate the energy-momentum of the three entering particles. However, since their sum is $P$, we need only to write it for two particles and we choose to write the variables for the two electrons. The first variable is for the $↑$ electron and the second one for the $↓$ electron. Hence, we have the three vertices $T_{3↑}(p, p')$, $T_{3↓}(p, p')$, and $T_{3h}(p, p')$. In order to have notations similar to the ones used in Sec. II, we include, for example, in the definition of $T_{3↑}(p, p')$ the free hole and free down electron propagators corresponding to the broken exciton, and similarly for $T_{3↓}(p, p')$ and $T_{3h}(p, p')$.

Proceeding as for Eq. (2) we can now write integral equations relating these vertices. Let us start with $T_{3↑}(p, p')$. In contrast to the first term of Eq. (2) there is no exchange term since the electrons have opposite spins. Hence, one possibility is that the hole interacts with the $↑$ electron, which is described by $T_2$, and then any process may happen, which is described by $T_{3↓}$. This is shown diagrammatically in Fig. 3(a). This is completely analogous to the second term of Eq. (2). However, another possibility is that the two electrons interact, as described by $T_2^e$, followed by all the processes described by $T_{3h}$. This is represented diagrammatically in Fig. 3(b). The corresponding equation for $T_{3↑}(p, p')$ gathering these two possible kinds of processes reads as

<figure><img src="image_4.png" /><figcaption>FIG. 3. The two contributions to the integral equation Eq. (9) for $T_{3↑}(p, p')$. (a) The $↑$ electron has its first interaction with the hole coming from the exciton. (b) The $↑$ electron has its first interaction with the $↓$ electron coming from the exciton. Same notations as in Fig. 2. The shaded circle indicates that the hole and the $↓$ electron are coming from the exciton, and that the first interaction in the diagram should not be between them. $T_{3↓}$ and $T_{3h}$ are the corresponding vertices given respectively by Eqs. (10) and (11). Note that there is no change of sign due to fermion line crossings, the ones appearing in the figures are for readability and can be removed by deforming appropriately the propagator lines.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 4. The two contributions to the integral equation Eq. (10) for $T_{3↓}(p, p')$. (a) The ↓ electron has its first interaction with the hole coming from the exciton. (b) The ↓ electron has its first interaction with the ↑ electron coming from the exciton. Same notations as in Fig. 2. The shaded circle indicates that the hole and the ↑ electron are coming from the exciton, and that the first interaction in the diagram should not be between them. $T_{3↑}$ and $T_{3↓}$ are the corresponding vertices given respectively by Eqs. (9) and (11).</figcaption></figure>

$$T_{3\uparrow}(p, p') = g_e(p')g_h(P-p-p') \sum_k \left[ T_2(P-p'; \mathbf{p}-r(\mathbf{P}-\mathbf{p}'), \mathbf{k}-r(\mathbf{P}-\mathbf{p}')) T_{3\downarrow}(k, p') \right. \\ \left. + T_2^e \left( p+p'; \mathbf{p}-\frac{\mathbf{p}+\mathbf{p}'}{2}, \mathbf{k}-\frac{\mathbf{p}+\mathbf{p}'}{2} \right) T_{3h}(k, p+p'-k) \right], \quad (9)$$

where we set $r = m_e/M$. For the case of electron-electron scattering, this ratio becomes merely 1/2.

Similarly, for $T_{3↓}(p, p')$ we have the possibility that the hole interacts with the ↓ electron, described by $T_2$, followed by all the processes corresponding to $T_{3↑}(p, p')$, as shown in Fig. 4(a). There is also the possibility shown in Fig. 4(b) of having the two electrons interacting as described by $T_2^e$, followed by $T_{3h}$ processes.

However, there is finally a process, not possible for $T_{3↑}(p, p')$, which is merely that the incoming exciton breaks to produce the outgoing ↑ electron while the incoming ↓ electron forms with the hole the outgoing dimer. Since some involved propagators are already factored out, this gives just an additional term $\delta_{p,0}g_h(P-p-p')$ if we take into account $\bar{p} = 0$. This leads to the equation

$$T_{3\downarrow}(p, p') = \delta_{p,0}g_h(P-p-p') + g_e(p)g_h(P-p-p') \sum_k \left[ T_2(P-p; \mathbf{p}'-r(\mathbf{P}-\mathbf{p}'), \mathbf{k}-r(\mathbf{P}-\mathbf{p}')) T_{3\uparrow}(p, k) \right. \\ \left. + T_2^e \left( p+p'; \mathbf{p}-\frac{\mathbf{p}+\mathbf{p}'}{2}, \mathbf{k}-\frac{\mathbf{p}+\mathbf{p}'}{2} \right) T_{3h}(k, p+p'-k) \right]. \quad (10)$$

Finally, we have to write a similar equation for $T_{3h}(p, p')$. We have the possibilities that either one of the electrons interacts with the hole, which is described by $T_2$, followed by either $T_{3↑}(p, p')$ or $T_{3↓}(p, p')$, as shown in Figs. 5(a) and 5(b). But there is again the simple case where the ↓ electron forms the outgoing exciton with the hole, the remaining ↑ electron giving the outgoing electron. This leads as above to

$$T_{3h}(p, p') = \delta_{p,0}g_e(p') + g_e(p)g_e(p') \sum_k \left[ T_2(P-p; \mathbf{p}'-r(\mathbf{P}-\mathbf{p}'), \mathbf{k}-r(\mathbf{P}-\mathbf{p}')) T_{3\uparrow}(p, k) \right. \\ \left. + T_2(P-p'; \mathbf{p}-r(\mathbf{P}-\mathbf{p}'), \mathbf{k}-r(\mathbf{P}-\mathbf{p}')) T_{3\downarrow}(k, p') \right]. \quad (11)$$

<figure><img src="image_6.png" /><figcaption>FIG. 5. The two contributions to the integral equation Eq. (11) for $T_{3h}(p, p')$. (a) The hole has its first interaction with the ↑ electron. (b) The hole has its first interaction with the ↓ electron. Same notations as in Fig. 2. The shaded circle indicates that the two electrons had repeated interaction, and that the first interaction in the diagram should not be between them. $T_{3↑}$ and $T_{3↓}$ are the corresponding vertices given respectively by Eqs. (9) and (10).</figcaption></figure>
There is in our problem an additional symmetry which we have not used, namely, the fact that the Hamiltonian is invariant with respect to the exchange of the two electrons. However, it is apparent in our equations. We see that, if we write the equation for the combination $T_{3+}(p, p') \equiv T_{3\uparrow}(p, p') + T_{3\downarrow}(p', p)$, it depends only on

$T_{3+}$ in the right-hand side. Similarly, the equation for $T_{3-}(p, p') \equiv T_{3\uparrow}(p, p') - T_{3\downarrow}(p', p)$ depends only on $T_{3-}$. With respect to the $T_{3h}$ contribution, we can put it under a more convenient form by making first, in Eqs. (9) and (10), the change of variable $k \to k + (p+p')/2$, and then, in the term coming from Eq. (10), the change $k \to -k$. This leads to

$$ T_{3\pm}(p, p') = \pm \delta_{p',0} g_h(P-p-p') + g_e(p') g_h(P-p-p') \sum_k \left[ \pm T_2(P-p'; \mathbf{p}+r\mathbf{p}', \mathbf{k}+r\mathbf{p}') T_{3\pm}(p', k) \right. \\ \left. + T_2^e \left( p+p'; \frac{\mathbf{p}-\mathbf{p}'}{2}, \mathbf{k} \right) T_{3\pm} \left( \frac{p+p'}{2} + k, \frac{p+p'}{2} - k \right) \right], \quad (12) $$

where we set $T_{3h\pm}(p, p') \equiv T_{3h}(p, p') \pm T_{3h}(p', p)$. We also use $T_2^e\{p+p'; -[(\mathbf{p}-\mathbf{p}')/2], -\mathbf{k}\} = T_2^e\{p+p'; [(\mathbf{p}-\mathbf{p}')/2], \mathbf{k}\}$ valid for a potential satisfying $V(-\mathbf{r}) = V(\mathbf{r})$. Finally, we use the fact that we restrict ourselves to the case $\mathbf{P}=\mathbf{0}$ to simplify the equation.

To close our set of equations we need an equation for $T_{3h\pm}(p, p')$ which is readily obtained from Eq. (11). We find

$$ T_{3h\pm}(p, p') = \delta_{p,0} g_e(p') + g_e(p) g_e(p') \sum_k T_2(P-p; \mathbf{p}' + r\mathbf{p}, \mathbf{k} + r\mathbf{p}) T_{3\pm}(p, k) \pm (p \leftrightarrow p'). \quad (13) $$

Hence, we simplify our problem from three coupled equations for the three vertices $T_{3\uparrow}$, $T_{3\downarrow}$, and $T_{3h}$ to two coupled equations for the vertices $T_{3\pm}$ and $T_{3h\pm}$. Naturally we can even introduce the expression for $T_{3h\pm}(p, p')$ into

the equation for $T_{3\pm}$ to obtain a single equation, but this is not particularly convenient. Clearly, we find the ground state in the symmetric subset $T_{3+}$ and $T_{3h+}$, and accordingly we restrict ourselves for simplicity to this case in the following, although the equations for the antisymmetric subset can be similarly obtained with a few changes of sign.

A further simplification appears in the equations when we notice that in the right-hand side of Eqs. (12) and (13) appear only the vertices summed over the frequency component $\omega_k$ of the variable $k$, namely, only $\sum_{\omega_k} T_{3+}(p', k)$ and $\sum_{\omega_k} T_{3h+}[(p+p')/2+k, (p+p')/2-k]$. This could be expected since $T_2$ depends separately on the wave vectors of the incoming and outgoing particles, but only on their total energy. Introducing the new functions $T_3(p, \mathbf{k}) = \sum_{\omega_k} T_{3+}(p, k)$ and $S_3(p, \mathbf{k}) = \sum_{\omega_k} T_{3h+}(p+k, p-k) = \sum_{\omega_k} T_{3h+}(p-k, p+k)$, we obtain the following equation for them by summing Eq. (12) over the frequency $\omega_{p'}$:

$$ T_3(p, \mathbf{p}') = \delta_{\mathbf{p}',0} g_h(P-p) + \sum_{\omega_{\mathbf{p}'}} g_e(p') g_h(P-p-p') \sum_{\mathbf{k}} \left[ T_2(P-p'; \mathbf{p}+r\mathbf{p}', \mathbf{k}+r\mathbf{p}') T_3(p', \mathbf{k}) \right. \\ \left. + T_2^e \left( p+p'; \frac{\mathbf{p}-\mathbf{p}'}{2}, \mathbf{k} \right) S_3 \left( \frac{p+p'}{2}, \mathbf{k} \right) \right]. \quad (14) $$

Furthermore, setting in Eq. (13) $p=Q-q$ and $p'=Q+q$ and summing over the frequency $\omega_q$, we obtain for $S_3$:

$$ S_3(Q, \mathbf{q}) = \delta_{Q,\mathbf{q}} g_e(2Q) + \sum_{\omega_q} g_e(Q-q) g_e(Q+q) \sum_{\mathbf{k}} T_2(P-Q+q; \mathbf{Q}+\mathbf{q}+r(\mathbf{Q}-\mathbf{q}), \mathbf{k}+r(\mathbf{Q}-\mathbf{q})) T_3(Q-q, \mathbf{k}) \\ + (\mathbf{q} \leftrightarrow -\mathbf{q}). \quad (15) $$

We finally proceed to perform the frequencies integration. As in Sec. II, this relies on the analytical properties of the various involved quantities. For example, one can show that, in the second term in the rhs of Eq. (14), the $\omega_{p'}$ integration reduces to an evaluation of the integrand at the pole $\omega_{p'} = \mathbf{p}'^2/2m_e$ of $g_e(p')$. Hence, the result depends only on the on-the-shell value of $T_3$, namely, $T(\mathbf{p}, \mathbf{p}') \equiv T_3(\{\mathbf{p}^2/2m_e, \mathbf{p}\}, \mathbf{p}')$.

Similar results hold for the third term in the rhs of Eq. (14), and for Eq. (15). This leads us to introduce similarly the on-the-shell value $S(\mathbf{Q}, \mathbf{q}) \equiv S_3(\{-E-\mathbf{Q}^2/2m_e, \mathbf{Q}\}/2, \mathbf{q})$. We give the details of this frequencies integration in Appendix A, because they are somewhat complicated by change of variables. This leads us to equations relating only $T(\mathbf{p}, \mathbf{p}')$ and $S(\mathbf{Q}, \mathbf{q})$. We obtain
$$ T(\mathbf{p}, \mathbf{p}') = -\frac{\delta_{\mathbf{p}', \mathbf{0}}}{E + \frac{\mathbf{p}^2}{2\mu} - \frac{1}{E + \frac{\mathbf{p}^2 + \mathbf{p}'^2}{2m_e} + \frac{(\mathbf{p} + \mathbf{p}')^2}{2m_h}} \times \sum_{\mathbf{k}} \left[ T_2 \left( \left\{ -E - \frac{\mathbf{p}'^2}{2m_e}, -\mathbf{p}' \right\}; \mathbf{p} + r\mathbf{p}', \mathbf{k} + r\mathbf{p}' \right) T(\mathbf{p}', \mathbf{k}) + T_2^e \left( \left\{ -E - \frac{(\mathbf{p} + \mathbf{p}')^2}{2m_h}, \mathbf{p} + \mathbf{p}' \right\}; \frac{\mathbf{p} - \mathbf{p}'}{2}, \mathbf{k} \right) S(\mathbf{p} + \mathbf{p}', \mathbf{k}) \right] \quad (16) $$

$$ S(\mathbf{Q}, \mathbf{q}) = -\frac{\delta_{\mathbf{Q}/2, \mathbf{q}}}{E + \frac{\mathbf{Q}^2}{2\mu} - \frac{1}{E + \frac{\mathbf{Q}^2}{2m_h} + \frac{\mathbf{Q}^2 + \mathbf{q}^2}{4m_e}} \times \sum_{\mathbf{k}} T_2 \left( \left\{ -E - \frac{(\mathbf{Q} + \mathbf{q})^2}{2m_e}, -\frac{(\mathbf{Q} + \mathbf{q})}{2} \right\} ; \frac{\mathbf{Q} - \mathbf{q}}{2} + r \frac{\mathbf{q}}{2} \right) T \left( \frac{\mathbf{Q} + \mathbf{q}}{2}, \mathbf{k} \right) + (\mathbf{q} \leftrightarrow -\mathbf{q}), \quad (17) $$

where we now use the explicit form of the propagators $g_e$ and $g_h$.

Although it does not look so simple, this set of equations is clearly the best we could hope for this problem. We have two vertices $T$ and $S$ instead of one because we take into account not only electron-hole interaction but also electron-electron interaction. Moreover, these quantities depend only on three variables, the modulus of each vector and the angle between them, which is expected since $T_2$ and $T_2^e$ depend on the entering and outgoing wave vectors in the general case we are dealing with.

It is interesting to see how these equations simplify to something similar to what we had in Sec. II when $T_2$ depends only on the total energy-momentum and its dependence on wave vectors can be neglected, as it is the case for the short-range interaction considered in Sec. II. Indeed, in this case we see that, in the right-hand sides of these equations, the summation over $\mathbf{k}$ introduces merely $\sum_{\mathbf{k}} T(\mathbf{p}', \mathbf{k})$, $\sum_{\mathbf{k}} S(\mathbf{p} + \mathbf{p}', \mathbf{k})$, and $\sum_{\mathbf{k}} T(\mathbf{Q}/2 + \mathbf{q}, \mathbf{k})$. Introducing $t(\mathbf{p}) = \sum_{\mathbf{k}} T(\mathbf{p}, \mathbf{k})$ and $s(\mathbf{Q}) = \sum_{\mathbf{k}} S(\mathbf{Q}, \mathbf{k})$, and summing Eqs. (16) and (17) over $\mathbf{p}'$ and $\mathbf{q}$, respectively, we see that we obtain a set of integral equations for $t(\mathbf{p})$ and $s(\mathbf{Q})$ which is quite simple, since $t(\mathbf{p})$ and $s(\mathbf{Q})$ depend actually only on the single variables $|\mathbf{p}|$ and $|\mathbf{Q}|$, respectively, and the angular integrations can be performed easily. We do not write them explicitly in the general case since we do not make use of them.

Let us now focus on the specific problem of this paper, namely, finding the ground-state energy of the three-body problem. This is obtained in the general case by making use of the fact that our three-body vertices diverge when the energy is equal to a bound-state energy, just in the same way as $T_2$ has poles when the energy is equal to a bound state energy. This implies that, when $E$ is equal to the ground-state energy, the homogeneous parts of Eqs. (16) and (17) have a solution. Hence, we are led to find the solutions of a set of homogeneous integral equations.

Before addressing this problem specifically in the case of the helium atom ground state, we close this section by

considering a very particular case which offers an interesting and convenient check of our method. This is the case where there is no electron-electron interaction, which means we take $T_2^e = 0$. This makes $S(\mathbf{Q}, \mathbf{q})$ irrelevant, and we have only to consider $T(\mathbf{p}, \mathbf{p}')$. Moreover, we consider the additional particular case where the hole mass is infinite $m_h \to \infty$, so that $\mu = m_e$ and $r = 0$. In this case, the ground-state energy is obvious. Indeed, since the hole mass is infinite it can be considered as a fixed impurity, and the electrons just feel the attractive potential of this impurity. Moreover, since they do not interact, we have just two independent one-body problems, one for each electron. Hence, the ground-state energy is merely the sum of the ground-state energy $E_0$ of each electron. The case of a short-range interaction is considered specifically in Appendix B. It is interesting because the integral equation resulting from Eq. (16) is a simple one-dimensional integral equation. Nevertheless, the solution corresponding to our physical situation is not at all trivial, but it does exist and corresponds to the correct ground-state energy.

In the case of the general interaction, we are again left with Eq. (16). Since $M = \infty$, the relative motion energy entering $T_2$ is $-(E + p'^2/2m_e) = -(2E_0 + p'^2/2m_e)$, and we have to make $\mu = m_e$, $m_h = \infty$, and $r = 0$ in Eq. (16). This leads us for the homogeneous integral equation to

$$ T(\mathbf{p}, \mathbf{p}') = -\frac{1}{E + \frac{\mathbf{p}^2 + \mathbf{p}'^2}{2m_e}} \sum_{\mathbf{k}} T_2 \left( -E - \frac{\mathbf{p}'^2}{2m_e}, \mathbf{p}, \mathbf{k} \right) T(\mathbf{p}', \mathbf{k}), \quad (18) $$

with the notation $T_2(\{\Omega, \mathbf{P}\}; \mathbf{p}, \mathbf{k}) \equiv T_2(\Omega_r, \mathbf{p}, \mathbf{k})$. As we detail in Appendix C, we find that

$$ T(\mathbf{p}, \mathbf{p}') = \frac{E_0 + \frac{\mathbf{p}^2}{2m_e}}{2E_0 + \frac{\mathbf{p}^2 + \mathbf{p}'^2}{2m_e}} \varphi_0(\mathbf{p}) \varphi_0(\mathbf{p}') \quad (19) $$

is solution, where $\varphi_0(\mathbf{p})$ is the single electron ground-state wave function corresponding to the binding energy $E_0$.
IV. WAVE FUNCTION

As we explain above, the ground-state energy of our three-body problem (or actually any eigenenergy) is obtained in the general case by requiring that the homogeneous parts of Eqs. (16) and (17) have a solution. Let us rewrite these homogeneous equations for clarity:

$$ T(\mathbf{p}, \mathbf{p}') = -\frac{1}{E + \frac{\mathbf{p}^2 + \mathbf{p}'^2}{2m_e} + \frac{(\mathbf{p} + \mathbf{p}')^2}{2m_h}} \times \sum_{\mathbf{k}} \left[ T_2\left(\left\{-E - \frac{\mathbf{p}'^2}{2m_e}, -\mathbf{p}'\right\}; \mathbf{p} + r\mathbf{p}', \mathbf{k} + r\mathbf{p}'\right) T(\mathbf{p}', \mathbf{k}) \right. \\ \left. + T_2^e\left(\left\{-E - \frac{(\mathbf{p} + \mathbf{p}')^2}{2m_h}, \mathbf{p} + \mathbf{p}'\right\}; \frac{\mathbf{p} - \mathbf{p}'}{2}, \mathbf{k}\right) S(\mathbf{p} + \mathbf{p}', \mathbf{k}) \right], \quad (20) $$

$$ S(\mathbf{Q}, \mathbf{q}) = -\frac{1}{E + \frac{\mathbf{Q}^2}{2m_h} + \frac{\mathbf{Q}^2 + 4\mathbf{q}^2}{4m_e}} \times \sum_{\mathbf{k}} T_2\left(\left\{-E - \frac{(\mathbf{Q} + \mathbf{q})^2}{2m_e}, -\frac{(\mathbf{Q} + \mathbf{q})}{2}\right\}; \frac{\mathbf{Q} - \mathbf{q}}{2} + r\left(\frac{\mathbf{Q} + \mathbf{q}}{2}\right), \mathbf{k} + r\left(\frac{\mathbf{Q} - \mathbf{q}}{2}\right)\right) \\ \times T\left(\frac{\mathbf{Q}}{2} + \mathbf{q}, \mathbf{k}\right) + (\mathbf{q} \leftrightarrow -\mathbf{q}). \quad (21) $$

Once we find the ground-state energy of our three-body problem, we may suspect that the corresponding wave function is related to the corresponding residue of our $T_3$ matrix. Actually, we have to handle carefully the frequency variables in order to find the proper relation between the wave function and the residue. Moreover, we have to take into account properly the fact that the matrices we use above do not correspond precisely to the full $T_3$ matrix.

Let us define a Green's function $\tilde{G}_3$ for the propagation of our three-body system. Precisely, we set

$$ \tilde{G}_3(t, \mathbf{k}_\uparrow, \mathbf{k}_\downarrow, \mathbf{k}_h, \mathbf{k}'_\uparrow, \mathbf{k}'_\downarrow, \mathbf{k}'_h) = -i \langle 0 | c_{\mathbf{k}_\uparrow}(t) c_{\mathbf{k}_\downarrow}(t) c_{\mathbf{k}_h}(t) c_{\mathbf{k}'_h}^\dagger(0) c_{\mathbf{k}'_\downarrow}^\dagger(0) c_{\mathbf{k}'_\uparrow}^\dagger(0) | 0 \rangle, \quad (22) $$

where the operators $c_{\mathbf{k}_\uparrow, \mathbf{k}_\downarrow, \mathbf{k}_h}(t)$ annihilate at time $t > 0$ the $\uparrow, \downarrow$ electrons and the $h$ hole created at time $t = 0$ by the operators $c_{\mathbf{k}_\uparrow, \mathbf{k}_\downarrow, \mathbf{k}_h, \mathbf{k}'_\uparrow, \mathbf{k}'_\downarrow}^\dagger(0)$ acting on vacuum $|0\rangle$. If we introduce the eigenstates $|n\rangle$ and eigenenergies $\mathcal{E}_n$ of the three-body Hamiltonian $H$, related by $H|n\rangle = \mathcal{E}_n|n\rangle$, we have for the Fourier transform $G_3(\omega)$ of $\tilde{G}_3(t)$:

$$ G_3(\omega, \mathbf{k}_\uparrow, \mathbf{k}_\downarrow, \mathbf{k}_h, \mathbf{k}'_\uparrow, \mathbf{k}'_\downarrow, \mathbf{k}'_h) = \sum_n \frac{\Phi_n(\mathbf{k}_\uparrow, \mathbf{k}_\downarrow) \Phi_n^*(\mathbf{k}'_\uparrow, \mathbf{k}'_\downarrow)}{\omega - \mathcal{E}_n + i0_+}, \quad (23) $$

where in the wave function $\Phi_n(\mathbf{k}_\uparrow, \mathbf{k}_\downarrow) = \langle 0 | c_{\mathbf{k}_\uparrow} c_{\mathbf{k}_\downarrow} c_{\mathbf{k}_h} | n \rangle$ we take into account that the total momentum of our three-body system is zero, which implies $\mathbf{k}_h = -(\mathbf{k}_\uparrow + \mathbf{k}_\downarrow)$. We see indeed that the residue corresponding to the pole $\mathcal{E}_n$ is the product of the wave functions $\Phi_n(\mathbf{k}_\uparrow, \mathbf{k}_\downarrow)\Phi_n^*(\mathbf{k}'_\uparrow, \mathbf{k}'_\downarrow)$.

On the other hand, we define the $T_3$ matrix with three different times, instead of a single one, by setting

$$ \begin{aligned} \tilde{T}_3(t_1, t_2, t_3, \mathbf{k}_\uparrow, \mathbf{k}_\downarrow, \mathbf{k}_h, \mathbf{k}'_\uparrow, \mathbf{k}'_\downarrow, \mathbf{k}'_h) \\ = (-i)^3 \langle 0 | c_{\mathbf{k}_\uparrow}(t_1) c_{\mathbf{k}_\downarrow}(t_2) c_{\mathbf{k}_h}(t_3) c_{\mathbf{k}'_h}^\dagger(0) c_{\mathbf{k}'_\downarrow}^\dagger(0) c_{\mathbf{k}'_\uparrow}^\dagger(0) | 0 \rangle. \end{aligned} \quad (24) $$

Its Fourier transform $T_3(\omega_1, \omega_2, \omega_3, \mathbf{k}_\uparrow, \mathbf{k}_\downarrow, \mathbf{k}_h, \mathbf{k}'_\uparrow, \mathbf{k}'_\downarrow, \mathbf{k}'_h)$ with respect to $t_1, t_2, t_3$ is basically the quantity we have dealt with from the beginning of the paper. Note, however, that for simplicity we do not try to properly describe the final state of our three-body problem, since when we consider the pole contribution at $\mathcal{E}_n$, the entering and the outgoing variables decouple completely. The variables we consider correspond actually to the entering variables of the scattering problem.

We recover $\tilde{G}_3(t)$ by taking $t_1 = t_2 = t_3 = t$ in $T_3(t_1, t_2, t_3)$. In Fourier transform this implies the relation

$$ \begin{aligned} G_3(\omega) &= \left(\frac{i}{2\pi}\right)^2 \int \left(\prod_i d\omega_i\right) \delta\left(\omega - \sum_i \omega_i\right) T_3(\omega_1, \omega_2, \omega_3) \\ &= \left(\frac{i}{2\pi}\right)^2 \int d\omega_1 d\omega_2 T_3(\omega_1, \omega_2, \omega - \omega_1 - \omega_2), \end{aligned} \quad (25) $$

where, for clarity, we do not write the momentum variables. This last relation is easily checked in the particular case where there is no interaction between the particles in our three-body system. One finds as expected that the three-body wave function is merely the product of the wave functions of each particle.

Switching back to our original variables, and taking again into account that the total momentum is zero, this means that we obtain the three-body wave function corresponding to the energy $\mathcal{E}_n$ as the residue at $\Omega = \mathcal{E}_n$ of
<figure><img src="image_7.png" /><figcaption>FIG. 6. The three contributions to the expression Eq. (27) for the connected part of the general vertex $\mathcal{T}_3$.</figcaption></figure>

$$ \left(\frac{i}{2\pi}\right)^2 \int d\omega_p d\omega_{p'} \mathcal{T}_3(p, p', P-p-p'), \quad (26) $$

with $P = \{\Omega, \mathbf{P}\}$ and $\mathbf{P} = \mathbf{0}$.

We now have to write $\mathcal{T}_3$ in terms of the vertices $T_{3\uparrow}$, $T_{3\downarrow}$, and $T_{3h}$, which we introduced in Sec. III. Naturally, in

doing this we disregard all the disconnected diagrams which do not contribute to the ground state we are looking for. In the definition of $T_{3\uparrow}$, $T_{3\downarrow}$, and $T_{3h}$, we had not included the free propagators of the incoming particles as well as the $T_2$ matrix corresponding to the first two interacting particles. We now have to write them explicitly to obtain $\mathcal{T}_3$, as it is shown in Fig. 6.

This leads to

$$ \begin{aligned} \mathcal{T}_3(p, p', P-p-p') = & g_e(p)g_e(p')g_h(P-p-p') \left[ \sum_{p'_1} T_2(P-p; \mathbf{p}' + r\mathbf{p}, \mathbf{p}_1' + r\mathbf{p}) T_{3\uparrow}(p, p_1') \right. \\ & \left. \qquad + \sum_{p_1} T_2(P-p'; \mathbf{p} + r\mathbf{p}', \mathbf{p}_1 + r\mathbf{p}') T_{3\downarrow}(p_1, p') \right. \\ & \left. \qquad + \sum_k T_2^e \left( p + p'; \mathbf{p} - \frac{1}{2}(\mathbf{p} + \mathbf{p}'), \mathbf{k} - \frac{1}{2}(\mathbf{p} + \mathbf{p}') \right) T_{3h}(k, p+p'-k) \right], \end{aligned} \quad (27) $$

where we take $\mathbf{P} = \mathbf{0}$ into account. We then introduced $T_{3\pm}(p, p') = T_{3\uparrow}(p, p') \pm T_{3\downarrow}(p', p)$ and we looked for a pole where $T_{3+}(p, p')$ diverges while $T_{3-}(p, p')$ does not. This implies $2T_{3\uparrow}(p, p') = 2T_{3\downarrow}(p', p) = T_{3+}(p, p')$ in the vicinity of the pole. Similarly, $T_{3h\pm}(p, p') = T_{3h}(p, p') \pm T_{3h}(p', p)$, where $T_{3h+}(p, p')$ diverges while $T_{3h-}(p, p')$ does not, which leads to $2T_{3h}(p, p') = 2T_{3h}(p', p) = T_{3h+}(p, p')$. This leads to

$$ \begin{aligned} 2\mathcal{T}_3(p, p', P-p-p') = & g_e(p)g_e(p')g_h(P-p-p') \sum_k \left[ T_2(P-p; \mathbf{p}' + r\mathbf{p}, \mathbf{k} + r\mathbf{p}) T_{3+}(p, k) \right. \\ & \left. \qquad + T_2(P-p'; \mathbf{p} + r\mathbf{p}', \mathbf{k} + r\mathbf{p}') T_{3+}(p', k) + T_2^e \left( p + p'; \frac{1}{2}(\mathbf{p} - \mathbf{p}'), \mathbf{k} \right) \right. \\ & \left. \qquad \times T_{3h} \left( \frac{1}{2}(p+p') + k, \frac{1}{2}(p+p') - k \right) \right], \end{aligned} \quad (28) $$

where we made the change $k \to k + (p+p')/2$ in the last term. We now make use of Eqs. (12) and (13), where the first term in the right-hand side is omitted since it is not divergent. When we calculate from these equations the combination $g_e(p)T_{3+}(p, p') + g_e(p')T_{3+}(p', p) + g_h(P-p-p')T_{3h+}(p, p')$, making appropriately the change of variable $k \to -k$ and using of $T_{3h+}(p, p') = T_{3h+}(p', p)$, we obtain twice the right-hand side of Eq. (28). This leads to

$$ \begin{aligned} 4\mathcal{T}_3(p, p', P-p-p') = & g_e(p)T_{3+}(p, p') + g_e(p')T_{3+}(p', p) \\ & + g_h(P-p-p')T_{3h+}(p, p'). \end{aligned} \quad (29) $$

To obtain the wave function, we still have from Eq. (26) to sum this quantity over $\omega_p$ and $\omega_{p'}$. However, from the definitions $T_3(p, \mathbf{k}) = \sum_{\omega_k} T_{3+}(p, k)$ and $S_3(p, \mathbf{k}) = \sum_{\omega_k} T_{3h+}(p+k, p-k) = \sum_{\omega_k} T_{3h+}(p-k, p+k)$, given above Eq. (14), we have

$$ \begin{aligned} \sum_{\omega_p, \omega_{p'}} g_e(p) T_{3+}(p, p') &= \sum_{\omega_p} g_e(p) T_3(p, \mathbf{p}'), \\ \sum_{\omega_p, \omega_{p'}} g_e(p') T_{3+}(p', p) &= \sum_{\omega_{p'}} g_e(p') T_3(p', \mathbf{p}), \end{aligned} \quad (30) $$

and
$$ \sum_{\omega_p, \omega'_p} g_h(P-p-p') T_{3h+}(p, p') = \sum_{\omega_q} g_h(P-q) S_3\left(\frac{q}{2}, \mathbf{k}\right), \quad (31) $$

where this last relation is obtained by going from the variables $p$ and $p'$ to the variables $q = p + p'$ and $k = (p - p')/2$, and using the definition of $S_3(q/2, \mathbf{k})$.

Hence, we are left with the calculation of

$$ \begin{aligned} & 4 \sum_{\omega_p, \omega'_p} T_3(p, p', P - p - p') \\ & = \sum_{\omega_p} g_e(p) T_3(p, \mathbf{p}') + \sum_{\omega'_p} g_e(p') T_3(p', \mathbf{p}) \\ & \quad + \sum_{\omega_q} g_h(P - q) S_3\left(\frac{q}{2}, \frac{\mathbf{p} - \mathbf{p}'}{2}\right). \end{aligned} \quad (32) $$

We have seen [see Eq. (A3)] that $T_3(p, \mathbf{p}')$ is analytical for $Im\omega_p < 0$. Going from summation to integration by $\sum_{\omega_p} \rightarrow i \int d\omega_p/(2\pi)$, and closing the integration contour on $\omega_p$ in the lower complex plane, only the pole from $g_e(p)$ contributes, and the result is just the on-the-shell value of $T_3(p, \mathbf{p}')$, namely, $T_3(\{\mathbf{p}^2/2m_e, \mathbf{p}\}, \mathbf{p}') = T(\mathbf{p}, \mathbf{p}')$. Similarly, the second term in Eq. (32) gives $T(\mathbf{p}', \mathbf{p})$. Finally [see Eq. (A4)], $S_3(q/2, (\mathbf{p} - \mathbf{p}')/2)$ is analytical for $Im\omega_q > 0$. Closing the integration contour in the $\omega_q$ upper complex plane, we obtain from the definition $S(\mathbf{Q}, \mathbf{q}) = S_3\{-E - \mathbf{Q}^2/2m_e, \mathbf{Q}\}/2, \mathbf{q}$ that the third term in Eq. (32) is $S(\mathbf{p} + \mathbf{p}', (\mathbf{p} - \mathbf{p}')/2)$.

Finally, we end up with the conclusion that the three-body bound-state wave function is given, from the solution $T(\mathbf{p}, \mathbf{p}')$ and $S(\mathbf{Q}, \mathbf{q})$ of Eqs. (20) and (21), by

$$ T(\mathbf{p}, \mathbf{p}') + T(\mathbf{p}', \mathbf{p}) + S\left(\mathbf{p} + \mathbf{p}', \frac{\mathbf{p} - \mathbf{p}'}{2}\right) \quad (33) $$

within a multiplicative constant, since the solutions $T$ and $S$ of the homogeneous equations are not normalized. This expression is, as expected, invariant under the exchange of $\mathbf{p}$ and $\mathbf{p}'$.

Note that, when we make the substitution $\mathbf{Q} = \mathbf{p} + \mathbf{p}'$ and $\mathbf{q} = (\mathbf{p} - \mathbf{p}')/2$ in Eq. (21) for $S(\mathbf{Q}, \mathbf{q})$, we obtain an expression which is identical to the $T_2$ term (that is the first term in the right-hand side) in Eq. (20) for $T(\mathbf{p}, \mathbf{p}') + T(\mathbf{p}', \mathbf{p})$. Hence, there is no need, in calculating the wave function Eq. (33), to evaluate Eq. (21) since this is already done when Eq. (20) is evaluated. In particular, if we come back to the case where the two electrons are non-interacting, in which case $T_2^e = 0$, the wave function from Eq. (33) is $2[T(\mathbf{p}, \mathbf{p}') + T(\mathbf{p}', \mathbf{p})]$. We obtain the same result (without the irrelevant factor of 2) if we argue that in this case $S$ is irrelevant from Eq. (20) and that only the $T$ terms should contribute in Eq. (33). From the explicit

solution Eq. (C2) which we have found in this case, we see that the ground-state wave function is just the product $\varphi_0(\mathbf{p})\varphi_0(\mathbf{p}')$ of the two single electron ground-state wave functions, as expected.

## V. COULOMB T MATRIX

We now come to our specific problem of handling the Coulomb potential. Since $T_2$ is directly linked to the Green's function, itself obtained by solving the Schrödinger equation, one expects to be able to write the Coulomb $T_2$ in terms of solutions of the Coulomb Schrödinger equation, namely hypergeometric functions. This can indeed be done, but this leads to expressions which are not so easy to handle numerically. On the other hand, it is obviously quite important to have a convenient expression for this $T_2$ in order to obtain a numerically efficient solution for our problem, which is one of our basic purpose. Fortunately, such an expression has been obtained by Schwinger [22]. We merely give here the results and our notations, which are necessary to carry out our helium calculation. However, since this Schwinger result is not so well known, we review briefly its derivation in Appendix D for completeness and give some additional details.

The Hamiltonian is

$$ H = -\frac{1}{2\mu} \Delta_{\mathbf{r}} - \epsilon_s \frac{Ze^2}{4\pi\epsilon r}, \quad (34) $$

where $e$ is the electronic charge and $\epsilon$ is the permittivity of the medium, be it vacuum or semiconductor. In case of the interelectronic repulsion, we have $\epsilon_s = -1$ and $Z = 1$, while for electron-hole attraction $\epsilon_s = 1$ and $Z$ depends on the “hole” charge, since we want to also consider the case of He where we have $Z = 2$ for the nucleus.

It is convenient to take half the Bohr radius $a_0 = 4\pi\epsilon/(2\mu Ze^2)$ as the unit of length. In the same way we take $1/(2\mu a_0^2)$ as the energy unit and we set $\omega = -\kappa^2/(2\mu a_0^2)$, where $\kappa$ will be real since we actually have to consider only negative values for $\omega$. Similarly, we express the wave vectors in terms of the unit $1/a_0$, which leads us to introduce reduced wave vectors by $\mathbf{k} = \mathbf{q}/a_0$, and so on for other wave vectors. The equation

$$ \sum_{\mathbf{k}''} \langle \mathbf{k} | (\omega - H) | \mathbf{k}'' \rangle G(\omega, \mathbf{k}'' | \mathbf{k}') = \delta_{\mathbf{k}, \mathbf{k}'} \quad (35) $$

for the Green's function becomes

$$ \begin{aligned} & - (\kappa^2 + \mathbf{q}^2)g(\kappa, \mathbf{q}, \mathbf{q}') + \epsilon_s \frac{1}{2\pi^2} \int d\mathbf{q}'' \frac{g(\kappa, \mathbf{q}'' | \mathbf{q}')}{(\mathbf{q} - \mathbf{q}'')^2} \\ & = (2\pi)^3 \delta(\mathbf{q} - \mathbf{q}'), \end{aligned} \quad (36) $$

where we introduce a reduced Green's function $g$ by $G(\omega, \mathbf{q}/a_0, \mathbf{q}'/a_0) = (2\mu a_0^5)g(\kappa, \mathbf{q}, \mathbf{q}')$, and we go from
discrete to continuous variables by $\sum_q \rightarrow \int dq/(2\pi)^3$ and $\delta_{q,q'} \rightarrow (2\pi)^3 \delta(q-q')$. Similarly, we introduce a reduced $T$ matrix by $T_2(\omega, q/a_0, q'/a_0) = (a_0/2\mu)t_2(\kappa, q, q')$. With these notations the expression for $t_2$, obtained from Schwinger's expression for the Green's function, reads

$$t_2(\kappa, \mathbf{q}, \mathbf{q}') = -\frac{4\pi\epsilon_s}{(\mathbf{q}-\mathbf{q}')^2} - \frac{2\pi}{\kappa} \frac{1}{(\mathbf{q}-\mathbf{q}')^2} I(\kappa, z), \quad (37)$$

where we set

$$\begin{aligned} I(\kappa, z) &= \int_0^1 du \frac{u^{-\epsilon_s/2\kappa}}{u+z(1-u)^2}, \\ z &= \frac{(\kappa^2 + \mathbf{q}^2)(\kappa^2 + \mathbf{q}'^2)}{4\kappa^2(\mathbf{q}-\mathbf{q}')^2}. \end{aligned} \quad (38)$$

The first term in Eq. (37) is merely the well-known Born approximation for $t_2$, the only term to survive in the limit $\kappa \to \infty$. Integrating by parts in $I(\kappa, z)$, Eq. (37) can also be conveniently rewritten as

$$t_2(\kappa, \mathbf{q}, \mathbf{q}') = -\frac{4\pi\epsilon_s z}{(\mathbf{q}-\mathbf{q}')^2} \int_0^1 du \frac{u^{-\epsilon_s/2\kappa}(1-u^2)}{[u+z(1-u)^2]^2}. \quad (39)$$

In our domain of interest we have $z \ge 1/4$, the minimum being reached when $\mathbf{q}$ and $\mathbf{q}'$ are antiparallel with $|\mathbf{q}| \cdot |\mathbf{q}'| = \kappa^2$. This implies that the poles $u_1$ and $u_2$ in

the integrand of $I(\kappa, z)$ are always complex conjugate. They are on the unit circle and they go from $u_1 = u_2 = -1$ for $z = 1/4$ to $u_1 = u_2 = 1$ for $z \to \infty$. Hence, they do not cause any problem in the numerical evaluation of $I(\kappa, z)$.

## VI. CASE OF THE HELIUM GROUND STATE

We now come to the explicit treatment of the helium atom ground state (with naturally only Coulomb interaction between particles retained in the Hamiltonian). Our basic purpose is to see how our method works in practice numerically and, in particular, to check numerically that it is exact. However, since the nucleus mass is very large compared to the electronic mass, we simplify a bit our practical task by taking it infinite, i.e., $m_h = \infty$. Since we only want to display an effective application of our method, this is an unimportant simplification. Naturally, when in a following paper we will consider the case of the trion in semiconductors, this simplification will be unacceptable since the hole mass is usually even lighter than the conduction band electronic mass. But it is easy to see that this does not bring in practice any sizable complication. The ground-state energy is known [4] in this case from variational calculations with an extremely high precision. Keeping the precision suitable for our purpose it is given by $E_0 = 2.903\,724$ a.u. $= 5.807\,448$ Ry.

Since $m_h = \infty$ implies $r = 0$ and $\mu = m_e$, Eqs. (20) and (21) simplify to

$$T(\mathbf{p}, \mathbf{p}') = -\frac{2m_e}{2m_e E + \mathbf{p}^2 + \mathbf{p}'^2} \sum_{\mathbf{k}} \left[ T_2\left(-E - \frac{\mathbf{p}'^2}{2m_e}, \mathbf{p}, \mathbf{k}\right) T(\mathbf{p}', \mathbf{k}) + T_2^e\left(-E - \frac{(\mathbf{p} + \mathbf{p}')^2}{4m_e}, \frac{\mathbf{p} - \mathbf{p}'}{2}, \mathbf{k}\right) S(\mathbf{p} + \mathbf{p}', \mathbf{k}) \right], \quad (40)$$

$$S(\mathbf{Q}, \mathbf{q}) = -\frac{4m_e}{4m_e E + \mathbf{Q}^2 + 4\mathbf{q}^2} \sum_{\mathbf{k}} T_2\left(-E - \frac{(\mathbf{Q} + \mathbf{q})^2}{2m_e}, \frac{\mathbf{Q} - \mathbf{q}}{2}, \mathbf{k}\right) T\left(\frac{\mathbf{Q}}{2} + \mathbf{q}, \mathbf{k}\right) + (\mathbf{q} \leftrightarrow -\mathbf{q}). \quad (41)$$

Here, we have already used explicitly the fact that $T_2(\{\Omega, \mathbf{P}\}; \mathbf{k}, \mathbf{k}') = T_2(\Omega - \mathbf{P}^2/2M, \mathbf{k}, \mathbf{k}')$ as indicated at the beginning of Sec. III. This yields $T_2(\{-E - \mathbf{p}'^2/2m_e, -\mathbf{p}'\}; \mathbf{p}, \mathbf{k}) = T_2(-E - \mathbf{p}'^2/2m_e, \mathbf{p}, \mathbf{k})$ since $M = m_e + m_h$ is infinite in this case. Similarly, $T_2^e(\{-E, \mathbf{p} + \mathbf{p}'\}; (\mathbf{p} - \mathbf{p}')/2, \mathbf{k}) = T_2(-E - (\mathbf{p} + \mathbf{p}')^2/4m_e, (\mathbf{p} - \mathbf{p}')/2, \mathbf{k})$, since here $M = m_e + m_h = 2m_e$.

We now make use of the expression found in Sec. V for the Coulomb $T_2(\omega, \mathbf{k}, \mathbf{k}')$. It is naturally quite convenient to use the same reduced units as in Sec. V, namely, take $a_0 = 4\pi\epsilon/(2\mu Ze^2)$ as unit of length, $1/(2\mu a_0^2)$ as energy unit, setting $\omega = -\kappa^2/(2\mu a_0^2)$, and express the wave vectors in terms of the unit $1/a_0$. These reduced units are for the electron-nucleus Coulomb problem, which means in our case that the reduced mass is merely the

electronic mass $\mu = m_e$, and $Z = 2$. However, we have also to use the solution of the electron-electron Coulomb problem to obtain $T_2^e$. For this case we have naturally to translate back the result of Sec. V in physical units, and then use the above reduced units to write the proper reduced expression. Finally, as in Sec. V, we set $T_2 = (a_0/2m_e)t_2$ and similarly $T_2^e = (a_0/2m_e)t_2^e$. In the following we use the same notations as above for the wave vectors, but they now have to be understood as being in reduced units. However, in order to avoid any confusion, we use small letters $t(\mathbf{p}, \mathbf{p}')$ and $s(\mathbf{Q}, \mathbf{q})$, instead of $T(\mathbf{p}, \mathbf{p}')$ and $S(\mathbf{Q}, \mathbf{q})$ to indicate that we work now with reduced units.

In this way, with $K^2 = 2m_e a_0^2 E$, Eqs. (40) and (41) become
$$t(\mathbf{p}, \mathbf{p}') = -\frac{1}{K^2 + \mathbf{p}^2 + \mathbf{p}'^2} \int \frac{d\mathbf{k}}{(2\pi)^3} \left[ t_2\left(\sqrt{K^2 + \mathbf{p}'^2}, \mathbf{p}, \mathbf{k}\right) t(\mathbf{p}', \mathbf{k}) + t_2^e\left(2Z\kappa_e, \frac{\mathbf{p}_-}{2}, \mathbf{k}\right) s(\mathbf{p}_+, \mathbf{k}) \right], \quad (42)$$

$$s(\mathbf{Q}, \mathbf{q}) = -\frac{2}{2K^2 + \mathbf{Q}^2 + 4\mathbf{q}^2} \int \frac{d\mathbf{k}}{(2\pi)^3} t_2\left(\sqrt{K^2 + \mathbf{Q}_+^2}, \mathbf{Q}_-, \mathbf{k}\right) t(\mathbf{Q}_+, \mathbf{k}) + (\mathbf{q} \leftrightarrow -\mathbf{q}), \quad (43)$$

where we use the abbreviations $\mathbf{p}_{\pm} = \mathbf{p} \pm \mathbf{p}', \mathbf{Q}_{\pm} = \mathbf{Q}/2 \pm \mathbf{q}$, and $\kappa_e = \sqrt{K^2/2 + \mathbf{p}_+^2/4}$, and we have from Sec. V:

$$t_2(\kappa, \mathbf{q}, \mathbf{q}') = -\frac{4\pi z}{(\mathbf{q} - \mathbf{q}')^2} \int_0^1 du \frac{u^{-1/2\kappa}(1 - u^2)}{[u + z(1 - u)^2]^2}, \qquad t_2^e(\kappa, \mathbf{q}, \mathbf{q}') = \frac{1}{Z} \frac{4\pi z_e}{(\mathbf{q} - \mathbf{q}')^2} \int_0^1 du \frac{u^{1/2\kappa}(1 - u^2)}{[u + z_e(1 - u)^2]^2}, \quad (44)$$

with, in $t_2(\kappa, \mathbf{q}, \mathbf{q}')$, $z = (\kappa^2 + \mathbf{q}^2)(\kappa^2 + \mathbf{q}'^2)/[4\kappa^2(\mathbf{q} - \mathbf{q}')^2]$, and, in $t_2^e(\kappa, \mathbf{q}, \mathbf{q}')$, $z_e = (\kappa_e^2 + \mathbf{q}^2)(\kappa_e^2 + \mathbf{q}'^2)/[4\kappa_e^2(\mathbf{q} - \mathbf{q}')^2]$, with $\kappa_e = \kappa/(2Z)$. One can show easily from Eqs. (42) and (43) themselves that $t(\mathbf{p}, \mathbf{p}')$ and $s(\mathbf{Q}, \mathbf{q})$ go very rapidly to zero when the modulus of any of the argument wave vectors go to infinity. We notice also that the calculation to be performed to obtain $s(\mathbf{Q}, \mathbf{q})$ in Eq. (43) is just the same as the one appearing for the first term in the rhs of Eq. (42) provided the substitution $\mathbf{p} \to \mathbf{Q}_-$ and $\mathbf{p}' \to \mathbf{Q}_+$ is made. The second term of Eq. (42) is also quite analogous to the first one.

From rotational invariance $t(\mathbf{p}, \mathbf{p}')$ depends only on the moduli $p$ and $p'$ of the two wave vectors, together with the angle $\alpha$ between them, so we may write $t(\mathbf{p}, \mathbf{p}') = t(p, p', \cos \alpha)$. Turning now to the practical evaluation of the first integral in the rhs of Eq. (42), we see that for fixed polar angle $\theta$ of $\mathbf{k}$ with respect to $\mathbf{p}'$ only the factor $t_2(\sqrt{K^2 + \mathbf{p}'^2}, \mathbf{p}, \mathbf{k})$ depends on the azimuthal angle $\varphi$ of $\mathbf{k}$ with respect to $\mathbf{p}'$. It turns out that this integration can be performed analytically from Eq. (44), so we are left with performing the $k$ and the $\theta$ integration numerically, with in addition the $u$ integration to be performed to obtain $t_2$ and $t_2^e$. The same point can be made for the second term in the rhs of Eq. (42) provided we replace $\mathbf{p}$ and $\mathbf{p}'$ by $\mathbf{p}_-$ and $\mathbf{p}_+$. Finally, we have noticed that the evaluation in Eq. (43) is related to the one in Eq. (42). Hence, in practice the integral equations Eqs. (42) and (43) are actually two-dimensional integral equations, which are fairly simple to handle numerically. The integration with respect to $\theta$ is quite conveniently performed with Gaussian integration, while an appropriate Simpson method is well suited for the $k$ integration. Some small problems linked to a refined treatment in the $u$ integration are considered in Appendix E.

Let us call $A$ the linear operator corresponding to the action of the rhs of Eqs. (42) and (43) on the two-dimensional column vector $(t, s)$. Solving Eqs. (42) and (43) is equivalent to finding an eigenvector of $A$ with the eigenvalue $\lambda = 1$. This problem has a solution only if $E$ corresponds to the energy of bound states of our three-body Hamiltonian. In particular, for very large $E$ (that is very large binding energy), which implies

large $\kappa$ and $\kappa_e$, $t_2$ reduces to the Born approximation and, from the prefactors in Eqs. (42) and (43), the operator $A$ goes to zero and all its eigenvalues are quite small. Hence, none of them can be equal to 1, and there is no state with very large binding energy, as expected. If we decrease $E$, the largest positive eigenvalue $\lambda_{\max}$ of $A$ will grow. When it reaches 1 we will have obtained the largest possible value for $E$ corresponding to an eigenstate. In other words, we will have the ground-state energy.

It is easy to obtain the largest eigenvalue of $A$ by applying iteratively $A$ to some convenient starting vector $(t_0, s_0)$. Indeed, iterating $n$ times is equivalent to applying the operator $A^n$ to $(t_0, s_0)$. But for large values of $n$, $A^n$ is dominated by its largest eigenvalue $\lambda_{\max}^n$ and is essentially equivalent to a projection on the corresponding eigenvector and multiplication by $\lambda_{\max}^n$. This allows us to conveniently identify $\lambda_{\max}$ and the corresponding eigenvector. Actually, this procedure works only if the spectrum of $A$ does not have nasty features, such as closely spaced largest and second largest eigenvalues, or large negative eigenvalues. Fortunately, we find that, in our case, this procedure works quite nicely. We find that in practice 20 iterations already give a satisfactory convergence for the precision we consider. For example, going up to 40 iterations does not bring any sizable change. It is also convenient, in order to find the ground-state energy, to start from the situation where the electron-electron interaction is zero (for which the answer is known) and crank it progressively to its actual value. In this way, at each stage, the range where the ground-state energy lies is fairly well known.

In practice we find that Gaussian integration is extremely efficient. Typically the precision for the ground-state energy increases exponentially with the number of Legendre polynomials used. For practical purposes, going up to $\ell = 5$ is already quite enough, although we use $\ell = 10$ in the results below. The limitation for precision comes mainly from the mesh we use for the $k$ integration. With our notation we have for the ground-state energy exactly $2K^2 = 1$ when the electrons are not interacting, while when they interact the known result [4] is $2K^2 = 0.72593$. Taking successively 10 points, 20 points, and 40 points for our $k$ mesh, we find numerically for the noninteracting electrons $2K^2 = 1.00742$,
$2K^2 = 1.00066$, and $2K^2 = 1.00007$, while for the interacting electrons we obtain $2K^2 = 0.73488$, $2K^2 = 0.72641$, and $2K^2 = 0.72604$. Hence, we see that we obtain in both cases the exact result with a precision which is typically $10^{-2}$, $10^{-3}$, and $10^{-4}$ successively. This is a quite clear numerical check that our method converges rapidly toward the exact result when the precision of the numerical procedure is increased. On the other hand, the calculation time increases markedly with improved precision. On our MacBook Pro computer, the 1% precision (which is in practice quite enough for the determination of the trion binding energy) takes typically 1 mn. We need 30 mn to obtain the $10^{-3}$ precision and 280 mn for $10^{-4}$. These times are for calculations starting without a priori information on the result and no refinement. Naturally they are markedly shortened as soon as some information from preceding calculations is used and/or the mesh is made more precise only near the end of the calculation when one looks for improved precision on the result.

Let us now turn to the wave function. We plot in Fig. 7, in the $(p, p')$ plane, the contour lines of the normalized

<figure><img src="image_8.png" /><figcaption>FIG. 7. Contour lines for the two-electron wave function. The wave function is normalized to 1 for $p = p' = 0$. The contour lines correspond to the respective values 0.9, 0.5, 0.1, and 0.01 for the wave function as indicated. The wave function depends also on the angle between **p** and **p'**, but the dependence is rather weak, and only the contour lines for the angles 0 (blue) and $\pi$ (red) are drawn, as indicated for 0.1 and 0.01. When the angle increases from 0 to $\pi$ the contour lines interpolate smoothly between these two limits. Note that the axes are for arctan $p$ and arctan $p'$, not $p$ and $p'$. The dashed lines correspond, for the same values 0.9, 0.5, 0.1, and 0.01 for the wave function, to the case where the two electrons are noninteracting, the result being obtained both analytically and numerically.</figcaption></figure>

wave function for values 0.9, 0.5, 0.1, and 0.01, the wave function being normalized to 1 for $p = p' = 0$. Naturally the wave function also depends on the angle between **p** and **p'**, but this dependence turns out to be rather weak. Hence, we plot only the results when **p** and **p'** are parallel and antiparallel. The correlation between the two electrons is seen to be stronger in this latter case. When the angle goes from 0 to $\pi$, the contour lines interpolate smoothly between these two limits. Naturally there is no angular dependence when **p** = 0 or **p'** = 0, so the angular dependence is strongest along the diagonal $p = p'$. For the sake of comparison we also plot the same contour lines in the case where the two electrons are noninteracting and uncorrelated. In this case the wave function is naturally known analytically, being the product of the two single electron wave functions. As expected, the electron-electron repulsion leads to an expansion of the wave function in direct space, and correspondingly to a contraction in $k$ space. This is indeed what we find. Finally, let us indicate that, in the case of noninteracting electrons, our numerical solution for the wave function is in excellent agreement with the known analytical expression.

## VII. CONCLUSION

In this paper, we present an exact general approach for the solution of the three-body problem for a general interaction, which happens to be simple and fast, and apply it to the case of the Coulomb interaction. Rather than starting with the Schrödinger equation for this problem, it rather makes use of a corresponding integral equation derived from the consideration of the scattering properties of the system, namely, when one body is scattered by the two-body system formed by the two other ones. In this way one makes full use of the solution of the two-body problem which appears through the corresponding $T_2$ matrix, and the interaction does not appear explicitly but only through this known $T_2$ matrix. We show that the frequencies can be eliminated and only on-the-shell evaluations of the involved vertices appear. When two bodies have the same mass and the same interaction with the third body, finding the ground state (or any bound state) of the three-body system amounts to finding for which energy two coupled three-dimensional linear integral equations have a solution. The wave function is directly obtained from the corresponding solution.

We apply this approach to the well-known helium atom ground-state problem, making use of the $T_2$ matrix for the Coulomb potential obtained by Schwinger. In this case, the linear integral equations turn out to be two dimensional. We obtain a perfect numerical agreement with the known result for the ground-state energy. We expect to apply this approach in the near future to other three-body problems of interest, and in particular to the trion problem in semiconductors.

Our method leads to integral equations instead of the standard partial differential equations found in the
Schrödinger equation. Hence, the mathematical formulation is qualitatively different and offers the possibility of new analyses of the problem, including analytical approaches. In a quantitative point of view, our method leads to a strong reduction in the numerical work, basically because we make use of the solution of the two-body problem by making use of $T_2$. Hence, a part of the work is already done from the start. We stress that this advantage applies for any shape of the interaction potential because if no convenient form for $T_2$ is available, it can be obtained by solving numerically the general equation Eq. (8), storing the results, and using them by appropriate interpolations. The method could even apply to the three-body problem between composite particles (provided naturally that there is appropriate knowledge of the interaction between these particles). Roughly speaking, our method leads typically to a division by a factor 2 of the number of degrees of freedom, which obviously is a huge simplification for the numerical work. In this way for a 2D problem one goes from 4 to 2 degrees of freedom. Similarly, one goes from 6 to 3 degrees of freedom for a

3D problem. However, the precise result naturally depends on the symmetries of the specific problem to be solved, and in our case of helium ground state the rotational invariance allows us to go down to 2 degrees of freedom. Although in our specific numerical application we assume an infinite mass for the helium nucleus, it is easy to see that taking general masses does not lead to appreciable complication. Finally, our method applies as well to excited states as to ground state. However, the numerical solution of the integral equations [Eqs. (20) and (21)] for the ground state may be somewhat easier, as we have seen in the helium case. In general, the standard numerical solution will be through discretization of these equations and use of the powerful numerical methods of linear algebra for their solution.

## APPENDIX A: FREQUENCY INTEGRATION

To make the equations more transparent, let us rewrite for a moment Eqs. (14) and (15) by displaying only the frequency variables. We have

$$T_3(\omega_p) = g_h(-E - \omega_p) + \sum_{\omega_{p'}} g_e(\omega_{p'}) g_h(-E - \omega_p - \omega_{p'}) \left[ T_2(-E - \omega_{p'}) T_3(\omega_{p'}) + T_2^e(\omega_p + \omega_{p'}) S_3 \left( \frac{\omega_p + \omega_{p'}}{2} \right) \right], \quad (A1)$$

$$S_3(\omega_Q) = g_e(2\omega_Q) + \sum_{\omega_q} g_e(\omega_Q - \omega_q) g_e(\omega_Q + \omega_q) T_2(-E - \omega_Q + \omega_q) T_3(\omega_Q - \omega_q) + (\mathbf{q} \leftrightarrow -\mathbf{q}). \quad (A2)$$

Making in the third term of the right-hand side of Eq. (A1) the change of variable $\omega_{p'} = \omega_x - \omega_p$, and in Eq. (A2) the change $\omega_q = \omega_Q - \omega_y$, we find

$$T_3(\omega_p) = g_h(-E - \omega_p) + \sum_{\omega_{p'}} g_e(\omega_{p'}) g_h(-E - \omega_p - \omega_{p'}) T_2(-E - \omega_{p'}) T_3(\omega_{p'}) \\ + \sum_{\omega_x} g_e(\omega_x - \omega_p) g_h(-E - \omega_x) T_2^e(\omega_x) S_3\left(\frac{\omega_x}{2}\right), \quad (A3)$$

$$S_3(\omega_Q) = g_e(2\omega_Q) + \sum_{\omega_y} g_e(\omega_y) g_e(2\omega_Q - \omega_y) T_2(-E - \omega_y) T_3(\omega_y) + (\mathbf{q} \leftrightarrow -\mathbf{q}), \quad (A4)$$

where naturally one should again understand $\sum_{\omega_x} \rightarrow i \int d\omega_x/(2\pi)$, and similarly for the frequency variables.

The second equation shows that $S_3(\omega_Q)$ is analytical in the upper $\omega_Q$ complex plane. On the other hand, $T_3(\omega_p)$ is analytical in the lower $\omega_p$ complex plane, since the three terms in the right-hand side of Eq. (A3) have this property. We can now make use of these properties to perform the frequency integration. In Eq. (A3), just as in Eq. (2), we see that in the second term in the right-hand side all the factors except $g_e(\omega_{p'})$ are analytical in the lower $\omega_{p'}$ complex plane. Closing the $\omega_{p'}$ integration contour by a semicircle at infinity in this half-plane, the only contribution in a residue integration

comes from the pole of $g_e(\omega_{p'})$ at $\omega_{p'} = \mathbf{p}'^2/2m_e$. Hence, we need in particular to evaluate $T_3(\omega_{p'})$ only on the shell. Proceeding in the same way in Eq. (A4), we see once again that only the on-the-shell value of $T_3(y)$ appears. Finally, for the integration of the third term of Eq. (A3), we can close the contour in the upper half $\omega_x$ complex plane where the only contribution comes from the pole of $g_h(-E - \omega_x)$ at $\omega_x = -E - \mathbf{x}^2/2m_h$, where $\mathbf{x} = \mathbf{p} + \mathbf{p}'$, the other factors being analytical functions. So only the value of $S_3(\omega_x/2)$ for this specific frequency is required.

Coming back to the full equations [Eqs. (14) and (15)], we perform the same change of variables as described
above and perform the frequency integration in the way we have indicated. Since only on-the-shell quantities come in, this leads us to write the equations for these quantities. We define $T(\mathbf{p}, \mathbf{p}') = T_3(\{\mathbf{p}^2/2m_e, \mathbf{p}\}, \mathbf{p}')$ and $S(\mathbf{Q}, \mathbf{q}) = S_3(\{-E - \mathbf{Q}^2/2m_e, \mathbf{Q}\}/2, \mathbf{q})$. This leads to the equations given in the text.

## APPENDIX B: ZERO ELECTRON-ELECTRON INTERACTION: SHORT-RANGE INTERACTION

As indicated in the main text, when the interaction is short-ranged, we are led to consider $t(\mathbf{p}) \equiv \sum_{\mathbf{k}} T(\mathbf{p}, \mathbf{k})$. We obtain easily from Eq. (16) an equation for $t(\mathbf{p})$. We find, now with the simpler notation $|\mathbf{p}| \equiv p$,

$$t(p) = -\frac{1}{E + \frac{p^2}{2m_e}} - \frac{1}{2\pi^2} \int_0^\infty dp' p'^2 \frac{t_2(-E - \frac{p'^2}{2m_e})}{E + \frac{p'^2+p^2}{2m_e}} t(p'). \quad (B1)$$

Here, we make explicit the fact that $T_2(\{\Omega, \mathbf{P}\})$ depends on the total energy-momentum $\Omega$ and $\mathbf{P}$ only through the relative motion energy $\Omega_r = \Omega - \mathbf{P}^2/2M$ by setting $T_2(\{\Omega, \mathbf{P}\}) \equiv t_2(\Omega_r)$. In our case, $M = \infty$, which leads to $T_2(\{-E - p'^2/2m_e, -\mathbf{p}'\}) = t_2(-E - p'^2/2m_e)$. For the short-range interaction, we have from Eq. (1) $t_2(\Omega) = (2\pi a/m_e)[a^{-1} - \sqrt{-2m_e\Omega}]^{-1}$.

We expect the ground-state energy to be the sum of the ground-state energy of each electron. In particular, for the short-range interaction, the ground-state energy is twice the energy $E_0 = 1/2m_e a^2$ of the bound state. So the homogeneous part of Eq. (B1) should have a solution for $E = 2E_0 = 1/m_e a^2$. We can make in this equation changes of function and variables appropriate to get rid of $m_e$ and $a$. But this is equivalent to taking $m_e$ and $a$ as units of mass and length. This leads us to conclude that the homogeneous integral equation,

$$t(p) = \frac{2}{\pi} \int_{0}^{\infty} dp' \frac{p'^2}{p^2 + p'^2 + 2} \frac{\sqrt{p'^2 + 2} + 1}{p'^2 + 1} t(p'), \quad (B2)$$

should have a solution. Although the physical problem and the corresponding ground-state energy are trivial, this is not the case for the corresponding integral equation Eq. (B2). Nevertheless, it is easily checked that this equation has the solution $t(p) = 1/(\sqrt{p^2 + 2} + 1)$, which we derive more systematically in the case of general interaction. Hence, we have checked that our method gives the correct ground-state energy for this very particular case. Let us just mention that this check can be extended to the case where the two electrons have different masses, where a similar but somewhat more involved solution can be found.

## APPENDIX C: ZERO ELECTRON-ELECTRON INTERACTION: GENERAL INTERACTION

In Eq. (18) the general expression of $T_2$ is obtained from Eq. (7) with $\omega = -(E + \mathbf{p}'^2/2m_e)$:

$$T_2(\omega, \mathbf{k}, \mathbf{k}') = \left(\omega - \frac{\mathbf{k}^2}{2m_e}\right) \sum_n \frac{\varphi_n(\mathbf{k})\varphi_n(\mathbf{k}')}{\omega - \mathcal{E}_n} \left(\omega - \frac{\mathbf{k}'^2}{2m_e}\right) - \left(\omega - \frac{\mathbf{k}^2}{2m_e}\right) \delta_{\mathbf{k}, \mathbf{k}'}, \quad (C1)$$

where we have expressed the Green's function in terms of the eigenfunctions and eigenenergies of the relative motion Hamiltonian $H\varphi_n(\mathbf{k}) = \mathcal{E}_n\varphi_n(\mathbf{k})$. By time-reversal invariance we can take the eigenfunctions as real. Since in our specific case $E = 2E_0$ only the ground state is involved, we may suspect that only the ground-state wave function $\varphi_0(\mathbf{k})$ appears in $T(\mathbf{p}, \mathbf{p}')$. Indeed, making use of the orthonormality of the eigenfunctions $\sum_{\mathbf{k}} \varphi_m(\mathbf{k})\varphi_n(\mathbf{k}) = \delta_{m,n}$, we find that

$$T(\mathbf{p}, \mathbf{p}') = \frac{E_0 + \frac{\mathbf{p}^2}{2m_e}}{2E_0 + \frac{\mathbf{p}^2 + \mathbf{p}'^2}{2m_e}} \varphi_0(\mathbf{p})\varphi_0(\mathbf{p}') \quad (C2)$$

is solution, as it is easily checked by carrying this expression into Eq. (18) with $E = 2E_0$.

Actually, for this case of a general interaction, we can extend the above argument to the case where the electrons are, respectively, in excited states $\varphi_1(\mathbf{p})$ and $\varphi_2(\mathbf{p})$, with energy $-E_1$ and $-E_2$. Since they do not interact we should find an excited bound state of our three-body system with energy $-(E_1 + E_2)$. This should also give rise to a divergence of our three-body vertex, so we should have a corresponding solution for Eq. (18) for $E = E_1 + E_2$. Indeed, we find that this equation has the solution

$$T(\mathbf{p}, \mathbf{p}') = \frac{(E_1 + \frac{\mathbf{p}^2}{2m_e})\varphi_1(\mathbf{p})\varphi_2(\mathbf{p}') + (E_2 + \frac{\mathbf{p}^2}{2m_e})\varphi_1(\mathbf{p}')\varphi_2(\mathbf{p})}{E_1 + E_2 + \frac{\mathbf{p}^2 + \mathbf{p}'^2}{2m_e}}, \quad (C3)$$

which generalizes Eq. (C2).

Finally coming back to the short-range interaction case, we can immediately have the expression for the eigenfunction of the single bound state by comparing the general expression Eq. (7) of $T_2$ in the vicinity of the pole $\omega = -1/(2m_e a^2) = -E_0$ with its specific expression given below Eq. (B1). This gives $\varphi_0(\mathbf{p}) = (8\pi/a)^{1/2}[p^2 + a^{-2}]^{-1}$. We can then make use of the general expression for the solution Eq. (C2) and of the definition $t(\mathbf{p}) = \sum_{\mathbf{k}} T(\mathbf{p}, \mathbf{k})$ to find the solution of Eq. (B2). In this way one recovers (for $a = 1, m_e = 1$)
the solution $t(\mathbf{p}) = 1/(1 + \sqrt{\mathbf{p}^2 + 2})$ already given below Eq. (B2).

## APPENDIX D: SCHWINGER GREEN'S FUNCTION FOR THE COULOMB POTENTIAL

The simplicity of the result is linked to the hidden symmetry of the Coulomb potential, which gives rise in classical mechanics to the existence of a special conserved quantity, the Lenz vector, and which has been used by Pauli in its operatorial solution of the Coulomb Schrödinger equation.

We rewrite the equation for the Green's function in reduced units:

$$ -(\kappa^2 + \mathbf{q}^2)g(\kappa, \mathbf{q}, \mathbf{q}') + \epsilon_s \frac{1}{2\pi^2} \int d\mathbf{q}'' \frac{g(\kappa, \mathbf{q}'', \mathbf{q}')}{(\mathbf{q} - \mathbf{q}'')^2} = (2\pi)^3 \delta(\mathbf{q} - \mathbf{q}'). \quad (\text{D1}) $$

We proceed now to a change of variables. Let $O$ be the origin of the three-dimensional $\mathbf{q}$ space and $Q$ the point with $OQ = \mathbf{q}$. We consider $Q$ as the stereographic projection of a point $X$ which is on a 4D-sphere with diameter $\kappa$, the 3D $\mathbf{q}$ space being tangent to this sphere at the origin $O$, and the pole $S$ of the stereographic projection being such that $OS$ is a diameter of the sphere; see Fig. 8. From elementary geometry we have $SX \cdot SQ = \kappa^2$. If $C$ is the center of the sphere, we set $CX = (\kappa/2)\xi$, where $\xi$ runs on a 4D sphere with unit radius $\xi^2 = 1$ and therefore carries the information on the orientation of $CX$. Our change of variables is from $\mathbf{q}$ to $\xi$.

The components of the 4D vector $\xi$ respectively parallel and perpendicular to the $\mathbf{q}$ plane are given by

$$ \xi_{\parallel} = \frac{2\kappa}{\kappa^2 + \mathbf{q}^2} \mathbf{q}, \quad \xi_{\perp} = \frac{\kappa^2 - \mathbf{q}^2}{\kappa^2 + \mathbf{q}^2}. \quad (\text{D2}) $$

If $Q$ and $Q'$ are two points in the $\mathbf{q}$ plane, with corresponding points $X$ and $X'$ on the sphere, we have in the triangle

<figure><img src="image_9.png" /><figcaption>FIG. 8. Schematic view of the change of variables from $\mathbf{q}$ to $\xi$.</figcaption></figure>

with sides $SXQ$ and $SX'Q'$, and common angle $\gamma$ between these two sides, $QQ'^2 = SQ^2 + SQ'^2 - 2SQ \cdot SQ' \cos \gamma$ and $XX'^2 = SX^2 + SX'^2 - 2SX \cdot SX' \cos \gamma$. Together with $SX \cdot SQ = SX' \cdot SQ' = \kappa^2$, this leads to $XX'^2 = \kappa^4 QQ'^2 / (SQ^2 \cdot SQ'^2)$, that is,

$$ (\xi - \xi')^2 = \frac{4\kappa^2}{(\kappa^2 + \mathbf{q}^2)(\kappa^2 + \mathbf{q}'^2)} (\mathbf{q} - \mathbf{q}')^2. \quad (\text{D3}) $$

Letting $\mathbf{q}' \to \mathbf{q}$ and $\xi' \to \xi$, this implies

$$ d\xi = \frac{2\kappa}{\kappa^2 + \mathbf{q}^2} dq \quad (\text{D4}) $$

for the corresponding infinitesimal lengths $dq$ and $d\xi$. Taking three such corresponding infinitesimal orthogonal variations for $\mathbf{q}$ and $\xi$, this leads to the relation between elementary volumes,

$$ d\xi = \left( \frac{2\kappa}{\kappa^2 + \mathbf{q}^2} \right)^3 d\mathbf{q}, \quad (\text{D5}) $$

where $d\xi$ is the elementary solid angle on the 4D unit sphere. This relation implies that the Jacobian of our transformation is $[2\kappa/(\kappa^2 + \mathbf{q}^2)]^3$. This has to be used in making the change of variables in the right-hand side of Eq. (D1):

$$ \delta(\mathbf{q} - \mathbf{q}') = \left( \frac{2\kappa}{\kappa^2 + \mathbf{q}^2} \right)^3 \delta(\xi - \xi'). \quad (\text{D6}) $$

Finally we make in Eq. (D1) the change

$$ g(\kappa, \mathbf{q}, \mathbf{q}') = - \frac{(2\pi)^3 (2\kappa)^3}{(\kappa^2 + \mathbf{q}^2)^2 (\kappa^2 + \mathbf{q}'^2)^2} \Gamma(\kappa, \xi, \xi'), \quad (\text{D7}) $$

which leads to

$$ \Gamma(\kappa, \xi, \xi') - \frac{\epsilon_s}{4\pi^2\kappa} \int d\xi'' \frac{\Gamma(\kappa, \xi'', \xi')}{(\xi - \xi'')^2} = \delta(\xi - \xi'). \quad (\text{D8}) $$

The integral is just a convolution product, and this equation is simply solved by an expansion in spherical harmonics. These are the 4D generalization [24] $Y_{n,j}(\xi)$ of the standard 3D spherical harmonics $Y_{lm}$. Actually, since all the quantities in Eq. (D8) depend only on the angle between the two involved directions, an expansion in the corresponding [24] 4D Legendre polynomials $P_n(\xi \cdot \xi')$ is enough. They are related to the spherical harmonics by
$$P_n(\xi \cdot \xi') = \frac{S_{D-1}}{N(D, n)} \sum_{j=1}^{N(D,n)} Y_{n,j}(\xi) Y_{n,j}^*(\xi'). \quad (D9)$$

Here, $S_{D-1}$ is the surface of the unit sphere in dimension $D$ ($S_2 = 4\pi$ and $S_3 = 2\pi^2$), $j$ is collectively for the azimuthal numbers necessary for the enumeration of the spherical harmonics with a given $n$, and $N(D, n) = [(2n + D - 2)/n]^{(n+D-3)}/(n-1)!$ is the degeneracy of the $n$ level, that is, the number of different spherical harmonics with a given $n$. For example, $N(3, n) = 2n + 1$ and $N(4, n) = (n + 1)^2$. The spherical harmonics are orthonormal:

$$\int d\xi Y_{n,j}^*(\xi) Y_{n',j'}(\xi') = \delta_{nn'} \delta_{jj'}, \quad (D10)$$

which implies

$$\int d\xi'' P_n(\xi \cdot \xi'') P_{n'}(\xi'' \cdot \xi') = \frac{S_{D-1}}{N(D, n)} P_n(\xi \cdot \xi') \delta_{nn'}. \quad (D11)$$

Moreover, they satisfy the closure relation:

$$\sum_{n,j} Y_{n,j}(\xi) Y_{n,j}^*(\xi') = \sum_{n} \frac{N(D, n)}{S_{D-1}} P_n(\xi \cdot \xi') = \delta(\xi - \xi'). \quad (D12)$$

Finally, the Legendre polynomials are linked to their generating function by

$$\frac{1}{(1 - 2rt + r^2)^{D/2 - 1}} = \sum_{n=0}^{\infty} \frac{D-2}{2n + D - 2} N(D, n) r^n P_n(t), \quad (D13)$$

which implies [24]

$$\frac{1 - r^2}{(1 - 2rt + r^2)^{D/2}} = \sum_{n=0}^{\infty} N(D, n) r^n P_n(t) \quad (D14)$$

obtained by multiplying Eq. (D13) by $r^{D/2-1}$ and taking the derivative with respect to $r$.

Taking $D=4$ and $r=1$ in Eq. (D13), together with $(\xi - \xi')^2 = 2(1 - \xi \cdot \xi')$, we have

$$\frac{1}{(\xi - \xi')^2} = \sum_{n=0}^{\infty} (n + 1) P_n(\xi \cdot \xi'). \quad (D15)$$

Inserting a Legendre polynomial expansion of $\Gamma(\kappa, \xi, \xi')$ in Eq. (D8),

$$\begin{aligned} \Gamma(\kappa, \xi, \xi') &= \sum_{n,j} \Gamma_n(\kappa) Y_{n,j}(\xi) Y_{n,j}^*(\xi') \\ &= \sum_n \frac{N(4, n)}{S_3} \Gamma_n(\kappa) P_n(\xi \cdot \xi'), \end{aligned} \quad (D16)$$

together with Eqs. (D15), (D11), and (D12), we obtain the solution of Eq. (D8) as

$$\Gamma_n(\kappa) = \frac{1}{1 - \frac{\epsilon_s}{2\kappa(n+1)}}. \quad (D17)$$

In particular, in the attractive case $\epsilon_s = 1$, the poles of $\Gamma_n(\kappa)$ are found for $2\kappa = 1/(n+1)$ with $n=0, 1, \dots$, which give the expected energies $\omega = -\kappa^2/(2\mu a_0^2) = -\mu(Ze^2)^2/[2(n+1)^2(4\pi\epsilon)^2]$ for the bound states, with the expected degeneracy $N(4, n) = (n+1)^2$.

It is possible to sum up the series Eq. (D16) for $\Gamma(\kappa, \xi, \xi')$ leading to a very nice closed form. We write

$$\begin{aligned} \Gamma_n(\kappa) &= \frac{1}{1 - \frac{\epsilon_s}{2\kappa(n+1)}} \\ &= 1 + \frac{\epsilon_s}{2\kappa(n+1)} + \frac{1}{(2\kappa)^2(n+1)} \frac{1}{n+1 - \frac{\epsilon_s}{2\kappa}}, \end{aligned} \quad (D18)$$

and make use in the last term of the integral representation:

$$\frac{1}{n + 1 - \frac{\epsilon_s}{2\kappa}} = \int_0^1 du \, u^{-\epsilon_s/2\kappa} u^n. \quad (D19)$$

Moreover, writing $1 - 2\xi \cdot \xi' u + u^2 = (1 - u)^2 + u(\xi - \xi')^2$, we use Eq. (D13) to sum up the corresponding series:

$$\sum_{n=0}^{\infty} \frac{u^n}{n+1} N(4, n) P_n(\xi \cdot \xi') = \frac{1}{(1-u)^2 + u(\xi - \xi')^2}. \quad (D20)$$

The first two terms in Eq. (D18) are summed through Eqs. (D12) and (D15). This leads us finally to

$$\begin{aligned} \Gamma(\kappa, \xi, \xi') = \delta(\xi - \xi') &+ \frac{\epsilon_s}{4\pi^2\kappa} \frac{1}{(\xi - \xi')^2} \\ &+ \frac{1}{8\pi^2\kappa^2} \int_0^1 du \frac{u^{-\epsilon_s/2\kappa}}{(1-u)^2 + u(\xi - \xi')^2}. \end{aligned} \quad (D21)$$

As pointed out by Schwinger, it is possible to change the integration contour to obtain, for any value of the reduced energy $\kappa$, a well-defined expression. However, in deforming the contour, care must be taken of the possible
contributions of the poles coming from the denominator. But we do not need to perform such a transformation.

Let us now come to the T matrix itself. It is convenient to also introduce reduced units in Eq. (7). Introducing, just as for the Green's function, a reduced T matrix by $T_2(\omega, \mathbf{q}/a_0, \mathbf{q}'/a_0) = (a_0/2\mu)t_2(\kappa, \mathbf{q}, \mathbf{q}'),$ Eq. (7) becomes

$$t_2(\kappa, \mathbf{q}, \mathbf{q}') = (\kappa^2 + \mathbf{q}^2)(\kappa^2 + \mathbf{q}'^2)[g(\kappa, \mathbf{q}, \mathbf{q}') - g_0(\kappa, \mathbf{q}, \mathbf{q}')], \quad (D22)$$

where $g_0(\kappa, \mathbf{q}, \mathbf{q}') = -(2\pi)^3\delta(\mathbf{q}-\mathbf{q}')/(\kappa^2 + \mathbf{q}^2)$. On the other hand, we rewrite Eq. (D21) in terms of the variable $\mathbf{q}$ by making use of Eqs. (D3), (D6), and (D7). As could be expected, the $\delta(\xi - \xi')$ term in Eq. (D21) cancels exactly with the $g_0$ term in Eq. (D22), and we get finally

$$t_2(\kappa, \mathbf{q}, \mathbf{q}') = -\frac{4\pi\epsilon_s}{(\mathbf{q}-\mathbf{q}')^2} - \frac{2\pi}{\kappa} \frac{1}{(\mathbf{q}-\mathbf{q}')^2} I(\kappa, z), \quad (D23)$$

where we have set

$$I(\kappa, z) = \int_0^1 du \frac{u^{-\epsilon_s/2\kappa}}{u+z(1-u)^2}, \quad z = \frac{(\kappa^2 + \mathbf{q}^2)(\kappa^2 + \mathbf{q}'^2)}{4\kappa^2(\mathbf{q}-\mathbf{q}')^2}. \quad (D24)$$

The first term in Eq. (D23) is merely the well-known Born approximation for $t_2$, the only term to survive in the limit $\kappa \to \infty$. Integrating by parts in $I(\kappa, z)$, Eq. (D23) can also be conveniently rewritten as

$$t_2(\kappa, \mathbf{q}, \mathbf{q}') = -\frac{4\pi\epsilon_s z}{(\mathbf{q}-\mathbf{q}')^2} \int_0^1 du \frac{u^{-\epsilon_s/2\kappa}(1-u^2)}{[u+z(1-u)^2]^2}. \quad (D25)$$

In our domain of interest we have $z \ge 1/4$, the minimum being reached when $\mathbf{q}$ and $\mathbf{q}'$ are antiparallel with $|\mathbf{q}| \cdot |\mathbf{q}'| = \kappa^2$. This implies that the poles $u_1$ and $u_2$ in the integrand of $I(\kappa, z)$ are always complex conjugate. They are on the unit circle and they go from $u_1 = u_2 = -1$ for $z = 1/4$ to $u_1 = u_2 = 1$ for $z \to \infty$. Hence, they do not cause any problem in the numerical evaluation of $I(\kappa, z)$. It is possible to express in general $I(\kappa, z)$ in terms of hypergeometric functions, but this is not particularly useful

for its numerical evaluation. A particular case is for $z = 1/4$, where one finds $I(\kappa, 1/4) = -2 - (\epsilon_s/\kappa) [\psi(1/2 - \epsilon_s/4\kappa) - \psi(-\epsilon_s/4\kappa)]$, with $\psi(x)$ being the digamma function.

In the case $\epsilon_s = 1$, the discrete spectrum for $1/2\kappa = n$ gives rise to poles in $t_2$, which appear from $I(\kappa, z)$ through the divergent behavior of the integrand for $u \to 0$. In particular, for $\kappa \to 1/2$, one finds easily $I(\kappa, z) \approx 1/[(2\kappa - 1)z]$, leading to $t_2(\kappa, \mathbf{q}, \mathbf{q}') \approx -4\pi/[(2\kappa - 1)(\mathbf{q}^2 + 1/4)(\mathbf{q}'^2 + 1/4)]$. This has to be compared to the expression of the Green's function, which in reduced units gives $g(\kappa, \mathbf{q}, \mathbf{q}') \approx t_2(\kappa, \mathbf{q}, \mathbf{q}')/[(\kappa^2 + \mathbf{q}^2)(\kappa^2 + \mathbf{q}'^2)] \approx 2\pi/[(\kappa^2 + 1/4)(\mathbf{q}^2 + 1/4)^2(\mathbf{q}'^2 + 1/4)^2]$. This agrees with the general expression $q_0(\mathbf{k})q_0(\mathbf{k}')/(\omega - \mathcal{E}_0)$ of the Green's function in the vicinity of this pole, which reads in reduced units $\Phi_0(\mathbf{q})\Phi_0(\mathbf{q}')/(-\kappa^2 + 1/4)$, and with the expression of the normalized ground-state wave function $\Phi_0(\mathbf{q}) = (2\pi)^{1/2}/(\mathbf{q}^2 + 1/4)^2$.

Let us finally note that it is possible to check analytically, with this expression of the ground-state wave function and the expression of $t_2$ given by Eq. (D23), that Eq. (C2) is indeed solution of Eq. (18), although the corresponding calculation is not that simple.

## APPENDIX E: DETAILS ON u INTEGRATION

There are two small problems arising in the $u$ integration, one for $u \to 0$ and the other for $u \to 1$. So it is better to split the integral into two integrals, in order to handle separately the $u \to 0$ and $u \to 1$ problems. The $u \to 0$ problem, which results from the somewhat singular behavior $u^{\mp 1/(2\kappa)}$, is easily settled by an appropriate change of variables. In the first term in the rhs of Eq. (42), the $u \to 1$ problem arises because, when $\mathbf{k} \to \mathbf{p}$, $z$ in Eq. (44) diverges and the $u$ integrand behaves as $1/(1-u)^2$, resulting in a divergent $u$ integral for $u \to 1$. The dominant behavior can be extracted and handled analytically. Ultimately this leads to an integrable logarithmic singularity at $k = p$ in the $k$ integration, which is nevertheless annoying numerically. We prefer to avoid this difficulty altogether by having an integrand which is exactly zero for $\mathbf{k} = \mathbf{p}$. This is done by subtracting and adding a same quantity in the integrand. Explicitly, the first term in the rhs of Eq. (42) becomes

$$\frac{1}{8\pi^2\kappa^2} \int d\Omega_k \int_0^\infty dk [k^2(k^2 + \kappa^2)t(p', k, \cos\theta) - p^2(p^2 + \kappa^2)t(p', p, \cos\alpha)] \int_0^1 du \frac{u^{-1/2\kappa}(1-u^2)}{[u(\mathbf{p}-\mathbf{k})^2 + z_0(1-u)^2]^2} \\ + \frac{1}{8\pi^2\kappa^2} p^2(p^2 + \kappa^2)t(p', p, \cos\alpha) \int d\Omega_k \int_0^\infty dk \int_0^1 du \frac{u^{-1/2\kappa}(1-u^2)}{[u(\mathbf{p}-\mathbf{k})^2 + z_0(1-u)^2]^2}, \quad (E1)$$

where, for simplicity and clarity, we do not write explicitly the analytical result of the $\varphi$ integration in the $d\Omega_k = \sin\theta d\theta d\varphi$ integration. We set $z_0 = (\kappa^2 + \mathbf{p}^2)(\kappa^2 + \mathbf{k}^2)/(4\kappa^2)$. In the second term, both the $k$ and the $\Omega_k$ integration can
be performed analytically, and one is left with only the u integration to be performed numerically. The second term in Eq. (42) and the rhs in Eq. (43) are handled in essentially the same way.

[1] D. Blume, *Few-Body Physics with Ultracold Atomic and Molecular Systems in Traps*, Rep. Prog. Phys. **75**, 046401 (2012).

[2] C. H. Greene, P. Giannakeas, and J. Perez-Rios, *Universal Few-Body Physics and Cluster Formation*, Rev. Mod. Phys. **89**, 035006 (2017).

[3] E. A. Hylleraas, *Neue Berechnung der Energie des Heliums im Grundzustande, Sowie des Tiefsten Terms von Ortho-Helium*, Z. Phys. **54**, 347 (1929).

[4] H. Nakashima and H. Nakatsuji, *Solving the Schrödinger Equation for Helium Atom and Its Isoelectronic Ions with the Free Iterative Complement Interaction (ICI) Method*, J. Chem. Phys. **127**, 224104 (2007).

[5] A. M. Frolov, *Bound-State Properties of Negatively Charged Hydrogenlike Ions*, Phys. Rev. A **58**, 4479 (1998).

[6] S. Chandrasekhar, *Some Remarks on the Negative Hydrogen Ion and Its Absorption Coefficient*, Astrophys. J. **100**, 176 (1944).

[7] M. A. Lampert, *Mobile and Immobile Effective-Mass-Particle Complexes in Nonmetallic Solids*, Phys. Rev. Lett. **1958** 450, 1.

[8] K. Kheng, R. T. Cox, Y. Merle dAubigné, F. Bassani, K. Saminadayar, and S. Tatarenko, *Observation of Negatively Charged Excitons X in Semiconductor Quantum Wells*, Phys. Rev. Lett. **71**, 1752 (1993).

[9] V. Huard, R. T. Cox, K. Saminadayar, A. Arnoult, and S. Tatarenko, *Bound States in Optical Absorption of Semiconductor Quantum Wells Containing a Two-Dimensional Electron Gas*, Phys. Rev. Lett. **84**, 187 (2000).

[10] F. Ferlaino, A. Zenesini, M. Berninger, B. Huang, H.-C. Nägerl, and R. Grimm, *Efimov Resonances in Ultracold Quantum Gases*, Few-Body Syst. **51**, 113 (2011).

[11] L. Pitaevskii and S. Stringari, *Bose-Einstein Condensation* (Oxford University Press, New York, 2003).

[12] S. Giorgini, L. P. Pitaevskii, and S. Stringari, *Theory of Ultracold Atomic Fermi Gases*, Rev. Mod. Phys. **80**, 1215 (2008).

[13] F. Werner and Y. Castin, *Unitary Quantum Three-Body Problem in a Harmonic Trap*, Phys. Rev. Lett. **97**, 150401 (2006).

[14] J. P. Kestner and L. M. Duan, *Level Crossing in the Three-Body Problem for Strongly Interacting Fermions in a Harmonic Trap*, Phys. Rev. A **76**, 033611 (2007).

[15] X. J. Liu, H. Hu, and P. D. Drummond, *Virial Expansion for a Strongly Correlated Fermi Gas*, Phys. Rev. Lett. **102**, 160401 (2009).

[16] X. J. Liu, H. Hu, and P. D. Drummond, *Three Attractively Interacting Fermions in a Harmonic Trap: Exact Solution, Ferromagnetism, and High-Temperature Thermodynamics*, Phys. Rev. A **82**, 023619 (2010). Together with the preceding works [13–15] this paper makes a very effective use of the fact that, in an isotropic harmonic trap, the one-particle eigenstates form a discrete set. Together with a partial analytical treatment this allows a quite efficient numerical solution of the problem. It would be interesting to apply to this specific case our approach of making full use of the *T* matrix of the two-body problem. We expect it could bring additional progress. However, the contact interaction, the isotropy, and the discreteness of the spectrum make the solution found by these papers so efficient that we do not expect an important quantitative progress.

[17] G. V. Skorniakov and K. A. Ter-Martirosian, *Three Body Problem for Short Range Forces. I. Scattering of Low Energy Neutrons by Deuterons*, Zh. Eksp. Teor. Fiz. **31**, 775 (1956) [Sov. Phys. JETP **4**, 648 (1957)].

[18] P. F. Bedaque and U. van Kolck, *Nucleon-Deuteron Scattering from an Effective Field Theory*, Phys. Lett. **B** **428**, 221 (1998).

[19] I. V. Brodsky, A. V. Klaptsov, M. Yu. Kagan, R. Combescot, and X. Leyronas, *Bound States of Three and Four Resonantly Interacting Particles*, JETP Lett. **82**, 273 (2005); *Exact Diagrammatic Approach for Dimer-Dimer Scattering and Bound States of Three and Four Resonantly Interacting Particles*, Phys. Rev. A **73**, 032724 (2006).

[20] M. Combescot, O. Betbeder-Matibet, and F. Dubin, *The Many-Body Physics of Composite Bosons*, Phys. Rep. **463**, 215 (2008).

[21] S. Y. Shiau, M. Combescot, and Y. C. Chang, *Trion Ground State, Excited States, and Absorption Spectrum Using Electron-Exciton Basis*, Phys. Rev. B **86**, 115210 (2012).

[22] J. Schwinger, *Coulomb Green's Function*, J. Math. Phys. (N.Y.) **5**, 1606 (1964); for the 2D version, see W. Dittrich, *The Coulomb Green's Function in Two Dimensions*, Am. J. Phys. **67**, 768 (1999).

[23] M. Combescot, R. Combescot, and F. Dubin, *Bose-Einstein Condensation and Indirect Excitons: A Review*, Rep. Prog. Phys. **80**, 066501 (2017).

[24] A very useful reference is C. R. Frye and C. Efthimiou, *Spherical Harmonics in p Dimensions*, arXiv:1205.3548.