text:
# Two-body correlations in 𝑁-body boson systems

## I Introduction

The average properties of an $`N`$-body system are often investigated in the mean-field approximation bay96 ; dal99 ; pet01 , where the wave function is a product of one-particle amplitudes. This excludes a priori effects of particle correlations, which often are responsible for decisive features or phenomena. A prominent example is the use of the bare nucleon-nucleon interaction in computations of nuclear ground state structure. The Hartree-Fock results are catastrophic, either producing unbound systems or collapsed point-like structures with infinite binding energy sie87 . The correlations must be dealt with here. One way is to use effective interactions and maintain the same Hilbert space of independent particles bra85 . One can get a long way by this procedure, but the particles are still not correlated in the wave functions, although reasonable sizes and binding energies are found bra85 ; sie87 ; dal99 . Any effect of correlations can therefore not be tested experimentally except as deviations from the mean-field results.

Another example, where correlations are essential, is the decay of Bose-Einstein condensates nie99 ; esr99 , which cannot be referred to as the lowest energy solution. Indeed, the $`N`$-body system has many lower lying states as immediately realized from the fact that two and more atoms form bound states. The energy is then already lower than the energy of the condensate, which therefore eventually decays into these lower lying structures. The decay can be either three-body recombination possibly enhanced by the presence of other particles nie99 ; esr99 ; bra02b , or macroscopic collapse into a structure of very small spatial dimension, which subsequently then recombine into more favorable bound states pit96 ; adh02b ; don01 .

The direct recombination process is increasingly probable with larger scattering length nie99 ; esr99 , and also the macroscopic collapse must increase strongly with scattering length sor02 . The limit of infinite scattering length is traditionally considered as difficult to solve as exemplified by the three-body system where the delicate Efimov states could occur efi70 ; fed93 . The structure of the $`N`$-body system in this limit is at least as difficult as the three-body problem nie01 , but the growing interest lie98 ; lie01 ; blu01 ; cow01  and the difficulties in this area demand new approaches.

The macroscopic collapse is conceptually very similar to the nuclear fission process where the (liquid) nucleus in a collective process is divided into two or more pieces. One difference is that in fission the fragments move away from each other, whereas (gaseous) Bose-Einstein condensates first collapse into a more dense state and after recombination into smaller subsystems the fragments move apart don01 . These collapse mechanisms can to a large extent be described in a mean-field picture boh98 ; adh02b , but especially for condensates the correlations could be very important first to establish the collective coordinate and the corresponding potential energy and second by influencing the collapse process itself sor01 ; sor02 .

The essential ingredients in a description of correlations in $`N`$-body systems are the techniques used to solve few-body problems car98 ; nie01 . Isolated two-body systems are easily solved and the key is very likely in handling of the three-body problem. This expectation arises since two particles plus all remaining particles effectively is a three-body problem, and this type of two-body correlations beyond the mean-field inside the $`N`$-body system is probably dominating. This is also the philosophy in the Faddeev fad61  and Yakubovsky yak67  equations, where the two-body amplitudes eventually are the basic quantities. The Yakubovsky reduction is rigorous, but very cumbersome in its full glory. However, such formulations may provide inspiration to practical approximations.

Unfortunately, it is impossible to include correlations directly on top of mean-field calculations with attractive zero range interactions phi97 ; fed01 . The Gross-Pitaevskii equation without a confining external trap can only be used for repulsion in a product wave function. The lowest state for an attraction would correspond to a divergent collapsed non-physical solution, but an additional external field is able to hold metastable solutions at larger distances. Skyrme Hartree-Fock is more sophisticated and an attraction leading to a bound system does avoid collapsed solutions bra85 . Including correlations triggers the Thomas effect tho35 , where collapsed three-body states would appear inside the many-body system fed01 . Renormalization can be invoked to cure these divergences and maintain the simplicity of the zero-range interaction fed01b ; ols01 , but correlations are still not included. Obviously a finite range interaction would also prevent the disaster when correlations are allowed.

To study correlations the form of the wave function must be flexible enough to include the corresponding degrees of freedom. As we learned from the mean-field approximation sizes and binding energies may be rather accurate even with an imprecise wave function neg82 ; sor01 . Thus for a specific purpose it seems possible to design a relatively simple wave function, which for other purposes may be a rather poor approximation, but precisely accounting for the desired degrees of freedom. A promising form of a correlated wave function suggested for nucleons rip84  was recently extended to more general systems bar99a . The simplest structure is clearly found in Bose-Einstein condensates. Here was recently introduced a formulation with generalized hyperspherical coordinates and an adiabatic expansion with the hyperradius as the adiabatic coordinate boh98 . Only the crude approximations of a zero-range interaction and the lowest (constant and thus non-correlated) hyperspherical angular wave function were used in boh98 .

In this paper we shall use the simplifications of the completely symmetric wave function. These systems reveal features only explainable by correlations and we can then already here obtain results of practical interest. The purpose of this paper is to formulate the details of a recently suggested sor01  method to include two-body correlations in an $`N`$-body system of identical bosons. We shall describe the qualitative features of the key quantities, give numerical examples, compare with established methods and available experimental data and provide predictions going beyond the present knowledge. The model is apparently rich with many unexplored possibilities perhaps especially by application to systems with large scattering length, to the dynamical evolution, and to other systems of different symmetries.

## II The $`N`$-body problem and the hyperspherical formulation

A formulation accounting for particle correlations in an $`N`$-body system must by definition go beyond the mean-field approximation. The Faddeev-Yakubovsky equations could be an appropriate starting point. However, first the form of the Hamiltonian describing the system must be decided. Then a convenient set of coordinates must be chosen in harmony with the formulation and the anticipated approximations. Derivations of suitable equations of motion are then possible and their properties can be investigated. This section describes how we choose two-body interactions, hyperspherical coordinates, a Faddeev-like decomposition of the wave function, use an adiabatic expansion, and assume at first only s-waves.

### II.1 The hyperspherical coordinates

The system of $`N`$ identical interacting bosons of mass $`m`$ may be described by $`N`$ coordinate vectors $`\stackrel{}{r}_i`$ and momenta $`\stackrel{}{p}_i`$, labeling the particles by the index $`i=1,\mathrm{},N`$. Here a more suitable choice of coordinates is the center of mass coordinates $`\stackrel{}{R}=_{i=1}^N\stackrel{}{r}_i/N`$, the $`N1`$ relative Jacobi vectors $`\stackrel{}{\eta }_k`$ with $`k=1,2,\mathrm{},N1`$ bar99a ; bar99b

$`\stackrel{}{\eta }_k`$ $`=`$ $`\sqrt{{\displaystyle \frac{Nk}{Nk+1}}}\left(\stackrel{}{r}_{Nk+1}{\displaystyle \frac{1}{Nk}}{\displaystyle \underset{j=1}{\overset{Nk}{}}}\stackrel{}{r}_j\right),`$ (1)
and their associated momenta. These Jacobi coordinates are illustrated for up to six particles in fig. 9a. We use the notation $`\eta _k|\stackrel{}{\eta }_k|`$, so $`\eta _{N1}`$ is proportional to the distance between particles $`1`$ and $`2`$, $`\eta _{N2}`$ is proportional to the distance between particle $`3`$ and the center of mass of $`1`$ and $`2`$, $`\eta _{N3}`$ is proportional to the distance between particle $`4`$ and the center of mass of the first three particles, etc.

Hyperspherical coordinates are now defined in relation to the Jacobi vectors. One length, the hyperradius $`\rho `$, is defined by

$`\rho _l^2{\displaystyle \underset{k=1}{\overset{l}{}}}\eta _k^2,\rho ^2\rho _{N1}^2={\displaystyle \frac{1}{N}}{\displaystyle \underset{i<j}{\overset{N}{}}}r_{ij}^2,`$ (2)
where $`r_{ij}|\stackrel{}{r}_i\stackrel{}{r}_j|`$. The $`N2`$ hyperangles $`\alpha _k[0,\pi /2]`$ for $`k=2,3,\mathrm{},N1`$ relate the length of the Jacobi vectors to the hyperradius via the definition

$`\mathrm{sin}\alpha _k{\displaystyle \frac{\eta _k}{\rho _k}}.`$ (3)
Since $`\rho _1=\eta _1`$ the variable $`\alpha _1=\pi /2`$ is superfluous, but is for convenience often included in the notation. Remaining are the $`2(N1)`$ angles $`\mathrm{\Omega }_\eta ^{(k)}=(\theta _k,\phi _k)`$ for $`k=1,2,\mathrm{},N1`$ defining the directions of the $`N1`$ $`\stackrel{}{\eta }_k`$-vectors, i.e. $`\theta _k[0,\pi ]`$ and $`\phi _k[0,2\pi ]`$. All angles are collectively denoted by $`\mathrm{\Omega }\{\alpha _k,\theta _k,\phi _k\},k=1,2,\mathrm{},N1`$. In total $`\mathrm{\Omega }`$ and $`\rho `$ amount to $`3(N1)`$ degrees of freedom and the center of mass coordinates $`\stackrel{}{R}`$ amount to three. These coordinates are connected by

$`{\displaystyle \underset{i=1}{\overset{N}{}}}r_i^2={\displaystyle \frac{1}{N}}{\displaystyle \underset{i<j}{\overset{N}{}}}r_{ij}^2+{\displaystyle \frac{1}{N}}\left({\displaystyle \underset{i=1}{\overset{N}{}}}\stackrel{}{r}_i\right)^2=\rho ^2+NR^2.`$ (4)
The total volume element is $`_{i=1}^Nd^3\stackrel{}{r}_i=N^{3/2}d^3\stackrel{}{R}_{k=1}^{N1}d^3\stackrel{}{\eta }_k`$, where the part depending on relative coordinates is $`_{k=1}^{N1}d^3\stackrel{}{\eta }_k`$. In hyperspherical coordinates this relative part becomes

$`{\displaystyle \underset{k=1}{\overset{N1}{}}}d^3\stackrel{}{\eta }_k=d\rho \rho ^{3N4}d\mathrm{\Omega }_{N1},`$ (7)
$`d\mathrm{\Omega }_k=d\mathrm{\Omega }_\alpha ^{(k)}d\mathrm{\Omega }_\eta ^{(k)}d\mathrm{\Omega }_{k1},`$
$`d\mathrm{\Omega }_\alpha ^{(k)}=d\alpha _k\mathrm{sin}^2\alpha _k\mathrm{cos}^{3k4}\alpha _k,`$
where $`d\mathrm{\Omega }_\eta ^{(k)}=d\theta _k\mathrm{sin}\theta _kd\phi _k`$ is the familiar angular volume element in spherical coordinates. The recursion stops at $`d\mathrm{\Omega }_1=d\mathrm{\Omega }_\eta ^{(1)}`$. Since the angle $`\alpha _{N1}`$ is related directly to a two-body distance, $`r_{12}`$, by $`\mathrm{sin}\alpha _{N1}=\eta _{N1}/\rho _{N1}=r_{12}/(\sqrt{2}\rho )`$, the volume element in eq. (7) related to this angle is especially important: $`d\mathrm{\Omega }_\alpha ^{(N1)}=d\alpha _{N1}\mathrm{sin}^2\alpha _{N1}\mathrm{cos}^{3N7}\alpha _{N1}`$.

The angular volume integrals can be computed sch68 , i.e. $`𝑑\mathrm{\Omega }_\alpha ^{(k)}=\sqrt{\pi }\mathrm{\Gamma }(3(k1)/2)/(4\mathrm{\Gamma }(3k/2))`$ and $`𝑑\mathrm{\Omega }_\eta ^{(k)}=4\pi `$, where $`\mathrm{\Gamma }(x)`$ is the gamma function.

Using eq. (7) an angular matrix element of an operator $`\widehat{O}`$ for fixed $`\rho `$ is then

$`\mathrm{\Psi }|\widehat{O}|\mathrm{\Phi }_\mathrm{\Omega }={\displaystyle 𝑑\mathrm{\Omega }_{N1}\mathrm{\Psi }^{}(\rho ,\mathrm{\Omega })\widehat{O}\mathrm{\Phi }(\rho ,\mathrm{\Omega })},`$ (8)
which in general is a function of $`\rho `$.

### II.2 Hamiltonian

We consider $`N`$ identical particles of mass $`m`$ interacting through short range two-body potentials. For completeness we will throughout the paper include an external trap confining the particles to a limited region of space through a simple three-dimensional harmonic oscillator potential $`m\omega ^2r_i^2/2`$. The trap is relevant for Bose-Einstein condensation but can easily be omitted. With a two-body central potential $`V_{ij}=V(r_{ij})`$ the total Hamiltonian is given by

$`\widehat{H}_{\text{total}}={\displaystyle \underset{i=1}{\overset{N}{}}}\left({\displaystyle \frac{\widehat{p}_i^2}{2m}}+{\displaystyle \frac{1}{2}}m\omega ^2r_i^2\right)+{\displaystyle \underset{i<j}{\overset{N}{}}}V(r_{ij}),`$ (9)
which, using eq. (4), is separable into a part only involving the center of mass coordinates and a part only involving relative coordinates. We can see this by subtracting

$`\widehat{H}_{\text{cm}}{\displaystyle \frac{\widehat{P}_R^2}{2M}}+{\displaystyle \frac{1}{2}}M\omega ^2R^2,`$ (10)
where $`\stackrel{}{P}_R_i\stackrel{}{p}_i`$ is the total momentum and $`M=Nm`$ is the total mass of the system, i.e.

$`\widehat{H}\widehat{H}_{\text{total}}\widehat{H}_{\text{cm}}={\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{\widehat{p}_i^2}{2m}}{\displaystyle \frac{\widehat{P}_R^2}{2M}}+`$
$`{\displaystyle \underset{i}{\overset{N}{}}}{\displaystyle \frac{1}{2}}m\omega ^2r_i^2{\displaystyle \frac{1}{2}}M\omega ^2R^2+{\displaystyle \underset{i<j}{\overset{N}{}}}V_{ij}.`$ (11)
Using eq. (4) and denoting the intrinsic kinetic energy operator by $`\widehat{T}_{i=1}^N\widehat{p}_i^2/(2m)\widehat{P}_R^2/(2M)`$ we can write

$`\widehat{H}=\widehat{T}+{\displaystyle \frac{1}{2}}m\omega ^2\rho ^2+{\displaystyle \underset{i<j}{\overset{N}{}}}V_{ij}.`$ (12)
In hyperspherical coordinates $`\widehat{T}`$ can be rewritten as bar99a ; bar99b

$`\widehat{T}={\displaystyle \frac{\mathrm{}^2}{2m}}\left[{\displaystyle \frac{1}{\rho ^{3N4}}}{\displaystyle \frac{}{\rho }}\rho ^{3N4}{\displaystyle \frac{}{\rho }}{\displaystyle \frac{\widehat{\mathrm{\Lambda }}_{N1}^2}{\rho ^2}}\right],`$ (13)
with the dimensionless angular kinetic energy operator $`\widehat{\mathrm{\Lambda }}_{N1}^2`$ recursively defined by

$`\widehat{\mathrm{\Lambda }}_k^2=\widehat{\mathrm{\Pi }}_k^2+{\displaystyle \frac{\widehat{\mathrm{\Lambda }}_{k1}^2}{\mathrm{cos}^2\alpha _k}}+{\displaystyle \frac{\widehat{l}_k^2}{\mathrm{sin}^2\alpha _k}},`$ (14)
$`\widehat{\mathrm{\Pi }}_k^2={\displaystyle \frac{^2}{\alpha _k^2}}+{\displaystyle \frac{3k6(3k2)\mathrm{cos}2\alpha _k}{\mathrm{sin}2\alpha _k}}{\displaystyle \frac{}{\alpha _k}},`$ (15)
where $`\mathrm{}\widehat{l}_k`$ is the angular momentum operator associated with $`\stackrel{}{\eta }_k`$. The recursion stops at $`\widehat{\mathrm{\Lambda }}_1^2=\widehat{l}_1^2`$. The angular kinetic energy operator is thus a sum of derivatives with respect to the various hyperspherical angles. Convenient transformations to get rid of first derivatives in eqs. (13) and (15) are

$`{\displaystyle \frac{2m}{\mathrm{}^2}}\widehat{T}_\rho \rho ^{(3N4)}{\displaystyle \frac{}{\rho }}\rho ^{3N4}{\displaystyle \frac{}{\rho }}=`$ (16)
$`\rho ^{(3N4)/2}\left({\displaystyle \frac{^2}{\rho ^2}}{\displaystyle \frac{(3N4)(3N6)}{4\rho ^2}}\right)\rho ^{(3N4)/2},`$
$`\widehat{\mathrm{\Pi }}_k^2=\mathrm{sin}^1\alpha _k\mathrm{cos}^{(3k4)/2}\alpha _k({\displaystyle \frac{^2}{\alpha _k^2}}{\displaystyle \frac{9k10}{2}}+`$
$`{\displaystyle \frac{(3k4)(3k6)}{4}}\mathrm{tan}^2\alpha _k)\mathrm{sin}\alpha _k\mathrm{cos}^{(3k4)/2}\alpha _k.`$ (17)
The Hamiltonian $`\widehat{H}`$ can now be rewritten as

$`\widehat{H}`$ $`=`$ $`\widehat{T}_\rho +{\displaystyle \frac{1}{2}}m\omega ^2\rho ^2+{\displaystyle \frac{\mathrm{}^2}{2m\rho ^2}}\widehat{h}_\mathrm{\Omega },`$ (18)
$`\widehat{h}_\mathrm{\Omega }`$ $``$ $`\widehat{\mathrm{\Lambda }}_{N1}^2+{\displaystyle \underset{i<j}{\overset{N}{}}}v_{ij},`$ (19)
where $`v_{ij}2m\rho ^2V_{ij}/\mathrm{}^2`$ is a dimensionless potential, $`\widehat{T}_\rho `$ is the radial kinetic energy operator, and $`\widehat{h}_\mathrm{\Omega }`$ is the dimensionless angular Hamiltonian. The intrinsic Hamiltonian, $`\widehat{H}`$, thus contains a part depending on $`\rho `$ and a part, $`\widehat{h}_\mathrm{\Omega }`$, depending on both $`\rho `$ (parametrically) and $`\mathrm{\Omega }`$.

### II.3 Equations of motion

Since the total Hamiltonian is given as $`\widehat{H}_{\text{total}}=\widehat{H}_{\text{cm}}+\widehat{H}`$, the total wave function for the $`N`$-particle system can without loss of generality be written as a product of a function, $`\mathrm{{\rm Y}}`$, depending only on $`\stackrel{}{R}`$ and a function, $`\mathrm{\Psi }`$, depending on $`\rho `$ and the $`3N4`$ angular degrees of freedom collected in $`\mathrm{\Omega }`$: $`\mathrm{{\rm Y}}(\stackrel{}{R})\mathrm{\Psi }(\rho ,\mathrm{\Omega })`$. The center of mass motion for the total mass $`M=Nm`$ is given by

$`\widehat{H}_{\text{cm}}\mathrm{{\rm Y}}(\stackrel{}{R})=E_{\text{cm}}\mathrm{{\rm Y}}(\stackrel{}{R}),`$ (20)
and as seen from eq. (10) the energy spectrum is that of a harmonic oscillator, i.e. $`E_{\text{cm},n}=\mathrm{}\omega (n+3/2)`$, where $`n`$ is a non-negative integer.

The relative wave function $`\mathrm{\Psi }(\rho ,\mathrm{\Omega })`$, obeying the stationary Schrödinger equation,

$`\widehat{H}\mathrm{\Psi }(\rho ,\mathrm{\Omega })=E\mathrm{\Psi }(\rho ,\mathrm{\Omega }),`$ (21)
is for each value of the hyperradius $`\rho `$ expanded as

$`\mathrm{\Psi }(\rho ,\mathrm{\Omega })=\rho ^{(3N4)/2}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}f_n(\rho )\mathrm{\Phi }_n(\rho ,\mathrm{\Omega }),`$ (22)
where the factor $`\rho ^{(3N4)/2}`$ is included to eliminate first derivatives in $`\rho `$, see eq. (16). Here the hyperradial wave functions, $`f_n(\rho )`$, are the expansion coefficients for fixed $`\rho `$ on the complete set of solutions $`\mathrm{\Phi }_n(\rho ,\mathrm{\Omega })`$ obtained by solving the angular eigenvalue equation:

$`(\widehat{h}_\mathrm{\Omega }\lambda _n)\mathrm{\Phi }_n(\rho ,\mathrm{\Omega })=0,`$ (23)
where $`\lambda _n`$ is the angular eigenvalue, which depends on $`\rho `$. We will usually apply the normalization $`\mathrm{\Phi }_n|\mathrm{\Phi }_m_\mathrm{\Omega }=\delta _{n,m}`$.

In complete analogy to the technique employed for $`N=3`$ jen97  we insert eq. (22) in eq. (21), use eqs. (18) and (23), and finally project the resulting equation onto the angular eigenfunctions $`\mathrm{\Phi }_n(\rho ,\mathrm{\Omega })`$. We then arrive at a set of coupled radial equations

$`\left({\displaystyle \frac{d^2}{d\rho ^2}}{\displaystyle \frac{2mE}{\mathrm{}^2}}+{\displaystyle \frac{\lambda _n(\rho )}{\rho ^2}}+{\displaystyle \frac{(3N4)(3N6)}{4\rho ^2}}+{\displaystyle \frac{\rho ^2}{b_\text{t}^4}}Q_{nn}^{(2)}(\rho )\right)f_n(\rho )={\displaystyle \underset{n^{}n}{}}\left(2Q_{nn^{}}^{(1)}(\rho ){\displaystyle \frac{d}{d\rho }}+Q_{nn^{}}^{(2)}(\rho )\right)f_n^{}(\rho ),`$ (24)
where the trap length is $`b_\text{t}\sqrt{\mathrm{}/(m\omega )}`$ and the coupling terms $`Q_{nn^{}}^{(i)}`$ ($`Q_{nn}^{(1)}=0`$) are defined as

$`Q_{nn^{}}^{(i)}(\rho ){\displaystyle \frac{\mathrm{\Phi }_n(\rho ,\mathrm{\Omega })\left|\left[\frac{}{\rho }\right]^i\right|\mathrm{\Phi }_n^{}(\rho ,\mathrm{\Omega })_\mathrm{\Omega }}{\mathrm{\Phi }_n(\rho ,\mathrm{\Omega })|\mathrm{\Phi }_n(\rho ,\mathrm{\Omega })_\mathrm{\Omega }}}.`$ (25)
The angular eigenvalues $`\lambda _n`$ enter these coupled equations as a radial potential. The total diagonal effective radial potential, $`U_n(\rho )`$, entering on the left hand side of eq. (24) is:

$`{\displaystyle \frac{2mU_n}{\mathrm{}^2}}{\displaystyle \frac{\lambda _n}{\rho ^2}}+{\displaystyle \frac{(3N4)(3N6)}{4\rho ^2}}+{\displaystyle \frac{\rho ^2}{b_\text{t}^4}}Q_{nn}^{(2)}.`$ (26)
This includes a $`\rho ^2`$-term due to the external harmonic field, a $`\rho ^2`$ centrifugal barrier-term due to the transformation of the radial kinetic energy operator, the angular potential $`\lambda _n`$, and the diagonal coupling term $`Q_{nn}^{(2)}`$.

If the non-diagonal coupling terms are neglected, i.e. the right hand side of eq. (24) vanishes, the equations simplify significantly to

$`\left({\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \frac{d^2}{d\rho ^2}}+U_n(\rho )E_{n,q}\right)f_{n,q}(\rho )=0.`$ (27)
We have here included another index $`q`$ to indicate a spectrum of energies and radial wave functions obtained for each angular solution $`n`$.

### II.4 The angular eigenvalue equation

The eigenvalue $`\lambda `$ from eq. (23) is the key quantity carrying essentially all information about the two-body interactions and therefore about possible correlations as well. The technique and approximations used to find $`\lambda `$ are then especially important.

#### II.4.1 Decomposition of the angular wave function

The angular eigenvalue is obtained by solving eq. (23). For each $`\rho `$ we first assume a decomposition of the wave function $`\mathrm{\Phi }`$ (omitting the index $`n`$) in additive components $`\mathrm{\Phi }_{ij}`$, i.e.

$`\mathrm{\Phi }(\rho ,\mathrm{\Omega })={\displaystyle \underset{i<j}{\overset{N}{}}}\mathrm{\Phi }_{ij}(\rho ,\mathrm{\Omega }),`$ (28)
where each term $`\mathrm{\Phi }_{ij}`$ is a function of $`\rho `$ and all angular coordinates $`\mathrm{\Omega }`$. This decomposition is in principle exact, since each term in itself is sufficient when all $`\mathrm{\Omega }`$ degrees of freedom are allowed. At first this ansatz seems clumsy by introducing an overcomplete basis. However, the indices $`i`$ and $`j`$ indicate special emphasis on the particle pair $`ij`$. The component $`\mathrm{\Phi }_{ij}`$ is expected to carry the information associated with two-body correlations of this particular pair. This has no significance before it is exploited in numerical techniques or subsequent approximations.

Rewriting the wave function obeying the Schrödinger equation as a sum of terms has been very successful in three-body computations. The advantage is that the correct boundary conditions are simpler to incorporate as expressed in the original formulation by Faddeev fad61  intended for scattering. Still mathematically nothing is gained or lost in this Faddeev-type of decomposition. For very weakly bound and spatially very extended three-body systems, s-waves in each of the Faddeev components are sufficient to describe the system. This is exceedingly pronounced for large scattering lengths where the delicate Efimov states appear jen97 ; nie01 .

The present $`N`$-body problem is of course in general more complicated. However, for dilute condensates essential similarities remain, i.e. the relative motion of two particles on average far from each other is most likely dominated by s-wave contributions. Each particle cannot detect any directional preference arising from higher partial waves. Only the monopole prevails. Implementing these ideas in the present context imply that each amplitude $`\mathrm{\Phi }_{ij}`$ for a fixed $`\rho `$ only should depend on the distance $`r_{ij}`$ between the two particles. For that purpose we define a two-index parameter $`\alpha _{ij}`$ by

$`\mathrm{sin}\alpha _{ij}{\displaystyle \frac{r_{ij}}{\sqrt{2}\rho }},`$ (29)
which is distinctively different from the $`\alpha _k`$’s of eq. (3). Thus we assume

$`\mathrm{\Phi }_{ij}(\rho ,\mathrm{\Omega })\varphi _{ij}(\rho ,\alpha _{ij}).`$ (30)
The boson symmetry implies that all the functions $`\varphi _{ij}`$ are equal and that we should not distinguish, so we therefore omit the indices. We then arrive at the angular wave function

$`\mathrm{\Phi }(\rho ,\mathrm{\Omega })={\displaystyle \underset{i<j}{\overset{N}{}}}\varphi (\rho ,\alpha _{ij})={\displaystyle \underset{i<j}{\overset{N}{}}}\varphi (\alpha _{ij}),`$ (31)
where we used $`\varphi _{ij}(\rho ,\alpha _{ij})=\varphi (\rho ,\alpha _{ij})\varphi (\alpha _{ij})`$ with omission of the coordinate $`\rho `$ in the last notation. The wave function in eq. (31) is symmetric with respect to interchange of two particles, $`ij`$, since $`\alpha _{ij}=\alpha _{ji}`$ and since terms like $`\varphi (\alpha _{ik})+\varphi (\alpha _{jk})`$ always appear symmetrically.

This ansatz of only s-waves dramatically simplifies the angular wave function. The original overcomplete Hilbert space is now reduced, so not every angular wave function can be expressed in this remaining basis. Thus rigorously the reduction resulted in an incomplete basis, but the degrees of freedom remaining in eq. (31) are expected to be precisely those needed to describe the main features of the condensate. The approximations are tailored to the problem under investigation.

#### II.4.2 Faddeev-like equations

Inserting the ansatz of eq. (31) along with eq. (19) in eq. (23) yields the angular equation

$`\left(\widehat{\mathrm{\Lambda }}_{N1}^2+{\displaystyle \underset{k<l}{\overset{N}{}}}v_{kl}\lambda \right){\displaystyle \underset{i<j}{\overset{N}{}}}\varphi _{ij}=0,`$ (32)
with $`\varphi _{ij}=\varphi (\alpha _{ij})`$. Rearrangement of summations leads to

$`{\displaystyle \underset{i<j}{\overset{N}{}}}\left[\left(\widehat{\mathrm{\Lambda }}_{N1}^2\lambda \right)\varphi _{ij}+v_{ij}{\displaystyle \underset{k<l}{\overset{N}{}}}\varphi _{kl}\right]=0.`$ (33)
For three particles the Faddeev equations are obtained by assuming that each term in the square brackets separately is zero jen97 . The same assumption for the $`N`$-particle system results in the $`N(N1)/2`$ Faddeev-like equations

$`\left(\widehat{\mathrm{\Lambda }}_{N1}^2\lambda \right)\varphi _{ij}+v_{ij}{\displaystyle \underset{k<l}{\overset{N}{}}}\varphi _{kl}=0,`$ (34)
which are actually identical due to symmetry. We shall not in the present paper rely on the validity of this assumption, but only use it to illustrate the procedure in the general discussion.

Choosing $`i=1`$ and $`j=2`$, with the ansatz for the wave function, eq. (31), the kinetic energy operator $`\widehat{\mathrm{\Lambda }}_{N1}^2`$ in eq. (14) reduce to $`\widehat{\mathrm{\Pi }}_{N1}^2`$, because $`\widehat{\mathrm{\Lambda }}_{N2}^2\varphi _{12}=0`$ and $`\widehat{l}_{N1}^2\varphi _{12}=0`$. Since $`\stackrel{}{\eta }_{N1}=(\stackrel{}{r}_2\stackrel{}{r}_1)/\sqrt{2}`$ and $`\rho _{N1}=\rho `$ we have $`\alpha _{N1}=\alpha _{12}`$ (compare eqs. (3) and (29)), so only derivatives with respect to $`\alpha _{12}`$ remains. Thus it is convenient to introduce the notation $`\widehat{\mathrm{\Pi }}_{12}^2\widehat{\mathrm{\Pi }}_{N1}^2`$.

In the sum over angular wave function components in eq. (34) only three different types of terms appear. Assuming $`i=1`$ and $`j=2`$ these types are classified by the set $`\{k,l\}`$ either having two, one, or zero numbers coinciding with the set $`\{1,2\}`$. Then eq. (34) is rewritten as

$`0=\left(\widehat{\mathrm{\Pi }}_{12}^2+v(\alpha _{12})\lambda \right)\varphi (\alpha _{12})+`$ (35)
$`v(\alpha _{12})\left({\displaystyle \underset{l=3}{\overset{N}{}}}\varphi (\alpha _{1l})+{\displaystyle \underset{l=3}{\overset{N}{}}}\varphi (\alpha _{2l})+{\displaystyle \underset{k3,l>k}{\overset{N}{}}}\varphi (\alpha _{kl})\right),`$
with $`v(\alpha _{ij})=2m\rho ^2V(\sqrt{2}\rho \mathrm{sin}\alpha _{ij})/\mathrm{}^2`$. Multiplying this equation from the left by $`\varphi (\alpha _{12})`$, followed by integration over all angular space except $`\alpha _{12}`$ results in an integro-differential equation in $`\alpha \alpha _{12}`$ of the form

$`\left(\widehat{\mathrm{\Pi }}_{12}^2+v(\alpha )\lambda \right)\varphi (\alpha )+v(\alpha )2(N2){\displaystyle 𝑑\tau \varphi (\alpha _{13})}`$
$`+v(\alpha ){\displaystyle \frac{1}{2}}(N2)(N3){\displaystyle 𝑑\tau \varphi (\alpha _{34})}=0.`$ (36)
Here $`d\tau d\mathrm{\Omega }_{N2}`$ is the angular volume element excluding the $`\alpha `$-dependence; the normalization is $`𝑑\tau =1`$. This projection leaves for every value of $`\alpha `$ only two different integrals due to symmetry between the first and second sum in eq. (35). Both of the remaining integrals can analytically be reduced to one dimension. The results, collected in appendix B, are denoted by

$`{\displaystyle 𝑑\tau \varphi (\alpha _{34})}=\widehat{R}_{34}^{(N2)}\varphi (\alpha ),`$ (37)
$`{\displaystyle 𝑑\tau \varphi (\alpha _{13})}=\widehat{R}_{13}^{(N2)}\varphi (\alpha ),`$ (38)
where $`\widehat{R}_{ij}^{(N2)}`$ is an operator acting on the function $`\varphi (\alpha )`$ resulting in a new function of $`\alpha `$. Mathematically $`\widehat{R}`$ resembles a rotation operator, hence the choice of notation. Eq. (II.4.2) can now be written as

$`0`$ $`=`$ $`(\widehat{\mathrm{\Pi }}_{12}^2+v(\alpha )\lambda +2(N2)v(\alpha )\widehat{R}_{13}^{(N2)}`$ (39)
$`+{\displaystyle \frac{1}{2}}(N2)(N3)v(\alpha )\widehat{R}_{34}^{(N2)})\varphi (\alpha ),`$
which is linear in the function $`\varphi `$.

#### II.4.3 Angular kinetic energy eigenfunctions

We first consider eq. (39) for non-interacting particles, i.e. $`v=0`$. Using the transformation in eq. (17) to get rid of first derivatives we find

$`({\displaystyle \frac{d^2}{d\alpha ^2}}+{\displaystyle \frac{(3N7)(3N9)}{4}}\mathrm{tan}^2\alpha `$
$`{\displaystyle \frac{9N19}{2}}\lambda )\stackrel{~}{\varphi }(\alpha )=0,`$ (40)
where $`\stackrel{~}{\varphi }(\alpha )`$ is defined as the reduced angular wave function

$`\stackrel{~}{\varphi }(\alpha )\mathrm{sin}\alpha \mathrm{cos}^{(3N7)/2}\alpha \varphi (\alpha ),`$ (41)
in perfect analogy to the transformation from radial to reduced radial wave function for the two-body problem. Since $`\varphi (\alpha )`$ for a physical state cannot diverge at $`\alpha =0`$ or $`\alpha =\pi /2`$, the boundary condition for the reduced angular wave function is $`\stackrel{~}{\varphi }(0)=\stackrel{~}{\varphi }(\pi /2)=0`$.

The (non-reduced) solutions to eq. (II.4.3) is given by the Jacobi polynomials abr

$`\varphi _K(\alpha )=P_\nu ^{(1/2,(3N8)/2)}(\mathrm{cos}2\alpha ),`$ (42)
where the hyperspherical quantum number $`K=2\nu =0,2,4,\mathrm{}`$ denotes the angular kinetic energy eigenfunction with $`\nu `$ nodes. The corresponding angular eigenvalues are $`\lambda _K=K(K+3N5)`$. The lowest eigenvalue is zero corresponding to a constant eigenfunction $`P_0^{(1/2,(3N8)/2)}=1`$.

In fig. 1a are shown examples of the reduced angular kinetic energy eigenfunctions for $`N=100`$ and the lowest three eigenvalues. The constant wave function, $`\varphi _{K=0}`$, is in the figure represented by $`\stackrel{~}{\varphi }_0(\alpha )=\mathrm{sin}\alpha \mathrm{cos}^{(3N7)/2}\alpha `$, where $`|\stackrel{~}{\varphi }_0|^2`$ then simply is the volume element in $`\alpha `$-space. The oscillations are located at relatively small $`\alpha `$-values. As seen in fig. 1b their amplitudes decrease as $`1/\sqrt{N}`$ due to the centrifugal barrier proportional to $`\mathrm{tan}^2\alpha `$ in eq. (II.4.3). Thus, as $`N`$ increases the probability becomes increasingly concentrated in a smaller and smaller region of $`\alpha `$-space around $`\alpha =0`$.

Some solutions may be spurious, i.e. each component $`\varphi `$ is non-vanishing but the full wave function $`\mathrm{\Phi }`$ in eq. (31) is identically zero. From eqs. (32) and (33) we see that the component $`\varphi `$ for a spurious solution is an eigenfunction of the angular kinetic energy operator. This special situation occurs for the $`K=2`$ eigenfunction in eq. (42), which satisfies

$`{\displaystyle 𝑑\tau \underset{i<j}{\overset{N}{}}\varphi _{K=2}(\alpha _{ij})}=0.`$ (43)
Solutions like $`\varphi _{K=2}`$ obtained by solving the full equation in eq. (33) are therefore independent of the interactions and the eigenvalue is independent of $`\rho `$.

### II.5 The $`\lambda `$-spectrum for large $`\rho `$

For large values of $`\rho `$ the short-range two-body potential $`v`$ with range $`b`$ is non-vanishing only when $`\alpha `$ is smaller than a few times $`b/\rho `$. For larger values of $`\alpha `$ the “rotation” terms, $`\widehat{R}\varphi `$, in the angular equation in eq. (39) can then be omitted.

For small values of $`\alpha `$ we find by substitution of $`r=\sqrt{2}\rho \mathrm{sin}\alpha `$ instead of $`\alpha `$ in eq. (39) that the remaining equation (without rotation terms) can be rewritten as

$`0=\left({\displaystyle \frac{\mathrm{}^2}{2\mu }}{\displaystyle \frac{d^2}{dr^2}}+V(r)E\right)u(r),`$ (44)
where $`2m\rho ^2E/\mathrm{}^2=\lambda +9N/29`$, $`u(\sqrt{2}\rho \mathrm{sin}\alpha )=\stackrel{~}{\varphi }(\alpha )`$, and the reduced mass $`\mu =m/2`$.

Each two-body bound state solution $`E<0`$ then corresponds to an eigenvalue $`\lambda `$ diverging towards $`\mathrm{}`$ as $`\rho ^2`$. Such solutions do not produce significant rotation terms, which is consistent with the omission in the derivation. The structure of the $`N`$-body system is given by the fully symmetrized wave function of two particles in the bound state, while all other particles are far away thus producing the large $`\rho `$.

All other solutions to eq. (39) without bound two-body states correspond to wave functions distributed over larger regions of $`\alpha `$-space. As the potentials then vanish for large $`\rho `$ we are left with the free solutions, i.e. the free spectrum of non-negative $`\lambda `$-values is obtained in this limit of large $`\rho `$ in addition to the discussed diverging eigenvalues for bound two-body states.

When the attractive potential contains precisely one bound two-body state with zero energy, the eigenvalue $`\lambda `$ approaches a negative constant for large $`\rho `$. This can be understood by considering a two-body state with energy slightly below zero, which forces $`\lambda `$ to diverge slowly as $`\rho ^2`$. On the other hand, if the two-body system is slightly unbound, $`\lambda `$ instead converges to zero, which is the lowest eigenvalue of the free solutions. Precisely at the threshold it seems that $`\lambda `$ cannot decide and therefore remains constant. Thus for infinite two-body s-wave scattering length one $`\lambda `$-value approaches a negative constant reached for large $`\rho `$, when the average distance between particles is much larger than the range of the interaction.

## III Angular variational equation

The angular equation is essential and solutions are not easily obtained. Proceeding with the Faddeev equations is one option, but we prefer first to derive the optimal angular equation within the Hilbert space defined by the form of the angular wave function in eq. (31). We also want to exploit that the two-body interaction is of very short range compared with the size of the system.

### III.1 Variational approach

The angular Schrödinger equation for fixed $`\rho `$ in eq. (23) and the ansatz of the wave function in eq. (31) allow us to express the eigenvalue as an expectation value, i.e.

$`\lambda ={\displaystyle \frac{\mathrm{\Phi }|\widehat{h}_\mathrm{\Omega }|\mathrm{\Phi }_\mathrm{\Omega }}{\mathrm{\Phi }|\mathrm{\Phi }_\mathrm{\Omega }}}={\displaystyle \frac{\varphi _{ij}|\widehat{h}_\mathrm{\Omega }|\mathrm{\Phi }_\mathrm{\Omega }}{\varphi _{ij}|\mathrm{\Phi }_\mathrm{\Omega }}}.`$ (45)
Since $`\widehat{h}_\mathrm{\Omega }`$ is invariant with respect to interchange of particles, the terms $`\varphi _{ij}|\widehat{h}_\mathrm{\Omega }|\mathrm{\Phi }_\mathrm{\Omega }=\varphi _{kl}|\widehat{h}_\mathrm{\Omega }|\mathrm{\Phi }_\mathrm{\Omega }`$ are identical and eq. (45) simplifies to

$`\lambda ={\displaystyle \frac{\varphi _{12}|\widehat{h}_\mathrm{\Omega }|\varphi _{kl}_\mathrm{\Omega }}{\varphi _{12}|\varphi _{kl}_\mathrm{\Omega }}}.`$ (46)
The total angular volume element is $`d\mathrm{\Omega }_{N1}=d\mathrm{\Omega }_\alpha ^{(N1)}d\mathrm{\Omega }_\eta ^{(N1)}d\mathrm{\Omega }_{N2}`$, see eq. (7). The integrands are independent of $`\mathrm{\Omega }_\eta ^{(N1)}`$ allowing to omit $`d\mathrm{\Omega }_\eta ^{(N1)}`$ from the integrations. Using eq. (8) we then obtain

$`{\displaystyle 𝑑\mathrm{\Omega }_\alpha ^{(N1)}\varphi _{12}^{}𝑑\mathrm{\Omega }_{N2}\left(\widehat{h}_\mathrm{\Omega }\lambda \right)\underset{k<l}{\overset{N}{}}\varphi _{kl}}=0.`$ (47)
The wave function component $`\varphi _{12}^{}`$ is now varied until the lowest eigenvalue is obtained. This gives the integro-differential equation

$`{\displaystyle 𝑑\mathrm{\Omega }_{N2}\underset{k<l}{\overset{N}{}}\left[(\widehat{\mathrm{\Lambda }}_{N1}^2\lambda )\varphi _{kl}+v_{kl}\underset{m<n}{\overset{N}{}}\varphi _{mn}\right]}=0,`$ (48)
where the unknown functions, $`\varphi _{ij}=\varphi (\alpha _{ij})`$, in fact all are the same identical functions of the different coordinates $`\alpha _{ij}`$. This result resembles the Faddeev-like equations of eq. (33). Many terms are identical, e.g. $`𝑑\mathrm{\Omega }_{N2}v_{12}\varphi _{34}=𝑑\mathrm{\Omega }_{N2}v_{12}\varphi _{56}`$, since particles $`1`$ and $`2`$ cannot distinguish between other pairs of particles, see appendix C.1 for the details. Collecting all terms yields

$`{\displaystyle 𝑑\mathrm{\Omega }_{N2}\left[\left(\widehat{\mathrm{\Pi }}_{12}^2+v_{12}\lambda \right)\varphi _{12}+G(\tau ,\alpha _{12})\right]}=0,`$ (49)
where $`\tau `$ denotes angular coordinates apart from $`\alpha _{12}`$, and

$`G(\tau ,\alpha _{12})={\displaystyle \frac{1}{2}}n_2\left(\widehat{\mathrm{\Pi }}_{34}^2+v(\alpha _{12})+v(\alpha _{34})\lambda \right)\varphi (\alpha _{34})`$
$`+{\displaystyle \frac{1}{2}}n_2v(\alpha _{34})\varphi (\alpha _{12})+2n_1v(\alpha _{13})\left(\varphi (\alpha _{12})+\varphi (\alpha _{23})\right)`$
$`+2n_1\left(\widehat{\mathrm{\Pi }}_{13}^2+v(\alpha _{12})+v(\alpha _{13})\lambda \right)\varphi (\alpha _{13})`$
$`+n_3\left(v(\alpha _{34})\left(\varphi (\alpha _{35})+\varphi (\alpha _{15})\right)+v(\alpha _{13})\varphi (\alpha _{45})\right)`$
$`+2n_2v(\alpha _{13})\left(\varphi (\alpha _{14})+\varphi (\alpha _{24})+\varphi (\alpha _{34})\right)`$
$`+2n_2v(\alpha _{34})\varphi (\alpha _{13})+{\displaystyle \frac{1}{4}}n_4v(\alpha _{34})\varphi (\alpha _{56}),`$ (50)
where $`n_i=_{j=1}^i(Nj1)`$, and $`\widehat{\mathrm{\Pi }}_{ij}^2`$ is defined from eq. (15) with $`k=N1`$ and with $`\alpha _k`$ replaced by $`\alpha _{ij}`$. In eq. (III.1) all terms depend at most on coordinates of the six particles $`1`$-$`6`$. The first three terms in eq. (49) do not depend on the integration variables $`\tau `$ leaving only $`G(\tau ,\alpha _{12})`$ for integration.

By appropriate choices smi77  of Jacobi systems, the relevant degrees of freedom can be expressed in terms of the five vectors $`\stackrel{}{\eta }_{N1},\mathrm{},\stackrel{}{\eta }_{N5}`$. One is the argument of the variational function and not an integration variable. The remaining twelve-dimensional integral is then evaluated with the corresponding volume element $`d\tau _{i=2}^5d\mathrm{\Omega }_\alpha ^{(Ni)}d\mathrm{\Omega }_\eta ^{(Ni)}`$ where the normalization is $`𝑑\tau =1`$. Then eq. (49) becomes

$`\left(\widehat{\mathrm{\Pi }}_{12}^2+v(\alpha _{12})\lambda \right)\varphi (\alpha _{12})+{\displaystyle 𝑑\tau G(\tau ,\alpha _{12})}=0,`$ (51)
where we used that the first terms are independent of the integration variables. Eq. (51) is a linear integro-differential equation in one variable containing up to five-dimensional integrals, see appendix C.2.

### III.2 Simple models for short-range interactions

The short-range two-body interaction with s-wave scattering length $`a_\text{s}`$ has in mean-field contexts dal99  been modelled by the three-dimensional $`\delta `$-function potential

$`V(\stackrel{}{r}_{kl})={\displaystyle \frac{4\pi \mathrm{}^2a_\text{s}}{m}}\delta (\stackrel{}{r}_{kl}).`$ (52)
With the constant angular wave function, $`\mathrm{\Phi }_{K=0}=_{i<j}^N\varphi _{K=0}(\alpha _{ij})`$, the expectation value of the angular Hamiltonian $`\widehat{h}_\mathrm{\Omega }`$ becomes

$`\lambda _{K=0}=\mathrm{\Phi }_{K=0}|\widehat{h}_\mathrm{\Omega }|\mathrm{\Phi }_{K=0}_\mathrm{\Omega }`$
$`=\mathrm{\Phi }_{K=0}\left|{\displaystyle \underset{k<l}{\overset{N}{}}}v_{kl}\right|\mathrm{\Phi }_{K=0}_\mathrm{\Omega },`$ (53)
without contribution from angular kinetic energy. With the $`\delta `$-interaction in eq. (52) we obtain (see also boh98 )

$`\lambda _{K=0}^\delta =\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{3N3}{2}\right)}{\mathrm{\Gamma }\left(\frac{3N6}{2}\right)}}N(N1){\displaystyle \frac{a_\text{s}}{\rho }}.`$ (54)
The $`\delta `$-interaction however does not contain the possibility of studying the short-range properties such as bound two-body systems and collapse of a condensate into clusters. An improvement is made by using a finite-range, but still short-range, potential. When $`\rho `$ is much larger than the potential range, $`b`$, we find an expectation value of the angular energy of the same form as $`\lambda _{K=0}^\delta `$ in eq. (54):

$`\lambda _{K=0}^{\text{finite}}\stackrel{\rho b}{}\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{3N3}{2}\right)}{\mathrm{\Gamma }\left(\frac{3N6}{2}\right)}}N(N1){\displaystyle \frac{a_\text{B}}{\rho }}.`$ (55)
The strength is now collected in the parameter $`a_\text{B}`$ instead of $`a_\text{s}`$, where

$`a_\text{B}{\displaystyle \frac{m}{4\pi \mathrm{}^2}}{\displaystyle d^3\stackrel{}{r}_{kl}V(\stackrel{}{r}_{kl})}`$ (56)
is the Born-approximation to the scattering length $`a_\text{s}`$.

However, in the opposite limit, when $`\rho b`$, the result is strongly dependent on the shape of the potential. For example, for a Gaussian potential $`V(r_{kl})=V_0\mathrm{exp}(r_{kl}^2/b^2)`$ with $`a_\text{B}=\sqrt{\pi }mb^3V_0/(4\mathrm{}^2)`$ we obtain

$`\lambda _{K=0}^{\text{finite}}\stackrel{\rho b}{}{\displaystyle \frac{4}{\sqrt{\pi }}}N(N1){\displaystyle \frac{a_\text{B}}{b}}\left({\displaystyle \frac{\rho }{b}}\right)^2.`$ (57)
As seen from these two limits there are some interesting scaling-properties for finite-range potentials. The angular eigenvalue at a given $`N`$-value depends only on $`a_\text{B}/b`$ and $`\rho /b`$. More specifically we find for a Gaussian potential

$`v_{kl}={\displaystyle \frac{2m\rho ^2V_0}{\mathrm{}^2}}e^{r_{kl}^2/b^2}={\displaystyle \frac{8a_\text{B}}{\sqrt{\pi }b}}\left({\displaystyle \frac{\rho }{b}}\right)^2e^{2(\rho /b)^2\mathrm{sin}^2\alpha _{kl}},`$ (58)
which implies that for a given $`a_\text{B}/b`$ the angular eigenvalue $`\lambda `$ is a function of $`\rho /b`$ only. Moreover the radial potential can be written as

$`{\displaystyle \frac{2mb^2U}{\mathrm{}^2}}={\displaystyle \frac{\lambda }{(\rho /b)^2}}+{\displaystyle \frac{(3N4)(3N6)}{4(\rho /b)^2}}+{\displaystyle \frac{(\rho /b)^2}{(b_\text{t}/b)^4}},`$ (59)
and the scaled energy, $`2mb^2E/\mathrm{}^2`$, is then for a given $`N`$-value a function of $`a_\text{B}/b`$ and $`b_\text{t}/b`$ only. These scaling properties are useful in model calculations.

The scattering lengths are apparently essential parameters. We show in fig. 2a the s-wave scattering length, $`a_\text{s}`$, as a function of the strength parameter $`|a_\text{B}|/b`$ for both the Gaussian potential and for an attractive box potential. When the numerical value of the parameter $`a_\text{B}`$ is increased from zero, the scattering length varies slowly and roughly linearly with $`a_\text{B}`$ for small $`a_\text{B}`$ until the value $`a_\text{B}^{(0)}`$, where $`a_\text{s}`$ diverges as a signal of the appearance of the first two-body bound state. The threshold value of $`a_\text{B}^{(0)}`$ is different for the Gaussian and for the box potentials, but $`a_\text{s}/a_\text{B}`$ as a function of $`a_\text{B}/a_\text{B}^{(0)}`$ results in virtually the same curves, see fig. 2b. This indicates that for simple potentials, the behaviour is approximately shape independent.

### III.3 Short-range approximation

The angular eigenvalue equation, eq. (51), simplifies in the limit when the two-body interaction range $`b`$ is much smaller than $`\rho `$. Then the integrals are either analytical or reduce to one-dimensional integrals. This substitution of a $`\delta `$-function is only allowed for the potentials appearing under the integrals. Otherwise the Thomas collapse occurs fed01 . Thus apart from the local terms containing $`v(\alpha )`$, the results only depend on $`a_\text{B}`$ defined in eq. (56). For example the $`𝑑\tau v(\alpha _{34})`$-term reduces to

$`{\displaystyle 𝑑\tau v(\alpha _{34})}v_1(\alpha )2\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{3N6}{2}\right)}{\mathrm{\Gamma }\left(\frac{3N9}{2}\right)}}{\displaystyle \frac{a_\text{B}}{\rho \mathrm{cos}^3\alpha }}.`$ (60)
Similarly the $`𝑑\tau v(\alpha _{13})`$-term reduces to

$`{\displaystyle 𝑑\tau v(\alpha _{13})}v_2(\alpha )`$
$`{\displaystyle \frac{8}{3\sqrt{3}}}\mathrm{cos}^{3N11}\beta _0\mathrm{\Theta }(\alpha <\pi /3)v_1(\alpha ),`$ (61)
where $`\mathrm{sin}\beta _0\mathrm{tan}\alpha /\sqrt{3}`$ and $`\mathrm{\Theta }`$ is the truth function. The remaining terms can in this limit be expressed through $`v_1(\alpha )`$, $`v_2(\alpha )`$, $`\widehat{R}_{ij}^{(k)}`$ from eqs. (37) and (38), and other related rotation operators $`\widehat{R}_{ijkl}^{(n)}`$. Corresponding definitions are given in appendix C.3.

The reductions can be understood qualitatively via the diagrams in fig. 3, which shows the geometry when the zero-range interaction is contributing to the integrals. As an example, in the integral $`𝑑\tau v(\alpha _{13})\varphi (\alpha _{34})`$, see fig. 3a, particles $`1`$ and $`3`$ must be close together as shown fig. 3b. Then the distance between particles $`3`$ and $`4`$ appearing in $`\varphi _{34}`$ is approximately equal to the distance between particles $`1`$ and $`4`$. Therefore $`𝑑\tau v(\alpha _{13})\varphi (\alpha _{34})𝑑\tau v(\alpha _{13})\varphi (\alpha _{14})`$.

For $`b\rho `$ we get almost independent of the potential

$`{\displaystyle 𝑑\tau G(\tau ,\alpha )}\left({\displaystyle \frac{n_2}{2}}v_1(\alpha )+4n_1v_2(\alpha )\right)\varphi (\alpha )`$
$`+{\displaystyle \frac{n_2}{2}}\widehat{R}_{34}^{(N2)}v\varphi (\alpha )+2n_1\widehat{R}_{13}^{(N2)}v\varphi (\alpha )`$
$`+{\displaystyle \frac{n_2}{2}}\left(\widehat{R}_{34}^{(N2)}\widehat{\mathrm{\Pi }}_{34}^2\varphi (\alpha )+(v(\alpha )\lambda )\widehat{R}_{34}^{(N2)}\varphi (\alpha )\right)`$
$`+2n_1\left(\widehat{R}_{13}^{(N2)}\widehat{\mathrm{\Pi }}_{13}^2\varphi (\alpha )+(v(\alpha )\lambda )\widehat{R}_{13}^{(N2)}\varphi (\alpha )\right)`$
$`+{\displaystyle \frac{n_4}{4}}v_1(\alpha )\widehat{R}_{34}^{(N3)}\varphi (\alpha )+n_3v_1(\alpha )\widehat{R}_{3435}^{(1)}\varphi (\alpha )`$
$`+n_3v_2(\alpha )\widehat{R}_{1345}^{(1)}\varphi (\alpha )`$
$`+n_3v_1(\alpha )\widehat{R}_{13}^{(N3)}\varphi (\alpha )+2n_2v_1(\alpha )\widehat{R}_{3413}^{(2)}\varphi (\alpha )`$
$`+2n_2v_2(\alpha )\left[2\widehat{R}_{1314}^{(2)}\varphi (\alpha )+\widehat{R}_{1324}^{(2)}\varphi (\alpha )\right].`$ (62)
In the extreme zero-range limit the two terms $`\widehat{R}_{ij}^{(N2)}v\varphi `$ are proportional to $`\varphi (0)`$. Using eqs. (60) and (61) we define $`\stackrel{~}{v}`$ and $`g`$ by

$$\stackrel{~}{v}(\alpha )\frac{n_2}{2}v_1(\alpha ),g(\alpha )\stackrel{~}{v}(\alpha )2n_1v_2(\alpha )$$
(63)
and rewrite eq. (III.3) as

$`{\displaystyle 𝑑\tau G(\tau ,\alpha )}`$ $``$ $`\stackrel{~}{v}(\alpha )\left(1+2g(\alpha )\right)\varphi (\alpha )+\widehat{R}\varphi (\alpha )`$ (64)
$`+\stackrel{~}{v}(\alpha )\left(1+g(\alpha )\right)\varphi (0),`$
where $`\widehat{R}\varphi (\alpha )`$ collectively denotes the remaining rotational terms. Collecting the terms from eq. (51) then leads to an integro-differential equation in one variable:

$`0`$ $`=`$ $`\left(\widehat{\mathrm{\Pi }}_{12}^2+v(\alpha )+\stackrel{~}{v}(\alpha )\left(1+2g(\alpha )\right)\lambda \right)\varphi (\alpha )`$ (65)
$`+\widehat{R}\varphi (\alpha )+\stackrel{~}{v}(\alpha )\left(1+g(\alpha )\right)\varphi (0).`$
Eq. (65) is a linear eigenvalue equation (in $`\varphi `$) for one variable $`\alpha `$. The eigenvalue $`\lambda (\rho )`$ is the key quantity in the much simpler radial equation.

If we neglect the rotational terms in $`\widehat{R}`$, assume $`\varphi (0)0`$, and use that $`g(\alpha )1`$ when $`N`$ is large, we get a solution $`\lambda \stackrel{~}{v}(0)`$. Since $`\stackrel{~}{v}(0)=\lambda _{K=0}^{\text{finite}}`$ at large $`N`$, this eigenvalue is the same as found in eq. (55) from the expectation value of a finite range potential in a constant angular wave function. Eq. (65) thus predicts an eigenvalue proportional to $`\rho ^1`$ with a proportionality factor given by eq. (55), where the determining parameter is $`a_\text{B}`$. However, as discussed in boh98  the angular eigenvalue for small scattering length $`a_\text{s}`$ is determined by $`a_\text{s}`$ only. This is for a sufficiently weak interaction the correct limit for $`a_\text{B}`$ as discussed in connection with fig. 2.

The repulsive effective potential in $`\stackrel{~}{v}(\alpha )`$ pushes the wave function into a narrow region at small $`\alpha `$ outside the repulsive core of $`v(\alpha )`$. Such a solution is often not a good approximation because the rotational terms are important. Also for attractive potentials the confinement to small $`\alpha `$ can not be achieved by the $`\stackrel{~}{v}(\alpha )`$ term. It is therefore crucial to include all the terms of eq. (65) as discussed in sor01 .

## IV Numerical illustrations

The method has to be tested by solving the derived equations for realistic parameter choices. We first discuss qualitatively which numerical technique we apply. As we shall explain the difficulties increase as $`N`$ and $`\rho `$ increase. We shall therefore concentrate on relatively small values of $`N`$, which happens to be a region of growing interest in the art of making Bose-Einstein condensates. The first physical results focus on the angular eigenvalues and the related wave functions. Then the fully defined radial equation and its solutions are discussed.

### IV.1 Method of solving

A usual method within the hyperspherical formalism is to expand the angular wave function on kinetic energy eigenfunctions, the so-called hyperspherical harmonics lin95 ; boh98 . However, since the hyperspherical harmonics contain oscillations at about $`\alpha 1/K`$, large $`K`$’s of the order of $`K_{\text{max}}\rho /b`$ are needed to describe potentials limited to $`\alpha <b/\rho `$. This $`K_{\text{max}}`$ becomes very large for the application to Bose-Einstein condensation, so a basis of hyperspherical harmonics is not suitable in this context. Thus using only $`K=0`$ for a zero-range interaction boh98  must be far from the optimal solution.

Instead we choose a basis of discrete mesh points distributed in $`\alpha `$-space $`\varphi (\alpha )\underset{¯}{\varphi }[\varphi (\alpha _1),\mathrm{},\varphi (\alpha _m),\mathrm{},\varphi (\alpha _M)]`$ to take into account the short range of the potential and to keep sufficient information about small $`\alpha `$. Derivatives are then written as finite differences koo90  and integrations like $`\widehat{R}\varphi (\alpha )`$ of eq. (65) can be expressed in matrix form, i.e. $`\widehat{R}\varphi (\alpha )\underset{¯}{\underset{¯}{R}}\underset{¯}{\varphi }`$.

Numerical computation of the integrals becomes increasingly difficult with decreasing interaction range. This can be understood in terms of the $`\alpha _{12}`$ coordinate, since the potential at a given $`\rho `$ and a given range, $`b`$, of the interaction, is confined to an $`\alpha _{12}`$-region of size $`\mathrm{\Delta }\alpha _{12}b/\rho `$, which for Bose-Einstein condensates easily becomes very small and thus cannot be handled directly numerically. In the following we shall use the equations obtained in the short-range approximation, where the difficulties are much smaller and the physics content is still maintained.

### IV.2 Angular solutions

The behaviour of the lowest eigenvalues $`\lambda `$ depend strongly on the potential. The characteristic feature is the large distance asymptotic behaviour, i.e. divergence as $`\rho ^2`$ corresponding to a bound two-body state (see eq. (44)) or convergence as $`\rho ^1`$ towards a finite value corresponding to the spectrum for free particles. The constant of proportionality to $`\rho ^1`$ is qualitatively recovered as the predicted boh98  dependence on $`a_\text{s}`$. The exception arises for infinite scattering length at the threshold for two-body binding, where one angular eigenvalue at large $`\rho `$ approaches a negative constant. This structure is illustrated in fig. 4, where the lowest angular eigenvalue is shown for various Gaussian potential strengths covering the region from unbound to bound two-body states.

For repulsive potentials the eigenvalues are positive and the lowest approaches zero from above. For weak attractions the lowest $`\lambda `$ is negative approaching zero from below. When the attraction can bind a pair of particles the divergence is seen while the second $`\lambda `$ for the same potential is positive and approaching zero from above in qualitative agreement with the lowest eigenvalue for a repulsive potential. The higher eigenvalues would then converge to $`K(K+3N5)`$ as $`1/\rho `$, where $`K=4,6,8\mathrm{}`$. The solution for $`K=2`$ is not allowed, corresponding to removal of the non-physical spurious solution of the Faddeev equations. Increasing the attraction to allow two bound two-body states would then shift the asymptotic spectrum such that one more eigenvalue diverge while the non-negative energy spectrum remains unchanged.

At each threshold for the appearance of a new bound two-body state one eigenvalue asymptotically approaches a negative constant. This especially interesting eigenvalue is responsible for the structure of the $`N`$-body system for very large scattering lengths.

The total angular wave functions are determined as a sum of the components. We show in fig. 5 an example of the reduced wave function for a potential with one bound two-body state. The amplitude increases with $`\rho `$ and concentrates at very small values of $`\alpha `$. This reflects the convergence towards the two-body bound state in agreement with the transformation $`r_{12}=\sqrt{2}\rho \mathrm{sin}\alpha `$. Recovering this behaviour numerically is essential, otherwise the large distance behaviour cannot be described.

The angular eigenfunction varies with the strength of the interaction. Examples of this variation are shown in fig. 6. The non-interacting wave function has only nodes at the endpoints. The repulsive case shows an oscillation, which lowers the angular energy due to the rotation terms sor01 . The fast change at small $`\alpha `$ is typical for interacting particles. The wave function for the excited state has an additional node. The corresponding lower-lying wave function is shown in fig. 5.

The wave function for infinite scattering length corresponds to an interaction where the two-body bound state is at the threshold for occurrence. This eigenfunction resembles those where a bound two-body state is present, compare with the results shown in fig. 5 with different scales on both axes. However, now (thick curve of fig. 6) the wave function is located at larger $`\alpha `$-values and a node is present in the tail at an intermediate $`\alpha 0.25`$.

The structure of the component of the angular wave function is further illustrated by the second moment defined by $`r_{12}^2_\varphi 2\rho ^2\varphi |\mathrm{sin}^2\alpha |\varphi `$. A number of these moments for different interactions are shown in fig. 7 as functions of $`\rho `$. For states obtained from repulsive potentials, moderately attractive potentials without bound two-body states, and for excited states of positive $`\lambda `$, this moment, $`r_{12}^2_\varphi `$, increases proportional to $`\rho ^2`$ for large $`\rho `$. This resembles the behaviour of the expectation value in the lowest angular state for a non-interacting system, $`K=0`$, where $`r_{12}^2_\varphi =2\rho ^2/(N1)`$. The qualitative explanation is that large $`\rho `$ implies the limit of a non-interacting spectrum with the corresponding wave functions. The particles exploit the possibility of being far from each other, since there is no energetic advantage of being close due to the lack of a bound two-body state.

In contrast a different behaviour is observed when the potential can bind two particles, i.e. $`r_{12}^2_\varphi `$ approaches a constant at large $`\rho `$. This can be understood as the angular equation in this limit approaches the two-body equation in eq. (44). The radial wave function in the zero-range limit converges to $`u(r)=e^{r/a_\text{s}}`$, where $`a_\text{s}`$ is the scattering length. The second moment is then found as $`u|r^2|u=a_\text{s}^2/2`$, which in the limit of large $`\rho `$ reproduce $`r_{12}^2_\varphi `$ for $`a_\text{s}/b=9.32`$ and $`a_\text{s}/b=5.31`$, see fig. 7.

This can be expressed in a different way: when a two-body bound state is present, the angular wave function is at increasing $`\rho `$ squeezed inside the potential, since the range in $`\alpha `$-space decreases proportional to $`\rho ^1`$, and we obtain $`\varphi |\mathrm{sin}^2\alpha |\varphi 1/\rho ^2`$. The distance between a pair of particles is therefore independent of $`\rho `$ at large values of $`\rho `$. This means that pairwise the two-body bound state is approached while all other particles are far away. The symmetrization does not affect this conclusion.

At the threshold for two-body binding, infinite scattering length, we again observe the intermediate behaviour resembling a logarithmic dependence in fig. 7.

### IV.3 Radial solutions

The most interesting angular eigenvalues either converge to zero or remain constant for large $`\rho `$. We therefore select an interaction where $`\lambda `$ approaches zero relatively slowly from below. This is then a weakly attractive potential without bound two-body states although not very far from this threshold. The parameter $`a_\text{s}/b=15.7`$ corresponds to this case. The resulting $`\lambda `$ shown in fig. 4 is used to compute the radial potential in eq. (26). The total potential is shown in fig. 8.

The radial potential always diverges at large $`\rho `$ due to the harmonic external field and at small $`\rho `$ due to the centrifugal barrier term. Thus there is always infinitely many bound states. For the moderate attraction used in fig. 8 the potential for the lowest angular eigenvalue has a global negative minimum at small $`\rho `$ separated by a barrier at intermediate $`\rho `$ from a local positive minimum at $`\rho \rho _\text{t}\sqrt{3N/2}b_\text{t}`$. The minimum at small $`\rho `$ disappears quickly with a slightly less attractive potential. It becomes deeper and wider for a more attractive potential, where the barrier at intermediate $`\rho `$ simultaneously is reduced and eventually disappears completely.

The radial potential corresponding to the second adiabatic angular potential is also shown in fig. 8. Since this $`\lambda `$ is positive for all $`\rho `$ an attractive pocket ($`U<0`$) in the radial potential cannot arise. It coincides with the lowest radial potential for large $`\rho `$ and due to the lack of attraction at small $`\rho `$ the potential therefore diverges to $`+\mathrm{}`$ for small $`\rho `$ without going through another minimum.

The coupling terms of eq. (25) contribute at most about 1 % compared with other terms of the full radial equation, eq. (24), and can be omitted. We therefore arrive at solving the uncoupled radial equation, eq. (27), the solutions of which are considered in the following.

The external field is negligible when $`\rho \sqrt{N}b_\text{t}`$ and the radial potential is therefore negative when $`\lambda +(3N4)(3N6)/4<0`$ and $`\rho `$ is sufficiently small. Then genuinely bound many-body states of negative energy are possible in our model without the influence of confinement from the trap. The radial equation corresponding to the parameters of the lowest potential shown in fig. 8 has only one negative-energy solution with the wave function located in the global negative minimum.

The first of the infinitely many excited states in this potential is located in the local minimum at larger $`\rho `$ created by competition between the centrifugal barrier and the external harmonic oscillator potential. This excited state is usually referred to as the Bose-Einstein condensate, because the corresponding wave function is similar to that obtained in experiments dal99 . It also resembles the wave function found from the purely repulsive radial potentials arising from both the second $`\lambda `$ with the present attractive interaction and from repulsive two-body interactions.

From eq. (4) we have $`\rho ^2=Nr_i^2NR^2=Nr_i^23b_\text{t}^2/2`$. Using eq. (2) we get $`2\rho ^2=(N1)r_{12}^2`$. The root mean square radius $`r_i^2^{1/2}`$ of the condensate state, $`f_{\text{ex}}`$, is then 2 % larger and the energy $`E_{\text{ex}}`$ is 3 % lower than the Gross-Pitaevskii results. The lowest radial state in $`U_1`$ has both energy and root mean square radius about 10 % larger than the Gross-Pitaevskii results.

The moderate attraction used to obtain the potential in fig. 8 is therefore seen to produce a lower-lying many-body bound state with an average distance between the particles about 37 times smaller than in the condensate state. The structure of this state could as well be characterized as a condensate (condensed $`N`$-body state), but it is much more unstable due to the many orders of magnitude larger density and the subsequent larger recombination probability nie99 . This ground state, $`f_{\text{gs}}`$, has no parallel in usual mean-field computations.

Even the state with $`E_{\text{gs}}<0`$, $`\rho ^2^{1/2}216b`$, has a root mean square distance between two particles, $`r_{12}^2^{1/2}70b`$, much larger than the interaction range $`b`$. This is a sufficiently large average distance between particles to assume that the short range details of the two-body potentials are unimportant. Increasing the attraction more negative-energy states appear in the attractive pocket inside the external trap. They occur at increasingly larger densities, but still low enough for the short-range details to be unimportant, see more details in sor02 .

At the two-body threshold, when the s-wave scattering length diverges, potentially infinitely many negative-energy many-body states could occur. However, the external harmonic oscillator field limits the number $`N_E`$ of these new negative-energy states, located inside the trap and outside the two-body potential, to approximately

$`N_E0.5N\mathrm{ln}\left({\displaystyle \frac{b_\text{t}}{37b}}\right),\text{for }b_\text{t}N|a_\text{s}|.`$ (66)
The number of these multi-particle Efimov states therefore scales linearly with $`N`$ and logarithmically with the ratio of trap length and interaction range, see also sor02 .

The stability of the $`N`$-body system as a coherent object (condensate) depends on the radial potential. Starting with a condensate in the minimum of $`f_{\text{ex}}`$ in fig. 8, decay occurs both through two- and three-body recombination processes and through macroscopical tunneling through the barrier to state(s) in the global minimum. These denser states recombine faster due to the smaller interparticle distances. Increasing the attraction, the rate for macroscopical tunneling increases due to a smaller barrier. Eventually, at large (infinite) scattering length, the barrier has vanished and the $`N`$-body system can contract freely, i.e. populate the many-body Efimov states. The time scale for this contraction is estimated to be about 0.1 ms sor02 , which can be compared with experiments don01  and mean-field computations adh02b . Thus sudden removal of the barrier over the first short period of time lead to spreading of the probability to the smaller distances where recombination then takes place and the condensate decays with a corresponding lifetime. This is qualitatively in complete agreement with the measured decay function don01 .

## V Summary and conclusions

We have formulated in details a new method to investigate two-body correlations within symmetric $`N`$-body systems. In contrast to the product assumption for the wave function in the mean-field approximation we use a Faddeev-type decomposition. The overcomplete basis is simplified by using in each of the two-body components only the lowest possible number of partial waves, i.e. s-waves. The allowed Hilbert space is then dramatically reduced, it is not complete, but hopefully precisely designed to treat the two-body correlations.

We use an adiabatic hyperspherical expansion with the hyperradius as the adiabatic coordinate. We derive the optimal angular equation for the two-body amplitude (the Faddeev component) arriving at an one-dimensional angular integro-differential equation. Its eigenvalues are closely related to the effective hyperradial potential, which receives contributions from the two-body interactions, the kinetic energy operator and from the external field confining the system. This potential diverges both at small and large hyperradii due to the centrifugal barrier and the external trap, respectively.

For repulsive two-body interactions only one minimum in the radial potential exists, but for moderately attractive two-body potentials two minima are present separated by a barrier. The minimum at shortest hyperradii is at negative potential and is therefore able to bind the system without help from the confining external field. For a state in this minimum the particles are on average still far outside the range of interaction with each other. The other minimum at larger hyperradii can also be pronounced enough to hold localized stationary states, which are similar to the solutions usually referred to as Bose-Einstein condensates.

As the attraction increases and approaches the threshold for two-body binding, i.e. the scattering length increases towards infinity, the intermediate barrier disappears and only the negative minimum at smaller hyperradii remains. However, it becomes both deeper and wider and as a consequence more bound states of negative energy appear. These many-body states have characteristic features of Efimov states. They have inter-particle distances much larger than the range of the two-body interaction and are therefore universal structures independent of details of the potential. Their number increases proportional to the number of particles and logarithmically with the ratio of trap length to interaction range.

In conclusion, we have formulated a method to treat two-body correlations in an $`N`$-body system. We applied the method to Bose-Einstein condensates and obtained a simple one-dimensional pictorial description in terms of an effective length coordinate. Macroscopic collapse is conjectured to proceed via the new universal Efimov states at intermediate hyperradii, which quickly, due to the much larger density, subsequently recombine into dimer or trimer states. This decay can be studied quantitatively in the present model. Another unique feature of the model is to provide a solution even in the case of very large scattering lengths, where the Gross-Pitaevskii equation breaks.

## Appendix A Coordinates

For use in the calculation of matrix elements different Jacobi trees can be chosen smi77 . The relevant ones in this context are shown in fig. 9.

The coordinates of the standard tree of fig. 9a are defined by

$`\stackrel{}{\eta }_{N1}`$ $`={\displaystyle \frac{1}{\sqrt{2}}}(\stackrel{}{r}_2\stackrel{}{r}_1),`$ (67a)
$`\stackrel{}{\eta }_{N2}`$ $`=\sqrt{{\displaystyle \frac{2}{3}}}\left(\stackrel{}{r}_3{\displaystyle \frac{1}{2}}(\stackrel{}{r}_2+\stackrel{}{r}_1)\right),`$ (67b)
$`\stackrel{}{\eta }_{N3}`$ $`=\sqrt{{\displaystyle \frac{3}{4}}}\left(\stackrel{}{r}_4{\displaystyle \frac{1}{3}}(\stackrel{}{r}_3+\stackrel{}{r}_2+\stackrel{}{r}_1)\right),`$
$`\mathrm{}`$
$`\stackrel{}{\eta }_1=`$ $`\sqrt{{\displaystyle \frac{N1}{N}}}\left(\stackrel{}{r}_N{\displaystyle \frac{1}{N1}}(\stackrel{}{r}_{N1}+\mathrm{}+\stackrel{}{r}_1)\right).`$ (67d)
In the (12)(34)-tree of fig. 9b two of the vectors are different from the standard tree:

$`\stackrel{}{\eta }_{N2}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}(\stackrel{}{r}_4\stackrel{}{r}_3),`$ (68a)
$`\stackrel{}{\eta }_{N3}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\stackrel{}{r}_4+\stackrel{}{r}_3\stackrel{}{r}_2\stackrel{}{r}_1).`$ (68b)
In the (123)(45)-tree of fig. 9c two of the vectors differ from the standard tree:

$`\stackrel{}{\eta }_{N3}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}(\stackrel{}{r}_5\stackrel{}{r}_4),`$ (69a)
$`\stackrel{}{\eta }_{N4}`$ $`=`$ $`\sqrt{{\displaystyle \frac{6}{5}}}\left({\displaystyle \frac{1}{2}}(\stackrel{}{r}_5+\stackrel{}{r}_4){\displaystyle \frac{1}{3}}(\stackrel{}{r}_3+\stackrel{}{r}_2+\stackrel{}{r}_1)\right).`$ (69b)
In the (12)(345)-tree of fig. 9d three of the vectors deviate from the standard tree:

$`\stackrel{}{\eta }_{N2}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}(\stackrel{}{r}_4\stackrel{}{r}_3),`$ (70a)
$`\stackrel{}{\eta }_{N3}`$ $`=`$ $`\sqrt{{\displaystyle \frac{2}{3}}}(\stackrel{}{r}_5{\displaystyle \frac{1}{2}}(\stackrel{}{r}_4+\stackrel{}{r}_3)),`$ (70b)
$`\stackrel{}{\eta }_{N4}`$ $`=`$ $`\sqrt{{\displaystyle \frac{6}{5}}}\left({\displaystyle \frac{1}{3}}(\stackrel{}{r}_5+\stackrel{}{r}_4+\stackrel{}{r}_3){\displaystyle \frac{1}{2}}(\stackrel{}{r}_2+\stackrel{}{r}_1)\right).`$ (70c)
In the (12)(34)(56)-tree of fig. 9e four of the vectors are different from the standard tree:

$`\stackrel{}{\eta }_{N2}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}(\stackrel{}{r}_4\stackrel{}{r}_3),\stackrel{}{\eta }_{N3}={\displaystyle \frac{1}{\sqrt{2}}}(\stackrel{}{r}_6\stackrel{}{r}_5),`$ (71a)
$`\stackrel{}{\eta }_{N4}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\stackrel{}{r}_4+\stackrel{}{r}_3\stackrel{}{r}_2\stackrel{}{r}_2),`$ (71b)
$`\stackrel{}{\eta }_{N5}`$ $`=`$ $`\sqrt{{\displaystyle \frac{4}{3}}}\left({\displaystyle \frac{\stackrel{}{r}_6+\stackrel{}{r}_5}{2}}{\displaystyle \frac{\stackrel{}{r}_4+\stackrel{}{r}_3+\stackrel{}{r}_2+\stackrel{}{r}_1}{4}}\right).`$ (71c)
Since only inter-relations between $`\stackrel{}{\eta }_{N1}`$, $`\stackrel{}{\eta }_{N2}`$, and $`\stackrel{}{\eta }_{N3}`$ are needed in evaluating the matrix elements, we will use the common notation:

$`\eta _{N1}`$ $`=`$ $`\rho \mathrm{sin}\alpha ,\eta _{N2}=\rho \mathrm{cos}\alpha \mathrm{sin}\beta `$ (72)
$`\eta _{N3}`$ $`=`$ $`\rho \mathrm{cos}\alpha \mathrm{cos}\beta \mathrm{sin}\gamma ,\stackrel{}{\eta }_k\stackrel{}{\eta }_l=\eta _k\eta _l\mathrm{cos}\theta _{k,l},`$ (73)
where $`\theta _{k,l}`$ is the angle between the $`k`$’th and $`l`$’th Jacobi vectors. We abbreviate $`\theta _{N1,N2}\theta _x`$, $`\theta _{N1,N3}\theta _y`$, and $`\theta _{N2,N3}\theta _z`$. The azimuthal angle $`\phi `$ determining the projection of $`\stackrel{}{\eta }_{N3}`$ onto the plane of $`\stackrel{}{\eta }_{N1}`$ and $`\stackrel{}{\eta }_{N2}`$ is defined such that

$`\mathrm{cos}\theta _z`$ $`=`$ $`\mathrm{sin}\theta _x\mathrm{sin}\theta _y\mathrm{cos}\phi +\mathrm{cos}\theta _x\mathrm{cos}\theta _y.`$ (74)
With $`\tau =\{\beta ,\gamma ,\theta _x,\theta _y,\phi \}`$ a matrix element of an arbitrary function $`f`$ of all the variables $`\alpha `$ and $`\tau `$ then becomes

$`{\displaystyle 𝑑\tau f(\alpha ,\tau )}`$ $`=`$ $`{\displaystyle \frac{_0^{\pi /2}𝑑\beta \mathrm{sin}^2\beta \mathrm{cos}^{3N10}\beta _0^{\pi /2}𝑑\gamma \mathrm{sin}^2\gamma \mathrm{cos}^{3N13}\gamma _0^\pi 𝑑\theta _x\mathrm{sin}\theta _x_0^\pi 𝑑\theta _y\mathrm{sin}\theta _y_0^{2\pi }𝑑\phi f(\alpha ,\tau )}{_0^{\pi /2}𝑑\beta \mathrm{sin}^2\beta \mathrm{cos}^{3N10}\beta _0^{\pi /2}𝑑\gamma \mathrm{sin}^2\gamma \mathrm{cos}^{3N13}\gamma _0^\pi 𝑑\theta _x\mathrm{sin}\theta _x_0^\pi 𝑑\theta _y\mathrm{sin}\theta _y_0^{2\pi }𝑑\phi }},`$ (75)
where the normalization is $`𝑑\tau =1`$. We need relations for interparticle distances and define $`\stackrel{}{\eta }_{ij}(\stackrel{}{r}_j\stackrel{}{r}_i)/\sqrt{2}`$ and the angle $`\alpha _{ij}`$ related to $`\eta _{ij}=\rho \mathrm{sin}\alpha _{ij}=r_{ij}/\sqrt{2}`$.

## Appendix B Integrals (Faddeev)

Eqs. (37) and (38) are evaluated as follows.

In the integral $`𝑑\tau \varphi (\alpha _{34})`$ a convenient choice of coordinates is the alternative Jacobi (12)(34)-tree of fig. 9b given by eqs. (68). The angle $`\alpha _{34}`$ is associated with the distance $`r_{34}=\sqrt{2}\eta _{34}`$ by the relation

$`\eta _{34}=\eta _{N2}=\rho \mathrm{cos}\alpha \mathrm{sin}\beta =\rho \mathrm{sin}\alpha _{34}`$ (76)
$`\mathrm{sin}\alpha _{34}=\mathrm{cos}\alpha \mathrm{sin}\beta .`$ (77)
The integrand, $`\varphi (\alpha _{34})`$, only depends on $`\alpha _{34}`$, which is a function of $`\alpha `$ and $`\beta `$. Therefore at fixed $`\alpha `$ eq. (75) reduces to

$`{\displaystyle 𝑑\tau \varphi (\alpha _{34})}={\displaystyle \frac{_0^{\pi /2}𝑑\beta \mathrm{sin}^2\beta \mathrm{cos}^{3N10}\beta \varphi (\alpha _{34})}{_0^{\pi /2}𝑑\beta \mathrm{sin}^2\beta \mathrm{cos}^{3N10}\beta }}`$
$`={\displaystyle \frac{4}{\sqrt{\pi }}}{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{3N6}{2}\right)}{\mathrm{\Gamma }\left(\frac{3N9}{2}\right)}}{\displaystyle _0^{\pi /2}}𝑑\beta \mathrm{sin}^2\beta \mathrm{cos}^{3N10}\beta \varphi (\alpha _{34})`$
$`\widehat{R}_{34}^{(N2)}\varphi (\alpha ),`$ (78)
where the last notation is convenient for repeated use.

To describe three particles in $`𝑑\tau \varphi (\alpha _{13})`$ simultaneously, Jacobi vectors of the standard tree are needed. The distance between particles $`1`$ and $`3`$ is related to the corresponding Jacobi vector

$`\stackrel{}{\eta }_{13}={\displaystyle \frac{1}{\sqrt{2}}}(\stackrel{}{r}_3\stackrel{}{r}_1)={\displaystyle \frac{1}{2}}\stackrel{}{\eta }_{N1}+{\displaystyle \frac{\sqrt{3}}{2}}\stackrel{}{\eta }_{N2},`$ (79)
The hyperangle $`\alpha _{13}`$, associated with the distance between particles $`1`$ and $`3`$, through $`\eta _{13}=r_{13}/\sqrt{2}=\rho \mathrm{sin}\alpha _{13}`$, is then

$`\mathrm{sin}^2\alpha _{13}={\displaystyle \frac{1}{4}}\mathrm{sin}^2\alpha +{\displaystyle \frac{3}{4}}\mathrm{cos}^2\alpha \mathrm{sin}^2\beta +`$
$`{\displaystyle \frac{\sqrt{3}}{2}}\mathrm{sin}\alpha \mathrm{cos}\alpha \mathrm{sin}\beta \mathrm{cos}\theta _x,`$ (80)
where $`\theta _x`$ is the angle between the Jacobi vectors $`\stackrel{}{\eta }_{N1}`$ and $`\stackrel{}{\eta }_{N2}`$. Note that $`\varphi (\alpha _{13})`$, through $`\alpha _{13}`$, for fixed $`\alpha `$ depends on $`\beta `$ and $`\theta _x`$, which leaves a two-dimensional integral. Therefore eq. (75) becomes

$`{\displaystyle 𝑑\tau \varphi (\alpha _{13})}`$ $`=`$ $`{\displaystyle \frac{_0^{\pi /2}𝑑\beta \mathrm{sin}^2\beta \mathrm{cos}^{3N10}\beta _0^\pi 𝑑\theta _x\mathrm{sin}\theta _x\varphi (\alpha _{13})}{_0^{\pi /2}𝑑\beta \mathrm{sin}^2\beta \mathrm{cos}^{3N10}\beta _0^\pi 𝑑\theta _x\mathrm{sin}\theta _x}}`$ (81)
$`=`$ $`{\displaystyle \frac{2}{\sqrt{\pi }}}{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{3N6}{2}\right)}{\mathrm{\Gamma }\left(\frac{3N9}{2}\right)}}{\displaystyle _0^{\pi /2}}𝑑\beta \mathrm{sin}^2\beta \mathrm{cos}^{3N10}\beta {\displaystyle _0^\pi }𝑑\theta _x\mathrm{sin}\theta _x\varphi (\alpha _{13}).`$
This integral can be reduced to one dimension by a partial integration. The final one-dimensional integral becomes

$`{\displaystyle }d\tau \varphi (\alpha _{13})={\displaystyle \frac{4}{\sqrt{3\pi }}}{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{3N6}{2}\right)}{\mathrm{\Gamma }\left(\frac{3N7}{2}\right)}}\mathrm{sin}^1\alpha \mathrm{cos}^{83N}\alpha [{\displaystyle _{(\alpha \pi /3)\mathrm{\Theta }(\alpha >\pi /3)}^{\pi /2|\pi /6\alpha |}}d\alpha _{13}\mathrm{cos}^{3N9}\gamma ^+\mathrm{sin}\alpha _{13}\mathrm{cos}\alpha _{13}\varphi (\alpha _{13})`$
$`{\displaystyle _0^{(\pi /3\alpha )\mathrm{\Theta }(\pi /3>\alpha )}}d\alpha _{13}\mathrm{cos}^{3N9}\gamma ^{}\mathrm{sin}\alpha _{13}\mathrm{cos}\alpha _{13}\varphi (\alpha _{13})]\widehat{R}_{13}^{(N2)}\varphi (\alpha ),`$ (82)
where $`\mathrm{sin}^2\gamma ^\pm =4(\mathrm{sin}^2\alpha +\mathrm{sin}^2\alpha _{13}\mathrm{sin}\alpha \mathrm{sin}\alpha _{13})/3`$, and $`\mathrm{\Theta }`$ is the truth function. The last notation is convenient for repeated use. The integration regions of the integrals are shown in fig. 10. At fixed $`\alpha `$ the range of $`\alpha _{13}`$ in the first integral is over regions I and II, while in the second it is only over region II. For $`N=3`$ the integrands of the two integrals are identical and thus the integration over region II in the first integral cancels the second integral.

## Appendix C Integrals (variational)

We first divide the integrals of eq. (48) into similar terms, then compute them in general, and finally in the zero-range limit.

### C.1 Counting terms

We have to evaluate the double sums of eq. (48) including the potential:

$`{\displaystyle \underset{k<l}{\overset{N}{}}}v_{kl}{\displaystyle \underset{m<n}{\overset{N}{}}}\varphi _{mn}.`$ (83)
Three types of terms occur, due to the fact that we vary the wave function component $`\varphi _{12}^{}`$ in eq. (47): the potential concerning particles $`1`$ and $`2`$, the potential concerning one of the particles $`1`$ or $`2`$ and a third particle and the potential concerning neither particle $`1`$ nor $`2`$, but a third and a fourth particle. We obtain

$`{\displaystyle \underset{k<l}{\overset{N}{}}}v_{kl}=v_{12}+{\displaystyle \underset{l=3}{\overset{N}{}}}v_{1l}+{\displaystyle \underset{l=3}{\overset{N}{}}}v_{2l}+{\displaystyle \underset{k3,l>k}{\overset{N}{}}}v_{kl}`$
$`v_{12}+2(N2)v_{13}+{\displaystyle \frac{1}{2}}(N2)(N3)v_{34},`$ (84)
where the arrow indicates the identity of the terms after integration over all angles except $`\alpha _{12}`$ (analogously to the steps leading up to eq. (II.4.2)). Treating each of these in the quadruple sum, where the repeated use of arrows ($``$) has the meaning given just above:

Fixing $`\varphi _{12}^{}`$ and $`v_{12}`$ yields three different terms:

$`v_{12}{\displaystyle \underset{m<n}{\overset{N}{}}}\varphi _{mn}=`$ (85)
$`v_{12}\left(\varphi _{12}+{\displaystyle \underset{n=3}{\overset{N}{}}}\varphi _{1n}+{\displaystyle \underset{n=3}{\overset{N}{}}}\varphi _{2n}+{\displaystyle \underset{m3,n>m}{\overset{N}{}}}\varphi _{mn}\right)`$
$`v_{12}\left(\varphi _{12}+2(N2)\varphi _{13}+{\displaystyle \frac{1}{2}}(N2)(N3)\varphi _{34}\right),`$
as shown in fig. 11.

Fixing $`\varphi _{12}^{}`$ and $`v_{13}`$ yields seven different terms. These can be identified in two steps, the first of which separates into four different sums:

$`v_{13}{\displaystyle \underset{m<n}{\overset{N}{}}}\varphi _{mn}=`$ (86)
$`v_{13}\left({\displaystyle \underset{n=2}{\overset{N}{}}}\varphi _{1n}+{\displaystyle \underset{n=3}{\overset{N}{}}}\varphi _{2n}+{\displaystyle \underset{n=4}{\overset{N}{}}}\varphi _{3n}+{\displaystyle \underset{m4,n>m}{\overset{N}{}}}\varphi _{mn}\right).`$
Each of these four terms are then identified as:

$`v_{13}`$ $`{\displaystyle \underset{n=2}{\overset{N}{}}}\varphi _{1n}=v_{13}\left(\varphi _{12}+\varphi _{13}+{\displaystyle \underset{n=4}{\overset{N}{}}}\varphi _{1n}\right)`$ (87)
$`v_{13}\left(\varphi _{12}+\varphi _{13}+(N3)\varphi _{14}\right).`$
The similarity of the terms in the first sum becomes apparent when carrying out the integration over $`d\mathrm{\Omega }_{N2}`$, e.g. $`𝑑\mathrm{\Omega }_{N2}v_{13}\varphi _{14}=𝑑\mathrm{\Omega }_{N2}v_{13}\varphi _{15}`$. The rest are found in similar ways:

$`v_{13}{\displaystyle \underset{n=3}{\overset{N}{}}}\varphi _{2n}=v_{13}\left(\varphi _{23}+{\displaystyle \underset{n=4}{\overset{N}{}}}\varphi _{2n}\right)`$
$`v_{13}\left(\varphi _{23}+(N3)\varphi _{24}\right),`$ (88)
$`v_{13}{\displaystyle \underset{n=4}{\overset{N}{}}}\varphi _{3n}v_{13}(N3)\varphi _{34},`$ (89)
$`v_{13}{\displaystyle \underset{m4,n>m}{\overset{N}{}}}\varphi _{mn}v_{13}{\displaystyle \frac{1}{2}}(N3)(N4)\varphi _{45}.`$ (90)
The resulting seven types are shown in fig. 12.

Fixing $`\varphi _{12}^{}`$ and $`v_{34}`$ yields six different terms, identified as follows. The first step is:

$`v_{34}{\displaystyle \underset{m<n}{\overset{N}{}}}\varphi _{mn}=v_{34}({\displaystyle \underset{n=2}{\overset{N}{}}}\varphi _{1n}+`$ (91)
$`{\displaystyle \underset{n=3}{\overset{N}{}}}\varphi _{2n}+{\displaystyle \underset{n=4}{\overset{N}{}}}\varphi _{3n}+{\displaystyle \underset{n=5}{\overset{N}{}}}\varphi _{4n}+{\displaystyle \underset{m5,n>m}{\overset{N}{}}}\varphi _{mn}).`$
In the next step the sums are treated:

$`v_{34}{\displaystyle \underset{n=2}{\overset{N}{}}}\varphi _{1n}=v_{34}\left(\varphi _{12}+\varphi _{13}+\varphi _{14}+{\displaystyle \underset{n=5}{\overset{N}{}}}\varphi _{1n}\right)`$
$`v_{34}\left(\varphi _{12}+2\varphi _{13}+(N4)\varphi _{15}\right),`$ (92)
$`v_{34}{\displaystyle \underset{n=3}{\overset{N}{}}}\varphi _{2n}=v_{34}\left(\varphi _{23}+\varphi _{24}+{\displaystyle \underset{n=5}{\overset{N}{}}}\varphi _{2n}\right)`$
$`v_{34}\left(2\varphi _{13}+(N4)\varphi _{15}\right),`$ (93)
$`v_{34}{\displaystyle \underset{n=4}{\overset{N}{}}}\varphi _{3n}=v_{34}\left(\varphi _{34}+{\displaystyle \underset{n=5}{\overset{N}{}}}\varphi _{3n}\right)`$
$`v_{34}\left(\varphi _{34}+(N4)\varphi _{35}\right),`$ (94)
$`v_{34}{\displaystyle \underset{n=5}{\overset{N}{}}}\varphi _{4n}v_{34}(N4)\varphi _{35},`$ (95)
$`v_{34}{\displaystyle \underset{m5,n>m}{\overset{N}{}}}\varphi _{mn}v_{34}{\displaystyle \frac{1}{2}}(N4)(N5)\varphi _{56}.`$ (96)
See the six types in fig. 13.

### C.2 Evaluation of terms

The term of fig. 11a is trivial since the integrand is independent of $`\tau `$. The terms of figs. 11b 11c, 12a, 12b, 13a, and 13c can be evaluated by eqs. (78) and (82).

The term of fig. 12c becomes with the use of the standard Jacobi tree of fig. 9a

$`{\displaystyle 𝑑\tau f(\alpha _{13})g(\alpha _{23})}={\displaystyle \frac{2}{\sqrt{\pi }}}{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{3N6}{2}\right)}{\mathrm{\Gamma }\left(\frac{3N9}{2}\right)}}{\displaystyle _0^{\pi /2}}𝑑\beta \mathrm{sin}^2\beta \mathrm{cos}^{3N10}\beta {\displaystyle _0^\pi }𝑑\theta \mathrm{sin}\theta f(\alpha _{13})g(\alpha _{23}),`$ (97)
$`\mathrm{sin}^2\alpha _{1(2)\mathrm{\hspace{0.33em}3}}={\displaystyle \frac{3}{4}}\mathrm{cos}^2\alpha \mathrm{sin}^2\beta +{\displaystyle \frac{1}{4}}\mathrm{sin}^2\alpha \pm {\displaystyle \frac{\sqrt{3}}{2}}\mathrm{cos}\alpha \mathrm{sin}\alpha \mathrm{sin}\beta \mathrm{cos}\theta .`$ (98)
The term of fig. 13f becomes with the use of the alternative (12)(34)(56)-tree of fig. 9e

$`{\displaystyle 𝑑\tau v(\alpha _{34})\varphi (\alpha _{56})}={\displaystyle \frac{2A_N}{\pi }}{\displaystyle 𝑑\beta \mathrm{sin}^2\beta \mathrm{cos}^{3N10}\beta 𝑑\gamma \mathrm{sin}^2\gamma \mathrm{cos}^{3N13}\gamma v(\alpha _{34})\varphi (\alpha _{56})},`$ (99)
$`\mathrm{sin}\alpha _{34}=\mathrm{cos}\alpha \mathrm{sin}\beta ,\mathrm{sin}\alpha _{56}=\mathrm{cos}\alpha \mathrm{cos}\beta \mathrm{sin}\gamma ,A_N(3N8)(3N10)(3N12).`$ (100)
The terms of figs. 12g, 13b, and 13d are evaluated using the (123)(45)- and (12)(345)-trees of figs. 9c and 9d, so

$`{\displaystyle 𝑑\tau I_5(\alpha ,\tau )}={\displaystyle \frac{A_N}{\pi }}{\displaystyle _0^{\pi /2}}𝑑\beta \mathrm{sin}^2\beta \mathrm{cos}^{3N10}\beta {\displaystyle _0^{\pi /2}}𝑑\gamma \mathrm{sin}^2\gamma \mathrm{cos}^{3N13}\gamma {\displaystyle _0^\pi }𝑑\theta _{x,z}\mathrm{sin}\theta _{x,z}I_5(\alpha ,\tau ),`$ (101)
where $`I_5(\alpha ,\tau )`$ can be either $`v(\alpha _{34})\varphi (\alpha _{35})`$ or $`f(\alpha _{13})g(\alpha _{45})`$. The relevant angles are $`\mathrm{sin}\alpha _{34}=\mathrm{cos}\alpha \mathrm{sin}\beta `$, $`\mathrm{sin}^2\alpha _{35}=\mathrm{cos}^2\alpha (3\mathrm{cos}^2\beta \mathrm{sin}^2\gamma +\mathrm{sin}^2\beta +2\sqrt{3}\mathrm{cos}\beta \mathrm{sin}\beta \mathrm{sin}\gamma \mathrm{cos}\theta _z)/4`$, $`\mathrm{sin}\alpha _{45}=\mathrm{cos}\alpha \mathrm{cos}\beta \mathrm{sin}\gamma `$, and $`\alpha _{13}`$ given by eq. (B). Note that the integral $`𝑑\tau v(\alpha _{34})\varphi (\alpha _{15})=𝑑\tau \varphi (\alpha _{13})v(\alpha _{45})`$.

The terms of figs. 12d, 12e, 12f, and 13e are evaluated using the standard Jacobi tree. Then eq. (75) reduces to

$`{\displaystyle 𝑑\tau f(\alpha _{13})g(\alpha _{i4})}`$ $`=`$ $`{\displaystyle \frac{A_N}{4\pi ^2}}{\displaystyle _0^{\pi /2}}𝑑\beta \mathrm{sin}^2\beta \mathrm{cos}^{3N10}\beta {\displaystyle _0^{\pi /2}}𝑑\gamma \mathrm{sin}^2\gamma \mathrm{cos}^{3N13}\gamma `$ (102)
$`{\displaystyle _0^\pi }𝑑\theta _x\mathrm{sin}\theta _x{\displaystyle _0^\pi }𝑑\theta _y\mathrm{sin}\theta _y{\displaystyle _0^{2\pi }}𝑑\phi f(\alpha _{13})g(\alpha _{i4});i=1,2,3.`$
The angles $`\alpha _{ij}`$ can be determined by $`\rho \mathrm{sin}\alpha _{ij}=\eta _{ij}`$ through the relations

$`\stackrel{}{\eta }_{13}={\displaystyle \frac{\sqrt{3}}{2}}\stackrel{}{\eta }_{N2}+{\displaystyle \frac{1}{2}}\stackrel{}{\eta }_{N1},\stackrel{}{\eta }_{14}=\sqrt{{\displaystyle \frac{2}{3}}}\stackrel{}{\eta }_{N3}+{\displaystyle \frac{1}{2\sqrt{3}}}\stackrel{}{\eta }_{N2}+{\displaystyle \frac{1}{2}}\stackrel{}{\eta }_{N1},`$ (103)
$`\stackrel{}{\eta }_{24}=\sqrt{{\displaystyle \frac{2}{3}}}\stackrel{}{\eta }_{N3}+{\displaystyle \frac{1}{2\sqrt{3}}}\stackrel{}{\eta }_{N2}{\displaystyle \frac{1}{2}}\stackrel{}{\eta }_{N1},\stackrel{}{\eta }_{34}=\sqrt{{\displaystyle \frac{2}{3}}}\stackrel{}{\eta }_{N3}{\displaystyle \frac{1}{\sqrt{3}}}\stackrel{}{\eta }_{N2}.`$ (104)
### C.3 Results in the $`\delta `$-limit

The integrals in the short-range limit, when the range $`b`$ of $`V(r_{ij})`$ is much smaller than the size scale $`\rho `$, are:

$`{\displaystyle 𝑑\tau v(\alpha _{34})\varphi (\alpha _{13})}v_1(\alpha )\widehat{R}_{3413}^{\text{(2)}}\varphi (\alpha ),`$ (105)
$`{\displaystyle 𝑑\tau v(\alpha _{34})\varphi (\alpha _{15})}v_1(\alpha )\widehat{R}_{13}^{(N3)}\varphi (\alpha ),`$ (106)
$`{\displaystyle 𝑑\tau v(\alpha _{34})\varphi (\alpha _{34})}=\widehat{R}_{34}^{(N2)}v\varphi (\alpha )v_1(\alpha )\varphi (0),`$ (107)
$`{\displaystyle 𝑑\tau v(\alpha _{34})\varphi (\alpha _{35})}v_1(\alpha )\widehat{R}_{3435}^{\text{(1)}}\varphi (\alpha ),`$ (108)
$`{\displaystyle 𝑑\tau v(\alpha _{34})\varphi (\alpha _{56})}v_1(\alpha )\widehat{R}_{34}^{(N3)}\varphi (\alpha ),`$ (109)
$`{\displaystyle 𝑑\tau v(\alpha _{13})\varphi (\alpha _{13})}=\widehat{R}_{13}^{(N2)}v\varphi (\alpha )v_2(\alpha )\varphi (0),`$ (110)
$`{\displaystyle 𝑑\tau v(\alpha _{13})\varphi (\alpha _{i4})}v_2(\alpha )\widehat{R}_{1314}^{\text{(2)}}\varphi (\alpha );i=1,3,`$ (111)
$`{\displaystyle 𝑑\tau v(\alpha _{13})\varphi (\alpha _{23})}v_2(\alpha )\varphi (\alpha ),`$ (112)
$`{\displaystyle 𝑑\tau v(\alpha _{13})\varphi (\alpha _{24})}v_2(\alpha )\widehat{R}_{1324}^{\text{(2)}}\varphi (\alpha ),`$ (113)
$`{\displaystyle 𝑑\tau v(\alpha _{13})\varphi (\alpha _{45})}v_2(\alpha )\widehat{R}_{1345}^{\text{(1)}}\varphi (\alpha ).`$ (114)
The integrals are given by

$`\widehat{R}_{ijkl}^{\text{(1)}}\varphi (\alpha ){\displaystyle \frac{4}{\sqrt{\pi }}}{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{3N9}{2}\right)}{\mathrm{\Gamma }\left(\frac{3N12}{2}\right)}}{\displaystyle _0^{\pi /2}}𝑑\gamma \mathrm{sin}^2\gamma \mathrm{cos}^{3N13}\gamma \varphi (\alpha _{kl}^0),`$ (115)
where $`\mathrm{sin}\alpha _{35}^0\sqrt{3}\mathrm{cos}\alpha \mathrm{sin}\gamma /2`$, $`\mathrm{sin}\alpha _{45}^0\mathrm{cos}\alpha \mathrm{cos}\beta _0\mathrm{sin}\gamma `$, $`\mathrm{sin}\beta _0\mathrm{tan}\alpha /\sqrt{3}`$, and

$`\widehat{R}_{ijkl}^{\text{(2)}}\varphi (\alpha ){\displaystyle \frac{2}{\sqrt{\pi }}}{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{3N9}{2}\right)}{\mathrm{\Gamma }\left(\frac{3N12}{2}\right)}}{\displaystyle _0^{\pi /2}}𝑑\gamma \mathrm{sin}^2\gamma \mathrm{cos}^{3N13}\gamma {\displaystyle _0^\pi }𝑑\theta _x\mathrm{sin}\theta _x\varphi (\alpha _{kl}^0),`$ (116)
$`\mathrm{sin}^2\alpha _{14}^0{\displaystyle \frac{1}{9}}\mathrm{sin}^2\alpha +{\displaystyle \frac{2}{3}}\mathrm{cos}^2\alpha \mathrm{cos}^2\beta _0\mathrm{sin}^2\gamma +{\displaystyle \frac{2\sqrt{2}}{3\sqrt{3}}}\mathrm{sin}\alpha \mathrm{cos}\alpha \mathrm{cos}\beta _0\mathrm{sin}\gamma \mathrm{cos}\theta _x,`$ (117)
$`\mathrm{sin}^2\alpha _{24}^0{\displaystyle \frac{4}{9}}\mathrm{sin}^2\alpha +{\displaystyle \frac{2}{3}}\mathrm{cos}^2\alpha \mathrm{cos}^2\beta _0\mathrm{sin}^2\gamma +{\displaystyle \frac{4\sqrt{2}}{3\sqrt{3}}}\mathrm{sin}\alpha \mathrm{cos}\alpha \mathrm{cos}\beta _0\mathrm{sin}\gamma \mathrm{cos}\theta _x,`$ (118)
$`\mathrm{sin}^2\alpha _{13}^0{\displaystyle \frac{1}{4}}\mathrm{sin}^2\alpha +{\displaystyle \frac{1}{2}}\mathrm{cos}^2\alpha \mathrm{sin}^2\gamma +{\displaystyle \frac{1}{\sqrt{2}}}\mathrm{sin}\alpha \mathrm{cos}\alpha \mathrm{sin}\gamma \mathrm{cos}\theta _x.`$ (119)
$`\widehat{R}_{ijkl}^{\text{(2)}}\varphi (\alpha )`$ appears as a two-dimensional integral but can be reduced to a one-dimensional integral, analogously to eq. (82).