Steps for the Solution of Faddeev Integral Equations in Configuration Space

Faddeev equations in configuration space for three-atom scattering processes that have previously been formulated in integral form allowing for additive and nonadditive forces, are now examined in terms of their numerical properties for a ”toy-model” case. The numerical implementation is based on a spectral decomposition in terms of Chebyshev polynomials. The potential for high accuracy of this method, of the order of 6 to 8 significant figures is one of the main motivations for the present investigation.The object of the equations are T−functions, that are the product of wave functions times potentials, that decay to zero in all directions. The driving and coupling terms are based on the two-body t−matrices, which describe the two-body correlations in each arrangement. The preferred form for both the driving term and the integral kernel terms is of a hybrid nature, the x−dependence is expanded into Chebyshev polynomials, and the y−dependence is given in terms of a Fourier series. The numerical properties of the driving term are examined in detail, and the integral kernel is left for a subsequent analysis.


Introduction
Our form of the equations have been modeled after the corresponding integral equations in momentum space [1], and a preliminary account is available [2], [3].A numerical evaluation of these equations has up to now not been implemented and progress in that direction for a simplified example will be described.An accuracy of 6 to 8 significant figures is envisaged by using a spectral expansion method in terms of Chebyshev polynomials [4], [5] for solving integral equations in configuration space.A basic ingredient is the integral of the two-body scattering t−matrix over a given function, which can now be obtained with a fast algorithm with an accuracy of more than eight significant figures [6].
The object of the equations are T −functions, T (x, y), that are the product of wave functions times potentials, that decay to zero in all directions.The numerical properties of the driving D(x, y) term of the coupled integral equations for the T -functions have now been established in the present investigation, and they serve as a guide for the choice of the algorithm for the solution of the overall coupled equations.The numerical evaluation of the latter is left for a future investigation.The most promising equation structure that emerges is a hybrid method, in which the behavior of T (x, y) in the x−variable is obtained by a spectral expansion into Chebyshev polynomials [4], while the behavior in the y−variable is expressed by a means of a Fourier integral representation.A Malfliet-Tjon potential is used for the two-body interaction, the three body potentials are set to zero, spin is ignored, all orbital angular momenta are set to zero, and the particles are taken to be identical.In arrangement 1 the x−vector is the displacement vector of particle 3 to particle 2, and the y−vector is the displacement vector from particle 3 to the center of mass of particles 2 and 3.
The basic equation to be solved in this toy model example is where D(x, y) is the driving term and K is the integration kernel, and where T (x, y) is the quantity to be calculated.The double integral will be denoted as K T in what follows.The directions of the vectors x and y have already been integrated out.We find that a convenient form for the driving term is where D F (q, x) can be expressed in terms of a function T F (q, x) whose x− dependence can be obtained very conveniently by means of expansions into Chebyshev polynomials.Here q 0 (q) is the momentum of the incident EPJ Web of Conferences (outgoing) particle on the target in the ground (excited ε t ) state of energy −|ε 0 |.The quantity T F (q, x) is obtained in terms of various nested integrals involving the two-body t−matrix at the energy ε t .The latter is the two-body energy of particles 2 and 3 obtained by subtracting from the total energy the kinetic energy of the particle 1 with momentum q.When q is near q 0 , ε t is close to −|ε 0 |, and T F has a pole.This, and other poles for excited bound states of the two-body system in the integral 2 can be handled by subtraction, as will be shown.
The purpose of the present report is to describe the procedure for obtaining numerical values for the function T F (q, x), and analyze it numerical properties so as to determine the most suitable ansatz for solving the integral equation (1), to be carried out in a future investigation.

The equations for D(x, y)
The driving term, given by Eq. (41) of Ref. [2] is composed of a direct term D c = t 1 P c Φ 1 and an exchange term D ac = t 1 P ac Φ 1 , where P c and P ac are the cyclic ant anti-cyclic permutation operators.The cyclic (or direct) term is where the vectors x ′′ and y ′′ are related to the vectors x ′ and y ′ according to For the anticyclic case t 1 P 13 Φ 1 one has Here Φ is the product of the ground state target wave function and the plane wave of the incident particle, where q 0 is the momentum of the particle 1 incident on the center of mass of particles 3 and 2. For a ℓ = 0 bound state of the target particle the wave function is given by where u 0 (x ′′ ) is the radial partial wave function.Further τ is the two-body T −matrix for particles 2 and 3, whose partial wave expansion is given by We take as zero the angular momenta in τ , and obtain where We now set to zero the following angular momenta e iq0•ρ0 → j 0 (q 0 |ρ 0 |).
Replacing the spherical Bessel function j 0 (z) by F 0 (z)/z, where F 0 (z) = sin(z) we obtain The integrations over the solid angle dΩ x ′ y ′ can be accomplished by making y ′ the polar axis and by using , dρ 0 increases positively as θ x ′ y ′ increases from 0 to π.Further using (ρ 2 ) 2 = − (4/3) (ρ 0 ) 2 + x ′ 2 + (4/3) y ′ 2 and defining the function ū0 (q 0 , x ′ , y ′ ) The integration over dΩ y ′ follows from ρ 3 dρ 3 = −y y ′ d(cos θ y,y ′ ), where ρ 3 = y − y ′ , with the result It is noteworthy that in the above result the dependence in y and y ′ is separable.After defining the function u S (q 0 , q, x ′ ) and 19 th International IUPAP Conference on Few-Body Problems in Physics one finally obtains The simple form of Eq. ( 18) is due to the separability of Eq. ( 15).As q changes the corresponding two-body energy ε q in τ 0 changes correspondingly according to where ε b is the binding energy of the two nucleons that form the target nucleus.In the appendix it is argued that the contribution from the anti-cyclic case to D is identical to that of the cyclic case, and hence 3 Numerical results for the driving term.
The two nested integrals ū0 (q 0 , x ′ , y ′ ), Eq. ( 14) and u S (q 0 , q, x ′ ), Eq. ( 16); and the solution of the L − S integral equation (17) for T F (q 0 , q, x), as described in Ref. [6], are carried out numerically using the spectral expansion techniques, first described in Ref. [4] and later used for several applications [7].A retrospective of this work is presented in Ref. [8].The bound state function u 0 (x), Eq. ( 8) is obtained with the iterative method described in Ref. [9].The Malfliet-Tjon potential, already multiplied by the factor 2µ/ 2 is of the form in units of f m −2 , as illustrated in Fig. (1).In order to suppress the singularity at the origin, the potential is replaced by a soft core at distances less than r cut = 0.2 f m given in terms of a second order polynomial in powers of r, matched to the potential in the vicinity of r cut as described previously.The resulting values of V at r = 0, 0.2, and 10 f m are approximately 231, 43.3 and −2.6 × 10 −7 f m −2 , respectively.The bound state occurs at an energy of −(0.1005443) 2 f m −2 .The bound state wave function u(r) has a value of 2.8 × 10 −9 at r = 140f m, but the calculation of the bound state energy is done out to a distance or 350 f m, using the procedure described in Ref. [9].The momentum of the incident particle, described by the y−coordinate, has the value of q 0 = 2 f m −1 .
The integration procedure consist in dividing the integration domain into partitions, expanding the integrand in each partition into a fixed set of Chebyshev polynomials (17 is standard), adaptively choosing the size of each partition such that the last three expansion coefficients have a magnitude less than the imposed accuracy parameter tol, and then performing the integrals over the Chebyshev polynomials using the standard matrices S L and S R [4].In the numerical applications [7] it was found that the accuracy of the final result was comparable to the value of tol, which in the present calculation has the value of 10 −8 .Allowing for a loss of accuracy of two significant figures, it is expected that the accuracy of T F (q 0 , q, x) is of the order of 1 : 10 −6 .Since the two-body τ −matrix has a pole at the bound state energy, which occurs when q = q 0 , special care has to be exercised in the vicinity of that point, as will be described in a forthcoming publication.
3.1 The x and q dependence of T F Numerical results for T F (x, q), Eq. ( 17) are illustrated in the figures below for several values of q, both smaller and larger than q 0 .The figures show considerable changes in T F as q varies, however, as q increases T F decreases rapidly.This is also illustrated in fig.(4), that shows the dependence of T F as a function of q for a fixed value of x = 1.0 f m The singularity at q = q 0 is clearly visible.The curve of T F (x, q) versus q for x = 0.5 f m is similar to that shown in Fig. (4), only its sign is opposite, and the magnitude is approximately a factor two larger.The rate of decrease of T F with increasing q can be seen from Fig. (5).For the expansion of T F (q, x) into Chebyshev polynomials it is advantageous to divide the domain of x into partitions, rather than expanding the whole domain all at once.This is illustrated in Fig. (6) for q = 1.The partitions are indicated in the figure.Already with 15 terms the size of the high order coefficients decreases to less than 10 −8 for each partition, which in turn is a measure of the truncation error of 05012-p.3Fig. 4. Dependence of TF on q for the fixed value of x = 1.0 f m.The subscript 1234 of TF indicates that this figure is a composite of four q−partitions, whose q−boundaries are: 1) from 0.01 to 1.00, 2) from 1.00 to 1.99, 3) from 2.01 to 2.4, and 4) from 2.4 to .4,all in units of f m −1 .The results for partitions 5 and 6, with 4 ≤ q ≤ 6, and 6 ≤ q ≤ 8, respectively, are displayed in Fig. (5).Fig. 6.The convergence of the Chebyshev expansion coefficients of TF (q = 1, x) for various sizes of the x−partitions.The global expansion in the whole interval 0 ≤ x ≤ 10 f m converges less rapidly than than the convergence in each of the five partitions into which the whole x−interval is divided, as described in the legend.

EPJ Web of Conferences
the expansion.That is not the case for the global expansion from x = 0 to x = 10.Expansions into other functions, such as Laguerre or Legendre polynomials or Fourier expansions converge much less slowly than the Chebyshev expansion.

The x and y dependence of the driving term
In order to obtain the y−dependence of D(x, y) it is necessary to carry out the integration sin(qy) T F (q, x) dq in Eq. ( 20).This is accomplished in terms of the analytic integrals over Chebyshev polynomials  20) over q.The result is obtained by applying Eq.( 21) to the four q−partitions 1 to 4 described in the caption to Fig. (4).The contribution from the pole singularity 1.99 ≤ q ≤ 2.01, shown in Fig. (8), is too small to be visible.
where the J is a Bessel function.So, for each fixed value of x one expands T F (q, x) into Chebyshev polynomials T n (q) and hence (23) The range of q is divided into partitions, and by means of Eqs.(21) the integrals over q are performed for each partition, and subsequently the results from each partition are added.The result for the fixed value of x = 1 is shown in Fig. (7) The following q−partitions were used: 0.01 → 1, 1 → 1.990, 2.01 → 2.4, 2.4 → 4. The singular partition 1.99 → 2.10 was left out in Fig. (7), and calculated separately.Its contribution is of the order of 10 −3 , which is too small to been seen in the figure.

Contribution from the pole term.
The contribution from the partition that contains the pole was calculated in two different ways.One, which uses a addition and subtraction method around the pole, similarly to what is done in the momentum representation, denoted as analytical method, and the (x, y), defined in Eq.( 24), with a numerical method based on an expansion of TF into Chebyshev polynomials, according to Eq.( 21).The good agreement between the two methods shows the power of Chebyshev expansions.
other, denoted as the numerical method, based on Eq. ( 21), is treated as described above.In that case the qdomains from 1.99 → 2.0 and from 2.0 → 2.01 are decomposed into 17 Chebyshev polynomials, each.Both methods gave very similar results, as shown in Fig. (8), which illustrates The good agreement shows that the "brute-force" numerical method, based on Eq.( 21) is robust, even in the vicinity of a singularity, but the more theoretical method, based on addition and subtraction of a known singularity, is preferable if high accuracy is desired.

The integral kernel K
The integral equation for the function T (x, y), according to Eq. (41) of Ref. [2] is The first term is the driving term and the second term contains the integral K = t 1 G 0 PT over the unknown function T , which will be the object of the present section.First the general structure of the nested integrals will be described, and the numerical realization is left to a separate publication.
Remembering that P contains a direct and an exchange term (cyclic and anti-cyclic permutations), making the transformation between different sets of Jacobi coordinates implied by P, and using the descriptions of the Green's function g 0 given by 05012-p.5

EPJ Web of Conferences
where F ℓ and G ℓ are the regular and irregular Riccati-Bessel functions for angular momentum number ℓ one obtains with ℓ = 0 the expression In the above dt = −d(cos θ x ′′ y ′′ ), the vectors ρ ′′ 2 , ρ ′′ 0 are defined in Section 2, and ρ ′′ 5 , ρ ′′ 4 are defined the same way for the anticyclic permutation.Similarly as was shown for the calculation of the driving term, the cyclic and anti cyclic parts are the same for the ℓ = 0 case, and hence the two integrals in the third line of Eq. ( 27) are equal.By transforming to the variable ρ ′′ 0 one finds dt = 2ρ ′′ 0 dρ ′′ 0 (3/4)x ′′ y ′′ ,and as a result the third line of Eq. ( 27) is The integration of T I (x ′′ , y ′′ ) over dx ′′ should include the Green's function g 0 (x ′ , x ′′ , ε q ) in Eq. ( 27) and is denoted as The integration of I 1 over dy ′′ should take into account the factor sin(qy ′′ ) which arises from the integration over the direction of q in the first line of Eq. ( 27) This quantity still has to be integrated over x ′ , with inclusion of the two-body τ 0 −matrix in Eq. ( 27) Combining all these expressions one finally obtains for This expression is of the same format as Eq. ( 20) for the driving term, rendering the overall equation for T This suggests that a convenient ansatz for T is in which case the two-variable equation for the unknown L would be The first term is the driving term, the second becomes an integral over T in view of the sequence of the nested integrals, Eqs. ( 28) through (33).If further the x−dependence of L(q, x) is expanded into Chebyshev polynomials where the x−domain is divided into partitions 1, 2, ..p max then, in view of Eq. (37), one obtains a set of coupled equations for the expansion coefficients a (p) n (q) that has to be solved numerically.Here p denotes the partition number, T (p) n (x) is the Chebyshev polynomial linearly transformed to the boundaries of each partition.The realization of such a scheme is left to a future study, however the structure of the equations to be solved (Eq.( 37)) is again of the hybrid form, in the variables x and q.

Summary and Conclusions
In this study a method of organizing the integrals needed to calculate both the driving term and the integral kernel of the Faddeev integral equation for the T −function was examined for a simplified "toy model".This study was motivated by the purpose of using the highly accurate spectral expansion method [4] for the numerical evaluations, all carried out in configuration space.It was found that both the driving term as well as the integral kernel are most suitably expressed in terms of a hybrid formulation in which the dependence on the x−variable of the T −function is best obtained in terms of an expansion into Chebyshev polynomials, together with a decomposition into partitions of the x−domain, while the y−dependence is best expressed in terms of a Fourier-type expression as an integral over a q−variable.The reason is based on the result that y−dependence of the driving term is highly oscillatory and decreases only slowly with distance, as illustrated in Fig. (7), while the Fourier coefficients for the driving term are not strongly oscillatory, and decrease reasonably rapidly with the q−variable, as illustrated in Fig. (4).The calculation of the driving term required three nested integrals, which were performed numerically, while the calculation of the integral kernel requires four nested integrals, to be carried out in a future study.The proposed set of of coupled equations to be solved for the Fourier expansion coefficients of the function T (x, y) were described near Eq.(37).The imaginary parts of the Greens functions were left out in the present exploratory study and the numerical evaluation of the integration kernel together with the solution of the coupled equations for the expansion coefficients will also be left for a subsequent study.
In conclusion, in order to demonstrate that our approach for solving the Faddeev integral equations in configuration space is feasible and can be rendered very accurate by means of spectral expansions into Chebyshev polynomials, the first step was to numerically evaluate and examine the driving term of the integral equation for the T −function for the "toymodel" case described here, and suggestions for how to formulate and solve the corresponding integral equation for the T −function have been presented.

Fig. 1 .
Fig. 1.The Malfliet-Tjon potential used throughout this text.It is already multiplied by the factor 2µ/ 2

Fig. 2 .
Fig.2.The dependence of TF on x for the values of q (in units of f m −1 ) indicated in the figure.The bound state pole occurs at q0.

Fig. 5 .
Fig. 5.The decrease of TF for large values of q, between 4 and 8 f m −1

Fig. 7 .
Fig.7.The function TF (x, y) for x = 1 obtained by integrating Eq.(20) over q.The result is obtained by applying Eq.(21) to the four q−partitions 1 to 4 described in the caption to Fig.(4).The contribution from the pole singularity 1.99 ≤ q ≤ 2.01, shown in Fig.(8), is too small to be visible.

Fig. 8 .
Fig. 8.Comparison of a rigorous numerical-analytical evaluation of T (pole) F