Few-body integrodifferential equation on Lagrange mesh

The two-variable integrodifferential equation for few-body systems is solved using the Lagrange-mesh method. The method transforms the equation into a system of algebraic equations that are solved as a non-symmetric matrix eigenvalue problem. Convergence properties of the solution to the integrodifferential equation in relation to the problem parameters is investigated. The accuracy of the converged solution is tested by calculating the binding energies and root-mean-square radii of selected few-body systems. The results are compared to those generated by other methods.


Introduction
The integrodifferential equation (IDE) approach is one of the several techniques devised to simplify the many-body Schrödinger equation for ease of numerical solution. The approach is based on the Faddeev formalism where the many-body Schrödinger wave function is decomposed into a sum of two-body amplitudes that depend on two global variables of the system [1]. The formalism also decomposes the Schrödinger equation into a set of coupled IDEs, for the amplitudes, that can be solved numerically with comparatively less complexity than the original Schrödinger equation. Such a decomposition explicitly accounts for two-body correlations in the system. All the IDEs are identical for systems of identical particles. Therefore, for such systems, only one IDE needs to be solved. The numerical solutions to the equations have been shown to be as accurate as other approaches to the solution of the many-body Schrödinger equation. However, unlike in the other approaches, the form of the IDEs in this approach does not depend on the size of the system.
The numerical solution of the IDEs has proved challenging because of their characteristic oscillatory non-local projection kernels. So far the two-variable integrodiferential equations have been solved using adiabatic approximation techniques [1,2,3,4] that decouples the equation into two single-variable equations. The challenge with the adiabatic approximation is that the energies of the systems are often extrapolated from the calculated ones [5]. This raises concerns, that have not yet been addressed, about the accuracy of the calculated amplitudes. Another, not so popular solution method for the IDEs, is the perturbation method [6,7], that solves the equations iteratively. In this method, the zeroth order (unperturbed) differential equations are separable in the two variables and can both be readily solved. Both the adiabatic approximation and the perturbation solutions of the IDEs become generally less reliable for few-body systems, and particularly for systems involving hard-core interactions [6]. In this paper the Lagrangemesh method [8] is explored for the direct solution of the few-body IDEs, for the first time. The Lagrange-mesh method is a variational method that employs basis functions that are constructed on a discretized problem domain. The basis functions so constructed have specialized properties, like orthonormality at the chosen discrete points, which leads to banded matrices for the Hamiltonian of the problem. This method was shown to generate converged solutions for the integral Schrödinger equation in momentum space [9], which involve non-local kernels similar to the few-body IDEs. The Lagrange-mesh method is widely discussed in the literature [8,9,10,11,12,13,14] and is here explored for solving the few-body IDEs.
This paper is organized as follows. In section 2 the Faddeev formalism is explained and the IDEs for few-body systems are outlined. Section 3 describes the construction of the Lagrange mesh matrix elements for the IDEs. Convergence tests of the numerical solutions are discussed in Section 4 and conclusions are reported in Section 5.

The few-body integrodifferential equation
Consider a system of A particles of equal mass m and position vector r i , i = 1, · · · , A, in which two-body interactions V ( r ij ), where r ij = r i − r j , are dominant. The Schröndiger equation for the many-body system has the form ⎡ where x ≡ ( r 1 , r 2 , · · · , r A ), E the energy and Ψ( x) the total wave function of the system. In the Faddeev formalism, the solution to this many-body problem is constructed by decomposing the many-body Schrödinger wave function in the form [15] Ψ where ψ ij ( x) are Faddeev two-body amplitudes. These amplitudes depend on all the degreesof-freedom of all the particles in the system and satisfy 1 2 A (A − 1) coupled equations of the form whereT represents the kinetic energy operator in Eq.
where F (r, z) are the reduced Faddeev two-body amplitudes and is the relative orbital angular momenta of the interacting pair. In Eq. (4), the kinetic energy for the center-of-mass of the is used. The weight function w 0 (z) = (1 − z) α 0 (1 + z) β 0 is defined by the parameters α 0 = (3A − 8)/2 and β 0 = + 1/2, which are specified by the problem of interest. The kernel f (z, z ) projects the angular variables for the other pairs of particles into that of the focus pair. This kernel is defined in terms of the Jacobi polynomials P α,β are related to pairs of particles that are connected and those that are not connected to the {ij} pair. For systems of particles with differing masses, the first term in the square brackets is slightly modified. Note that (4) represent a set of coupled equations for the amplitudes with different values of . This set of equations is reduced to one IDE by setting = 0. The resulting equation is then modified by adding the zeroth (K = 0) and subtracting the first (K = 1) order hypercentral potential multipoles, from the equation. The zeroth order multipole accounts for the addition of the effects of higher orbital angular momentum [16], while the first order multipole accounts for the removal of the effects of the spurious amplitudes [5]. Only the effects of spurious amplitudes are considered in the present work. As a result, the two-variable s−state intgrodifferential equation solved is [15] and the subscripts on the amplitude are now omitted since they are redundant. In anticipation of later discussions, equation (9) is cast in the eigenproblem form whereĤ For bound-state solutions to (11), which are the focus of this paper, specific boundary conditions need to be satisfied by the amplitudes F (r, z quantum mechanics involving the requirement for a realistic description of the physical properties of the system. The form of the interaction potentials in the equation determines the type of boundary conditions to be satisfied by the amplitudes. The asymptotic short-and long-range radial forms of the amplitudes, for bound states, are En r (13) where E n is the energy of one of the states of the system. These asymptotic forms translate into the boundary conditions which are consistent with suppressing the singularity of (12). The amplitudes are also required to vanish at all other singularities of the Hamiltonian. Since the angular-related part of the kinetic energy operator does not appear to have singularities, the boundary conditions are adopted. That is, if the potential does not have singularities in the z domain, then F (r, ±1) are just functions of r. In the next section the numerical solution of (11) is constructed using the Lagrange-mesh method.

Lagrange mesh matrix elements
The Lagrange-mesh method is widely discussed in the literature and the reader is referred to [8] for details on this method. In this method, the problem domain [a, b] is segmented into intervals defined by a by a quadrature where g(x) is a real function, while x i and λ i are the quadrature abscissas and weights, respectively. The method is characterized by basis functions f j (x) that have the properties [8] where V (x) is a real function and f (k) ≡ d k f/dx k . These properties reduce dynamical equations for quantum mechanical problems to a simple set of algebraic equations that are readily solved. A simple example application of the Lagrange functions is provided by the one-dimensional Schrödinger equation where V (x) is a scalar potential, ψ(x) the wave function, and E the energy of a given system. The Lagrange-mesh solution to this equation is constructed by approximating the wave function as where c j are variational parameters and f j (x) the N Lagrange functions defined on a grid a ≤ x j ≤ b (j=1,2,3,. . . ,N). Taking the matrix elements of (20) with the expansion (21) leads to This equation simplifies, thanks to the properties of the Lagrange functions, to the set of algebraic equations [8] for the variational parameters c i . Note that the matrix elements in (23) require only the evaluation of the known functions f j (x) and V (x) at the known points x i . The Lagrangemesh method is similarly applied to (11), however, this equation presents non-Hermitian and non-local operators, which lead to a non-symmetric matrix eigenvalue problem. The amplitudes in (11) are approximated with the series where C ij are variational parameters, R i (r) and U j (z) bases functions that satisfy the boundary conditions (14) and (15) where h σ K is the normalization coefficient of and w L (r) = r σ e −r the weight function for the Laguerre polynomials L σ K (r) of order K. The parameter σ is usually chosen so that (25) displays the required boundary behaviour of the exact wave function at small r. However, here σ = 0. The regularized Lagrange-Jacobi functions have the form [17] where μ is the regularizing parameter and γ = α + β + 1. The values of the parameters ( α, β ) are not necessarily the same as those specified in the previous section. However, when this is the case, then the corresponding polynomials P α,β N (z) constitute the eigenfunctions of the angular component of the kinetic energy operator. In general, these parameters need to be optimized for a given potential.
Using the expansion (24) in (11) with (25) and (26), and constructing the matrix elements for the equation, leads to the matrix eigenvalue problem for the variational parameters, where λ i = w J i /w J (z i ) (w J i are the quadrature weights), h a scaling parameter and The matrix elements T r ij are those given in [8] while T z ij are taken from [18]. It should be noted that the matrix elements in the eigenvalue problem (27) are functions of the quadrature weights and abscissas, as well as the free parameters h, α, β and μ. Convergence properties of the solutions to (27) with variations in the bases sizes are discussed in the next section.

Convergence tests
The evaluation of the matrix elements discussed in the previous section require the specification of the parameters {h, μ, α, β}, as well as the roots of the Laguerre and Jacobi polynomials. The mesh was based on Gaussian quadrature, which were determined using standard routines [19] with a tolerance of 3 × 10 −14 . The regularisation parameter μ is determined by the boundary conditions (15). The value μ = − 1 2 was adopted. The scaling parameter h is determined by the hyperradial range of the total potential of the system [20]. Each particle in the system experiences an average hyperradial potential described by the hypercentral potential V 0 (r). Figure 1 show the three-body to six-body hypercentral potentials described by (8). As can be seen in the figure, the range of this hypercentral potential depends on the size of the system. In this work, the scaling parameter was chosen such that h = r m /x N where r m is the range of the hypercentral potential and x N the last zero of the Laguerre polynomial. The range r m = 20 fm was used for the three-body and four-body, while r m = 30 fm was used for the five-body and six-body systems. The parameters of the functions U j (z) were set to (α, β) = (α 0 , β 0 ), which are related to the eigenfunctions of the angular kinetic energy operator. Test applications of the matrix elements were conducted with the central Volkov potential [21] and the Ali-Bodmer potential [22] for three-, four-, five-and six-body systems. The variation of the ground-state energies of the systems with the bases sizes were investigated. The investigations were restricted only to the cases of equal radial and angular bases size N = K. In addition, a restriction of M ≥ N , on the number of terms M in the series in (6), was imposed, based on tests results of the projection function on polynomial and exponential functions.

Three-body systems
In the first case, the three-body s-state IDE was solved for the 3n ( 2 /m = 41.47 MeV · fm −2 ) and 3α ( 2 /m = 10.366 MeV · fm −2 ) systems. MeV for the 3n system and to −4.141 MeV for the 3α system. These results coincide with those reported in the literature [23,24]. The IDE results, on the other hand, converge very slowly to −8.111 MeV and −3.302 MeV for the Volkov and Ali-Bodmer potentials, respectively. These results are about 5% and 20% less than those of the Faddeev for the two potentials. There are several possible causes of these discrepancies, one of which is the effects of the spurious potential [25].

Many-body systems
To investigate the effects of the spurious potential on the results, ground state energies for A-body systems, where A = 3, 4, 5, 6, were calculated with the spurious potential eliminated. The E IDE results are shown in table 2 and table 3 for the Volkov and Ali-Bodmer potential, respectively. The energies for all the systems converge quite slowly, and, as indicated in [25], do not appear to be significantly affected by the spurious potential. There is a general discrepancy of about 11% in the Volkov results compared to the s-wave results of the hyperspherical harmonics method [23] for the same potential. However, the hyperspherical harmonics results are obtained with bases sizes in the orders of millions, for six-body system, for example. The quality of the wave functions is often tested by calculating root-mean-square radii of the systems studied. The results for the root-mean-square radii of three-body and four-body systems are presented in table 4 for the Volkov and Ali-Bodmer potentials. The results for the 3α system are about 12% higher than those reported in the literature [24].

Conclusion
The Lagrange-mesh method was explored for solving the few-body IDEs. Matrix elements for the non-Hermitian and non-local operators of the equations were constructed using regularized Lagrange-Laguerre and Lagrange-Jacobi functions. The resulting non-symmetric matrix was first solved with standard nuclear potentials for three-body systems. To investigate the effect of the non-symmetric form of the eigenproblem on the convergence of the numerical solution,   Table 3. Same as in Table 2 for the Ali-Bodmer potential. the non-symmetric form for the IDE is not responsible for the slow convergence of the solution. Therefore, the convergence problems could emanate from the form of the basis functions and/or the implementation of boundary conditions of the problem. The IDE was also solved for the ground-state energies of A-body systems where A = 3, 4, 5, 6, with the contribution of spurious potential removed. Ground state energies and root-meansquare radii of the A = 3, 4 systems were calculated. The results showed very slow convergence with the basis sizes, to values ∼ 11% less than those obtained using the hyperspherical harmonics method. Both the ground-state energies and root-mean-square radii of the systems were poorly determined, reflecting the possible poor quality of the numerical wave functions. The convergence of the IDE Lagrange-mesh solutions, albeit slow, and the rapid convergence of the Faddeev Lagrange-mesh solutions, indicate that the method is a promising tool to generate direct solutions to the IDE. Further investigations on the effectiveness of the Lagrange-mesh solution of the IDE are necessary. Work is underway to improve the numerical treatment of the few-body IDE.