The Faddeev-Yakubovsky symphony

TL;DR

This paper summarizes the derivation of the Faddeev-Yakubovsky equations in configuration space for systems of three to five interacting particles, providing a foundational overview of their formulation.

Contribution

It offers a concise summary of the main steps involved in deriving the Faddeev-Yakubovsky equations for N=3, 4, and 5 particles, clarifying their development.

Findings

Clarified the derivation process of Faddeev-Yakubovsky equations.

Provided a structured overview for N=3, 4, 5 particle systems.

Facilitated understanding of multi-particle quantum systems.

Abstract

We briefly summarize the main steps leading to the Faddeev-Yakubovsky equations in configuration space for N=3, 4 and 5 interacting particles.

Tables

Table 1. Table 1: Scalability of the Faddeev-Yakubovsky scheme as a function of the number of interacting particles N.: ne is the number of equations,(equal to the number of FY components), neid the same number in case of identical particles, n a subscript 𝑛 𝑎 n_{a} the number of amplitudes in the partial wave expansion.

N	Ne	Neid	$n_a$
2	1	1	2
3	3	1	$\sim {10}^2$
4	18	2	$\sim {10}^4$
5	180	5	$\sim {10}^6$
6	2700	15
	…
N	$\frac{N!\left(N-1\right)!}{2^{N-1}}$	Int$\left(\frac{2\left(N-1\right)!}{{\left(\pi /2\right)}^N}\right)$

Equations

$$(E-{\mathcal{H}}_0)\psi =V\psi {\mathcal{H}}_0=-\sum _{i=1}^N\frac{{\hslash }^2}{2m_i}{\Delta }_{{\vec{r}}_i}V=\sum _{i<j=1}^N{V}_{ij}=V_{12}+V_{23}+V_{23}$$

$$\Psi ({\vec{r}}_1,{\vec{r}}_2,{\vec{r}}_3)\sim {\chi }_e\left(\frac{{\vec{r}}_2+{\vec{r}}_3}{2}-{\vec{r}}_1\right)u_e({\vec{r}}_2-{\vec{r}}_3),$$

$$\Psi ({\vec{r}}_1,{\vec{r}}_2,{\vec{r}}_3)\sim {\chi }_i\left(\frac{{\vec{r}}_2+{\vec{r}}_3}{2}-r_1\right)u_i^{\ast }({\vec{r}}_2-{\vec{r}}_3),$$

$$\Psi ({\vec{r}}_1,{\vec{r}}_2,{\vec{r}}_3)\sim {\chi }_r\left(\frac{{\vec{r}}_1+{\vec{r}}_3}{2}-{\vec{r}}_2\right)u_r^{\ast }({\vec{r}}_1-{\vec{r}}_2),$$

$$\Psi ({\vec{r}}_1,{\vec{r}}_2,{\vec{r}}_3)\sim \frac{e^{ik\rho }}{\sqrt{\rho }};{\rho }^2=({\vec{r}}_2-{\vec{r}}_3{)}^2+{\left(\frac{{\vec{r}}_2+{\vec{r}}_3}{2}-{\vec{r}}_1\right)}^2.$$

$$T=V+V{G}_0(E)T,G_0(E)=\frac{1}{E-H_0+iϵ},$$

$$V_{23}({\vec{k}}_1,{\vec{p}}_1,{\vec{k}}_1^{\prime },{\vec{p}}_1^{\prime })=v_{23}({\vec{k}}_1,{\vec{k}}_1^{\prime })\delta ({\vec{p}}_1-{\vec{p}}_1^{\prime }),$$

$$T$$

$$+$$

$${\vec{x}}_{\alpha }$$

$${\vec{y}}_{\alpha }$$

$$M\vec{R}=m_1{\vec{r}}_1+m_2{\vec{r}}_2+m_3{\vec{r}}_3;M=m_1+m_2+m_3.$$

$$\psi ({\vec{r}}_1,{\vec{r}}_2,{\vec{r}}_3)=\Psi ({\vec{x}}_i,{\vec{y}}_i)e^{i\vec{P}\cdot \vec{R}}.$$

$$[E-H_0]\Psi ({\vec{x}}_i,{\vec{y}}_i)=V\Psi ({\vec{x}}_i,{\vec{y}}_i),$$

$$H_0=-({\Delta }_{{\vec{x}}_i}+{\Delta }_{{\vec{y}}_i}).$$

$${\vec{x}}_{\alpha }$$

$${\vec{y}}_{\alpha }$$

$$c_{\alpha \beta }=-\sqrt{\frac{{\mu }_{\alpha }{\mu }_{\beta }}{(1-{\mu }_{\beta })(1-{\mu }_{\alpha })}},s_{\alpha \beta }={ϵ}_{\alpha \beta }\sqrt{1-c_{\alpha \beta }^2},$$

$$\Psi =\sum _{i<j=1}^N{\Phi }_{ij}={\Phi }_{12}+{\Phi }_{13}+{\Phi }_{23},$$

$${\Phi }_{ij}=G_0(E)V_{ij}\Psi ;G_0(E)=(E-H_0{)}^{-1},$$

$$(E-H_0-V_{ij}){\Phi }_{ij}=V_{ij}\sum _{kl\ne ij}{\Phi }_{kl},$$

$$\begin{array}[]{l}\Phi_{1}(x_{1},y_{1})\cr\Phi_{2}(x_{2},y_{3})\cr\Phi_{3}(x_{3},y_{3})\end{array}=\begin{array}[]{lclrr}u_{1}(x_{1})&\left[\sin{q_{1}y}\right.&+&\tan(\delta)\cos{q_{1}y_{1}}&\left.\right]/(x_{1}y_{1}),\cr&0,&&&\cr&0,&&&\cr\end{array}$$

$$\begin{array}[]{l}\phantom{\sqrt{q_{1}\over q_{2}}K_{12}\cos{q_{2}y_{2}}}\Phi_{1}(x_{1},y_{1})\cr\phantom{\sqrt{q_{1}\over q_{2}}K_{12}\cos{q_{2}y_{2}}}\Phi_{2}(x_{2},y_{3})\cr\phantom{\sqrt{q_{1}\over q_{2}}K_{12}\cos{q_{2}y_{2}}}\Phi_{3}(x_{3},y_{3})\end{array}=\begin{array}[]{lclrr}u_{1}(x_{1})&\left[\sin{q_{1}y}\right.&+&K_{11}\cos{q_{1}y_{1}}&\left.\right]/(x_{1}y_{1}),\cr u_{2}(x_{2})&&&\sqrt{q_{1}\over q_{2}}K_{12}\cos{q_{2}y_{2}}/(x_{2}y_{2}),&\cr u_{3}(x_{3})&&&\sqrt{q_{1}\over q_{3}}K_{13}\cos{q_{3}y_{3}}/(x_{3}y_{3}).&\cr\end{array}$$

$${\Phi }_i({\vec{x}}_i,{\vec{y}}_i)=\sum _{{\alpha }_i}\frac{{\phi }_{{\alpha }_i}(x_i,y_i)}{x_i{y}_i}{Y}_{{\alpha }_i}({\hat{x}}_i,{\hat{y}}_i),$$

$$Y_{{\alpha }_i}({\hat{x}}_i,{\hat{y}}_i)=\sum _{m_x{m}_y}<l_x{m}_x;l_y{m}_y|l_x{l}_x;LM>Y_{l_x{m}_x}({\hat{x}}_i)Y_{l_y{m}_y}({\hat{y}}_i),$$

$$\left[q^2+{\partial }_x^2+{\partial }_y^2-v_{\alpha }^e(x,y)\right]{\phi }_{\alpha }(x,y)$$

$$v_{\alpha }^e(x,y)$$

$$(q^2+{\partial }_x^2+{\partial }_y^2-v_{\alpha }){\varphi }_{\alpha }(x,y)=v_{\alpha }\sum _{\beta =1}^{n_a}{c}_{\alpha \beta }{\int }_{-1}^{+1}du\frac{xy}{x^{\prime }{y}^{\prime }}{\varphi }_{\beta }(x^{\prime },y^{\prime }),$$

$$2{x^{\prime }}^2(x,y,u)$$

$$2{y^{\prime }}^2(x,y,u)$$

$$(E-H_{0})\Psi=V\Psi;\qquad V=\sum_{i<j=1}^{4}V_{ij}=V_{12}+V_{13}+V_{14}+V_{23}+V_{24}+V_{34},\$$

$$H_0=-\left({\Delta }_{\vec{x}}+{\Delta }_{\vec{y}}+{\Delta }_{\vec{z}}\right),$$

Full text

1 Introduction

The Faddeev-Yakubovsky (FY) equations constitute a rigorous mathematical formalism for solving the non relativistic N-body problem of Quantum Mechanics. They were motivated by the drawbacks of the fundamental Schrödinger equation in determining the physical solutions, mainly those related with the scattering states, of a 3-body problem. FY equations provide a framework enabling us to implement additional constrains to the solutions of the Schrödinger equation, which guarantee their physical meaning.

FY equations should in principle allow us to obtain the exact solutions of this fundamental equation for an arbitrary number of interacting particles with the only input of the interparticle potential. The term ”exact” must be here understood in the numerical sense, i.e. allowing numerical methods that could provide accurate enough solutions in a controlled converging scheme. Only when such a possibility is ensured one would be able to disentangle the consequences of a deficient interaction from the consequences of an approximate solution, as was the main aim in the historically first physical problem at the origin of these theoretical developments: the three-nucleon problem. In practice, and more than 50 years after the publication of Faddeev’s seminal work Fad_JETP39_1960 ; Fad_61 ; Fad_65 ; Faddeev_Merkuriev_QSTh_1993 ; Yakubovsky_SJNP5_1967 , the number of interacting particles for which the FY equations have been properly solved remains severely limited to N$\leq$5.

The formulation of FY equations can be considered as a starting point of a very active domain of theoretical physics, the so called Few-Body Physics, where ”Few” emphasizes that obtaining exact solutions is a very difficult task and one might never outreach systems with small number of particles, when aiming the complete solution of the Quantum Mechanical problem. That is if we include the full complexity of the scattering solutions.

We will present in this contribution a brief summary of the main steps in the formulation of these fundamental equations, essentially in configuration space, starting from the pioneering work of Faddeev up to the recent solutions of the 5-body problem Rimas_5N_FBS_2018 ; Rimas_5N_PRC97_2018 ; LHC_PLB791_2019 . Our main focus will be set to explain the generalisation procedure when the number of interacting particles is increased, describing in parallel an efficient and unified numerical method and pointing out the natural limitations of this approach.

2 Shortcomings of Schrödinger’s equation in describing the 3-body scattering

The problem faced by the Schrödinger’s equation to describe $N\geq 3$ interacting particle has been realized since the very beginning of Quantum Mechanics, can be summarized in the two following complementary aspects.

When formulated in configuration space, the Schrödinger equation for N=3 particles interacting via pairwise interactions $V_{ij}$

[TABLE]

is a second order partial derivative equation for a function $\psi(\vec{r}_{1},\vec{r}_{2},\vec{r}_{3})$ defined on $\mathbf{R}^{9}$. Its solution is physical only once $\psi$ satisfy some appropriate boundary conditions (BC).

However, the scattering of one particle on a 2-body bound state, say 1+(23) reaction, may give rise to many asymptotic channels, simultaneously coexisting, which have different behaviours. It turns to be impossible to implement the ensemble of these complex asymptotic structures in a single function. One can expect, for instance:

• •

In the elastic channel $1+(23)\to 1+(23)$ an asymptotic behavior like

[TABLE]

where $u_{e}$ is the bound state wave function of particle (23).

Depending on the incoming energy, this channel coexists with as many inelastic channels (eventually infinite) as excited states of the (23) subsystem are accessible at this energy. All of them will have similar, but different, boundary conditions like

[TABLE]

where $u^{*}_{i}$ represents an excited state wave function.

• •

In the $1+(23)\to 3+(12)$ rearrangement channel, one expects a solution having the form

[TABLE]

where $u^{*}_{r}$ is the bound state wave function of particles (12).

Corresponding expressions are expected also for the $1+(23)\to 2+(13)$ rearrangement channels.

• •

For the $1+(23)\to 1+2+3$ break-up channel one could guess

[TABLE]

This list still misses electromagnetic processes like radiative capture as well as relativistic effects like creation or annihilation of particles.

In a consistent theory, it would be highly desirable (i) to deduce these different asymptotic behaviors from the equation itself and (ii) uncouple and impose them naturally as a boundary condition in the corresponding regions of configuration space. This is however not possible with eq. (1). It is furthermore clear, that the Schrödinger equation is not flexible enough to account and treat such a rich variety of physical situations which are present already for N=3.

On another hand, when formulated in momentum space, the Schrödinger equation, as well as its equivalent Lippmann-Schwinger (LS) equation for the 3-body T-matrix,

[TABLE]

does not need any kind of BC. However the potential $V$ in momentum space is plagued with delta-like singularities on each parwise term $V_{ij}$, like for instance

[TABLE]

which propagate in all orders of the perturbative expansion of (6)

[TABLE]

thus breaking the compactness of the integral kernel Fad_65 ; Faddeev_Merkuriev_QSTh_1993 . Because of that, the mathematical foundations of the theory, mainly the existence and uniqueness of the solutions, are lost and we are left with an inconsistent ”theoretical no man’s land”.

3 Faddeev equations

In order to edge these problems, L.D. Faddeev derived in 1960s a set of equations, equivalent to (1), which constitutes a rigorous mathematical framework for describing the variety of physical situations involved in N=3. This seminal result was later used by S.P. Merkuriev to derive the boundary conditions in configuration space Merkuriev_71 ; Merkuriev_74 what allowed the first solutions in the case of 3-nucleon problem with realistic interactions GL_PRL29_1972 ; GLM_PRL33_1974 ; MGL_AP99_1976 .

A preliminary step consists in isolating the intrinsic dynamics of the three-body Hamiltonian (1). This is achieved by introducing the Jacobi coordinates

[TABLE]

where $(\alpha\beta\gamma)$ is a circular permutation of (123), $\mu_{s,t}$ denotes the reduced mass of the system formed by particle clusters $s$ and $t$, $\vec{R}_{s}$ the center of mass of cluster $s$. An arbitrary mass $m_{0}$ is introduced to fix the length scale. They are supplemented with the center of mass coordinate

[TABLE]

In terms of them, the solution of (1), $\psi$ factorizes into an intrinsic part $\Psi$ and a plane wave for the center of mass motion

[TABLE]

The intrinsic wave function $\Psi$, the only one in which we will be interested hereafter, is a solution of

[TABLE]

with (in ${\hbar^{2}/m_{0}}$ units)

[TABLE]

This trivial geometrical operation not only reduces the dimensionality of the problem, from on $\mathbf{R}^{9}$ to $\mathbf{R}^{6}$, but is the only way to properly disentangle the intrinsic 3-particle energy the from its center of mass excitations. Notice that there are 3 different sets of Jacobi coordinates (12) and that the form of intrinsic free Hamiltonian (16) is independent of a particular choice. They are related each other by orthogonal transformations

[TABLE]

with

[TABLE]

and $\epsilon_{\alpha\beta}=(-1)^{\beta-\alpha}sig{(\beta-\alpha)}$.

The central idea of Faddeev was to split the wave function in a sum of three terms, the so-called Faddeev components (FC),

[TABLE]

defined by

[TABLE]

and so each of them associated with an interacting pair $V_{ij}$.

The ensemble of FC obey the system of three coupled equations

[TABLE]

known as Faddeev equations for the 3-body problem, and provide a solution of the corresponding Schrödinger equation. Its is indeed straightforward to see that if $\Phi_{ij}$ represent a solution of eqs. (22), their sum (20) provides a solution of eq. (1).

Notice that, although using relations (19), each Faddeev component can be expressed in any of the Jacobi coordinate sets ($x_{\alpha},y_{\alpha}$), it has a natural expression in the one on which the corresponding pair potential has the simplest form, $V_{ij}(x_{k})$. This justifies the often used one-index notation for the components $\Phi_{i}\equiv\Phi_{jk}$, valid only for N=3.

The main advantage of the Faddeev equations is that they provide three functions, on which we can impose the proper boundary conditions corresponding to the different channels. In the asymptotic region of the configuration space ($x_{k}\to\infty$) the three Faddeev equations, coupled by short range interactions $V_{ij}(x_{k})$, decouple from each other. This fact allows to associate to each FC a well defined asymptotic structure which can be used as boundary conditions for the system of differential equations (22). For instance for the 1+(23)$\to$ 1+(23) S-wave elastic scattering one can prove MGL_AP99_1976 that

[TABLE]

whereas when all rearrangements channels are open one must impose

[TABLE]

We have used here the one-index notation and denoted by $q_{i}$ the momentum of each channel.

More interestingly, and this was the main Faddeev’s result, the T-matrix perturbative expansion (10) can be reordered in such a way that all delta functions disappear leading to compact equations.

To solve in practice equations (22) one can perform a partial wave expansion of each FC in its proper set of Jacobi coordinates

[TABLE]

where $Y_{\alpha_{i}}$ denote the bipolar spherical harmonics

[TABLE]

and $\alpha_{i}=\{l_{x}l_{y}LM\}$ the set of intermediate quantum numbers in the angular (and eventually spin and isospin) couplings. After inserting (25) in (22) and projecting the angular part, one is led with a system of two-dimensional coupled integrodifferential equations having the the form

[TABLE]

with the effective potential

[TABLE]

The integral kernels $h^{(3)}_{\alpha\beta}$ constitutes the key ingredient of the calculation. Their precise expressions can be found in Gloeckle_Book_83 ; 4B_h .

Faddeev equations have a very simple form for 3 identical particles in the so called S-wave approximation, i.e. where all the angular momenta are set to zero:

[TABLE]

with $\alpha,\beta=1,\ldots,n_{a}$,

[TABLE]

and $c_{\alpha\beta}$ are numerical coefficients which take the values: $n_{a}$=1 and $c$=2 for a 3 boson system with total angular momentum L=0, $n_{a}$=1 and $c$=-1 for a 3 fermion system with total spin J=3/2, $n_{a}$=2 and $c_{11}$=$c_{22}$=1/2, $c_{12}$=$c_{21}$=3/2 for a 3 fermion system with total spin J=1/2.

Before closing this section, several remarks are in order:

We have conscientiously ignored until now, the existence of three-body terms in the potential. They play an important role in the description of nuclear systems but are not relevant in our discussion. The interested reader can found a sound presentation in Gloeckle_Book_83 . 2. 2.

One of the essential properties of the Faddeev equations in view of fixing the boundary condition their decoupling in the asymptotic regions of the configuration space. This is ensured under the assumption of short range pairwise interactions. However in presence of Coulomb forces this decoupling is not guaranteed, at least in the above presented original form, and the whole formalism could seem questionable. This problem was solved by Merkuriev Merkuriev_80 ; Merkuriev_81 by introducing an artificial splitting of the Coulomb potential into a short and long range parts by means of a smooth cut-off function. The long range part is kept in the left hand side of the Faddeev equations to define the Coulomb asymptotes and the short range parts appear in the right hand side ensuring the decoupling of the equations. Merkuriev approach has proven to be very successful in several nuclear physics problems like p-d reactions Papp as well as handling purely atomic problems KCG_PRA46_1992 ; Valdes ; Yakovl . 3. 3.

In parallel with Faddeev work, P.H. Noyes and collaborators had independently proposed, and properly solved for the bound as well as for elastic scattering cases, the differential form of the same equations Noyes_Fideldey_1968 . They are sometimes denoted as Faddeev-Noyes equations, specially in their simplified S-wave form (28). 4. 4.

It is worth noticing that, even for the three-body case, there remain many unresolved problems, for instance those related with the presence of an infinite number of the open asymptotic channels. In particular, the very challenging problem of atomic anti-hydrogen production slow antiproton collisions with hydrogen atoms due to the presence of large number of the open asymptotic channels.

The Faddeev solution of the quantum mechanical 3-body problem was of paramount importance in theoretical physics and definitely set the foundations for an ab initio solution of the many particle systems.

Nevertheless it remains a very particular case which somehow hides the complexity of the N-body problem as well as the road for its generalization to larger N. In the Faddeev equations, 3 is the number of particles, the number of interacting pairs, the number of Jacobi sets, as well as the number of ways for breaking (or building!) a 3 body cluster (123) into (or from) a lower rank sub-clusters Faddeev_Birmingham_1969 . This later characteristic, of particular relevance in the N-body case, can indeed occur according to the 3 different ”partitions” (12)3,(13)2,(23)1 which have in fact the same topological properties, i.e. that can be related to each other by permutation operators. It is convenient to represent these different topological types in form of ”trees” diagrams, which in the N=3 case reduces to the single one represented in Fig. 1.

Thus, in this case, the number of particles encodes in fact very different properties and make difficult to disentangle the different role they play in the general theory. This is however not the case in general as it will be illustrated in the following sections.

4 Faddeev-Yakubovsky equations for N=4

Few years after Faddeev’s pioneer work, Yakubovsky Yakubovsky_SJNP5_1967 proposed a consistent scheme to build a set of equations which should allow the solution of the general N-body problem. His demonstration was based on induction with respect to the number of particles and this requires the previous solution of the N-1, N-2, …, 2 problems. After several unfruitful attempts Weinberg_63_I_II , Yakubovsky’s result was a real tour de force to get rid of the $\delta$-like singularities in the perturbative expansion of the N-body T-matrix, which disappear after N-2 iterations of the proposed equations.

Rather than reproducing this general approach, which remain quite abstruse even after some pedagogical efforts NY_1980 , we will detail in this section – without any formal demonstration – the main steps in deriving the differential equations for the N=4 case.

The starting point is the intrinsic Schrödinger’s equation for N=4:

[TABLE]

where the 4-body free Hamiltonian $H_{0}$ has the same spherical form (16) as for N=3 case

[TABLE]

expressed in some set of Jacobi coordinates to be detailed later.

The corresponding Faddeev equations (21) and (22) can be written in close analogy with the N=3 case, in terms of the wave function partition

[TABLE]

However, contrary to N=3, an additional decomposition of each component (30) is needed to represent all different asymptotes for the non interacting particles. For instance $\Phi_{12}$, associated with the (12) interacting pair, and fulfilling the Faddeev equation

[TABLE]

is split into 3 components belonging to two different types

[TABLE]

defined in terms of the interacting 2-body Green function ”embedded in the 4-body space”

[TABLE]

by

[TABLE]

They satisfy

[TABLE]

By repeating this procedure to each of the six Faddeev components (30)

[TABLE]

one obtains the set of 18 Faddeev-Yakubovsky (FY) components defined by

[TABLE]

and satisfying a system of 18 coupled equations

[TABLE]

These equation were first formulated in momentum space by Yakubovsky Yakubovsky_SJNP5_1967 , their differential form with the corresponding boundary conditions was elaborated some time later by Merkuriev and collaborators in MY_1982 ; MY_1982_1983 ; MYG_NPA_1984 . Notice that in terms of the total wave function the FY components read

[TABLE]

Each FY component of $\Phi_{ij,k}^{l}$ or $\Phi_{ij,kl}$ type111Usually denoted by $K_{ij,k}^{l}\equiv\Phi_{ij,k}^{l}$ and $H_{ij,kl}\equiv\Phi_{ij,kl}$ corresponds to a topologically different way of breaking a 4-body cluster into its subsystems up to the point where a single interacting pair remains unbroken. These sequences of breaking an N-body cluster – denoted in the literature ”complete partition chains” – are illustrated in Fig. 2 for the pair (12).

The upper left figure represents - from bottom to the top in chronological order - a ”tree” where a four particle system (1234) by separating the particle 4 is first broken into 3+1 structure [(123)4]. Then cluster (123) is decomposed by separating the particle 3, thus building the 4-particle structure [(12)34]. This corresponds to the K-like FY component $\Phi_{12,3}^{4}$ and to the partition chain (1234)$\supset$(123)4$\supset$(12)34.

In the upper right panel, an alternative partition strategy is demonstrated, where a four particle system (1234) is first broken by separating pairs (12) and (34), getting [(12)(34)] partition. In the next step the (34) pair is decomposed [(12)34]. It corresponds to the H-like FY component $\Phi_{12,34}$ and to the partition chain: (1234)$\supset$(12)(34)$\supset$(12)34.

.

By permuting particle indexes, one obtains 6x2=12 K-like components and 3x2=6 H-like, which give the 18 FY coupled equations (44). In case of identical particles permutation symmetry allows to express all 12 K-like components from any of them (idem for H-like components); thus one may reduce the problem to two coupled equations for the two FYC of fig Fig. 2. The Jacobi coordinates, associated with each FY component ($\vec{x}_{K},\vec{y}_{K},\vec{z}_{K}$) or ($\vec{x}_{H},\vec{x}_{H},\vec{x}_{H}$), are schematically represented in the lower panel of Fig. 2. They are denoted generically by $\vec{q}$ and have the general form given in (12), i.e.

[TABLE]

where $\vec{R}_{i}$ denotes the center of mass positions of the two clusters (eventually single particles) connected by each coordinate and $\mu_{st}$ their reduced mass.

A practical solution of the FY equations can be achieved by means of a partial wave expansion of the FY components in analogy with the 3-body case (25). One ends with a system of tridimensional integrodifferential equations of the form

[TABLE]

with $\Delta_{x,y,z}=\partial^{2}_{x}+\partial^{2}_{y}+\partial^{2}_{z}$ and the effective potential

[TABLE]

and $h^{(4)}_{\alpha\beta}$ some integral kernels which can be found in 4B_h . For the case of 4 identical bosons in the S-wave approximation former equations take the simple form MYG_NPA_1984

[TABLE]

where

[TABLE]

During the last twenty years these equations have been used to solve very diverse bound state and scattering problems related with nuclear or atomic physics, including the break-up reactions 4B_Nogga ; 4B_Uzu ; 4B_Deltuva ; 4B_Laza . The interested reader can find in these references a detailed description of the proper techniques to implement the boundary conditions and the associated numerical tools.

5 Faddeev-Yakubovsky equations for N=5

As in the previous cases, the solution of the corresponding Schrödinger equations starts with the Faddeev-like components associated to each interacting pair

[TABLE]

according to (21) and (22). The number of interacting pairs, equal to the number of Faddeev equations (22), is now 5x4/2=10.

Each of these ten Faddeev components $\Psi_{ij}$ is decomposed into six 4-body like Yakubovsky components

[TABLE]

defined, as in the N=4 case (41), by

[TABLE]

In its turn, the 4-body components are further decomposed in terms of five independent 5-body FY components – denoted following Rimas_5N_FBS_2018 ; Rimas_5N_PRC97_2018 by K,H,T,S,F – in the following way:

[TABLE]

with

[TABLE]

Each of these 5-body FY components corresponds to a topologically independent way of breaking a 5-body cluster. At each step one of the ”interacting” clusters is decomposed into two pieces, giving rise to a corresponding FY equation. The procedure is repeated until there remains only a single ”interacting” pair. The possible partition chains for a 5-body system are represented in the upper part of Fig. 3. On the left side one may identify K-like and H-like structures, which are associated with 4+1 particle structures. The three components on the right hand side (T-like, S-like and F-like) represent 3+2 particle structures.

The lower part of Fig. 3 represents the Jacobi coordinates associated to each FY component. They have the general form (47) and are denoted by

[TABLE]

By permuting the particle index one can see that there are 60 different K-type amplitudes, and 30 for each of the other types: H,T,S,F. Finally, this result into the set of 180 FY equations, first derived in Sasakawa_5b ,

[TABLE]

with $\psi$’s given by (58).

Each FY component $F_{a}$=(K,H,T,S,F) is expressed in its own Jacobi set and take now the values on ${\mathbf{R}}^{12}$ and expanded in spherical harmonics for each angular variable:

[TABLE]

where $Y_{\alpha}$ is a generalized ”quadripolar harmonic” accounting for the angular momentum, and eventually spin and isospin, couplings

[TABLE]

and $\alpha$ labels the set of quantum numbers involved in the intermediate couplings.

After projecting the angular part, one is left with a set of coupled four-dimensional integro-differential equations for the reduced radial amplitudes $f_{a,\alpha}$ of the form

[TABLE]

with $\Delta_{xyzw}=\partial^{2}_{x}+\partial^{2}_{y}+\partial^{2}_{z}+\partial^{2}_{w}$ and the effective potential

[TABLE]

These equations have been solved for the first time in a recent work devoted to study the n-4He scattering Rimas_5N_PRC97_2018 , more recently they were applied to compute the complex energies of the ${}^{5}H$ resonant states LHC_PLB791_2019

6 Numerical methods

In the last years we have developed a numerical protocol to solve the FY equations in configuration space for N=3,4 and 5 particles. For the sake of simplicity we will illustrate it in the case N=3; the close analogy we kept in the final equations (27), (49) and (67), make the generalization rather straightforward.

In N=3, we wish to determine on a two dimensional domain $D=[0,x_{n_{x}}]\times[0,y_{n_{y}}]$ the partial wave amplitudes of the FY components which are solution of eq. (27). Our unique approximation is the assumption that the solution we are looking for can be locally expanded in terms of some polynomial basis:

[TABLE]

We used two kind of such bases: the so called splines $f_{i}(x)\equiv S_{i}(x)$ (cubic or quintic) and the Lagrange interpolating functions $f_{i}(x)\equiv L_{i}(x)$. In what follows we will particularise the case of cubic splines. The interested reader can find a detailed explanation of both choices in Sect. 2.9 of Ref. Rimas_DHDR_2019 .

The cubic spline functions $S_{i}(q)$ are associated to each variable $q=x,y$, they are constructed upon a ($n_{q}$+1)-point grid $G_{q}=\{q_{0},q_{1},\ldots,\ldots q_{n}\}$ defined on each of the intervals $[q_{0}=0,q_{n}]$ of the resolution domain $D$.

Two splines are associated to each grid point $q_{i}$: $S_{2i}$ and $S_{2i+1}$. Both have a finite support in the two consecutive intervals $[q_{i-1},q_{i}]\cup[q_{i},q_{i+1}]$, are piecewise cubic polynomials on each of them and have $C^{1}$ matching between them. They have the useful properties:

[TABLE]

By inserting (68) into eq. (27) and validating on an ensemble of $N_{x}\times N_{y}\equiv(2n_{x}+2)\times(2n_{y}+2)$ well chosen points $\{\bar{x}_{i}\bar{y}_{j}\}$ for each FC one obtains a linear system allowing to determine the unknown coefficients of the expansion (68):

[TABLE]

with $i,i^{\prime}=0,\ldots 2n_{x}+1$ and $j^{\prime}=0,\ldots 2n_{y}+1$ and $\alpha,\alpha^{\prime}=1,\ldots n_{a}$ where $n_{a}$ is the number of partial wave amplitudes of the Faddeev component $\phi_{a}$. Usually, by considering the fact that Faddeev components are regular at the origin, the $i=0$ and $j=0$ splines might be neglected as they will require coefficients with null values. Then the dimension of the linear system for N=3 is $d=3\times n_{a}\times(2n_{x}+1)(2n_{y}+1)$.

The matrix elements of the left hand side are generalized by

[TABLE]

and those of the right hand side by

[TABLE]

In the scattering problems, the boundary conditions introduce an inhomogeneous term in the right hand side of (69) and one is left with solving a linear system, generalized as

[TABLE]

where $x$ holds for the unknown coefficients $c_{\alpha ij}$. However the bound state problem is an homogeneous one

[TABLE]

where both $\lambda$ and $x$ are unknowns. Both problems can be unified by using the inverse iteration method. It consist in solving iteratively the inhomogeneous system

[TABLE]

starting with a trial value $\lambda_{0}$ and an initial guess $x^{0}$. One can show that under some conditions the series $x^{(0)},x^{(1)},x^{(2)},\ldots$ converges towards the eigenvector of $A$ whose eigenvalue is closest to the trial value $\lambda_{0}$. By doing so we are always reduced to solving an inhomogeneous linear system.

In view of solving very large linear systems

[TABLE]

the direct methods (Gauss elimination, LU decomposition,…) are not applicable for they require storage capabilities beyond the present technology and turns to be very slow. We use alternative methods, based also on iterative procedures which, starting from an ”educated guess” $x^{(0)}$, generate a series $x^{(i)}$ which minimize the residual $r^{(i)}=\mid Ax^{(i)}-b\mid$ until we consider it to be ”small enough”. They are based on the matrix-vector multiplication operations which require only to store a small amount of data, which allows reconstruction of the matrix A elements on the fly.

There are several families of iterative algorithms Saad_2003 . The so called Bi-Conjugate Gradient Stabilized (BICGSTAB) is a very robust one that we have extensively used.

One of the limitation of the iterative methods is however the number of iteration required until convergence. Even in the 3-body problem with $d\sim 10^{5}$ this number is prohibitive. To overcome this difficulty one uses the preconditioning technique, which consist in finding an approximation of the inverse matrix $A^{-1}$, say $\hat{A}^{-1}$, and solve instead of (72) the equivalent system

[TABLE]

We have systematically used this technique taking as a preconditioning matrix the one appearing in the left-hand-side of equations, the matrix $\hat{A}=L$. Its inversion can be performed exactly by means of the so called ”tensor trick”, introduced by Kok . It is based on the fact that the matrix $L$ is a sum of 4 terms having a tensorial structure $L_{i}=P_{i}\otimes Q_{i}$. We will illustrate this procedure in the case of one single amplitude

[TABLE]

One easily identifies the following tensorial structure

[TABLE]

where the factors

[TABLE]

have dimensions usually equivalent to the square root of the original matrix $L$. Let us rewrite $L$ it in the form

[TABLE]

and diagonalize in $\mathbf{C}$

[TABLE]

with U unitary and D’s diagonal matrices. We obtain

[TABLE]

which can be written in the form

[TABLE]

The inverse matrix is then easily obtained

[TABLE]

with

[TABLE]

This procedure, that can be generalized to N=4 Rimas_PhD_2003 and N=5 case, allows us to obtain a satisfactory preconditioning of the linear system (72) and a final solution with a reasonably small number of iterations ($\sim 20$).

The numerical protocol presented above involves as unique assumption the local expansion of the solution (68) in terms of cubic (eventually quintic) polynomials. Its accuracy relies only on the number of grid points ($n_{x},n_{y},\ldots$) for each coordinate $q=x,y,\ldots$ as well as in the number of partial wave amplitudes ($n_{a}$) on which the Faddeev-Yakubovsky components are expanded. The convergence of the results is studied by systematically increasing these numbers.

7 Summary

We have attempted to describe, in a unified scheme of increasing complexity, the Faddeev-Yakubovsky equations for the N=3, 4 and 5-body problems in configuration space, as well as the numerical methods allowing its solution.

They are based on a recursive splitting of the intrinsic N-body wave function, solution of the Schrödinger equation, in to as many components as there exist independent ways (or ”complete partition chains”) to break fully connected N-particle system into disconnected ones with a single interacting particle pair remaining. They are represented by the ”tree” diagrams of Figs. 1, 2 and 3. To each diagram is associated a function, named Faddeev-Yakubovsky component, which takes the values on $\mathbf{R}^{3(N-1)}$. The pioneering work of L.D. Faddeev for N=3 Fad_JETP39_1960 and its extension to an arbitrary N achieved by O.A. Yakubovsky Yakubovsky_SJNP5_1967 , both from the Leningrad/Saint Petersburg University, formulated the equations fulfilled by the ensemble of these components.

The original works of L.D. Faddeev and O.A. Yakubovsky were formulated in momentum space, in terms of T-matrix partitions, and solved under this form by Glöckle and collaborators for N=3 and N=4 Gloeckle_et_al The formulation in configuration space was due to S.P. Merkuriev and S. Yakovlev in a fruitful collaboration with the Grenoble theory group of the former Institut des Sciences Nucléaires Merkuriev_71 ; Merkuriev_74 ; GLM_PRL33_1974 ; MGL_AP99_1976 ; MY_1982 ; MY_1982_1983 ; MYG_NPA_1984 .

In configuration space, the Faddeev-Yakubovsky equations result into a system of partial derivative equations coupling the ensemble of the Faddeev-Yakubovsky components. After the partial wave expansion is performed, they turn into a coupled system of integro-differential equations on $\mathbf{R}^{N-1}$ with some smooth (N-2)-dimensional integral kernels which can be solved by standard linear algebra methods. Their generalization is natural, as it can be seen from (27), (49) and (67)

The Faddeev-Yakubovsky equations provide a mathematically rigorous approach for the full solution of the N-body problem. However its scalability with the number of interacting particles is a serious drawback that dramatically limits its range of applicability. The fact that the formal object to study is a function with arguments on ${\mathbf{R}}^{3(N-1)}$ is a first sign of it, but not the only one. The situation, discussed in some detail in Rimas_5N_FBS_2018 , is well illustrated in Table 1 where the dimension of the N-body solution is detailed in terms of number of equations, partial wave amplitudes, and the linear algebra problem.

We would like to notice that other rigorous schemes have been proposed for the solution of the N-body problem. One of them is the AGS equations AGS_1967 which have provided accurate results for the 3- and 4-nucleon problem Fonseca_Deltuva .

It is worth emphasizing also that such a rigorous mathematical schemes are not necessary when dealing with bound states or simple $1+(N-1)$ elastic scattering processes. The Schrödinger equation can then be directly solved by several methods, like GFMC GFMC , with Hyperspherical Harmonics Mario ; Pisa or Gaussian basis Emiko , NCSM NCSM ; dt_Nature_2019 , Lorentz Integral Transform LIT_94 which produces also very accurate results, in some cases well beyond the technical capabilities of the Faddev-Yakubovsky approach (see Orlandini_2017 for a more detailed review). However, the Faddeev-Yakubovsky partition of the wave function is interesting to increases the numerical convergence of the results or even unavoidable for an appropriate implementation of the boundary conditions Emiko ; Pisa_Scat .