Hyperspherical harmonics expansion on Lagrange meshes for bosonic systems in one dimension

A one-dimensional system of bosons interacting with contact and single-Gaussian forces is studied with an expansion in hyperspherical harmonics. The hyperradial potentials are calculated using the link between the hyperspherical harmonics and the single-particle harmonic-oscillator basis while the coupled hyperradial equations are solved with the Lagrange-mesh method. Extensions of this method are proposed to achieve good convergence with small numbers of mesh points for any truncation of hypermomentum. The convergence with hypermomentum strongly depends on the range of the two-body forces: it is very good for large ranges but deteriorates as the range decreases, being the worst for the contact interaction. In all cases, the lowest-order energy is within 4.5$\%$ of the exact solution and shows the correct cubic asymptotic behaviour at large boson numbers. Details of the convergence studies are presented for 3, 5, 20 and 100 bosons. A special treatment for three bosons was found to be necessary. For single-Gaussian interactions, the convergence rate improves with increasing boson number, similar to what happens in the case of three-dimensional systems of bosons.


Introduction
The expansion in hyperspherical harmonics (HH) of a many-body wave function is a powerful scheme that allows the Schrödinger equation to be solved in a systematic way by reducing it to a system of coupled equations of only one dynamic variable for any number of particles N [1]. This method is a standard one to describe few-body physics in three dimensions. The hyperspherical expansion in one dimension for a three-body problem has been considered in Refs. [2,3] for contact two-body potentials. Reference [2] also comments on the validity of the lowest-order approximation of the HH expansion at large N . The authors point out that the N -dependence of the lowest-order energy is quadratic while the exact solution gives a cubic N -dependence. This seems to contradict more recent observations made for three-dimensional systems that the HH expansion for a single-Gaussian two-body potential converges better with larger N [4,5]. This means that the lowest-order energy should get closer to the exact solution with increasing N , which was checked numerically in Ref. [5]. Since the contact interaction can be seen as a limit of a single-Gaussian potential with a very short range, the same argument should hold for the contact interaction as well. Therefore, the statement of Ref. [2] regarding the quadratic versus cubic dependence of the lowest-order HH and exact energies needs to be reconsidered.
In this paper, the HH expansion in one dimension is revisited and formulated in terms of an expansion on a single-particle oscillator basis [6,7]. The HHs are expressed as functions of symmetrized oscillator states. The matrix elements of the potential between these HHs are obtained from integral transforms for N > 3. The resulting coupled hyperradial equations depend on the hyperspherical quantum number K and must be truncated at some value K max . For a given truncation K max , they are solved with the Lagrange-mesh method, an approximate variational technique which has the simplicity of a mesh calculation [8][9][10]. For a number of applications, this method provides accurate results with small numbers of basis functions [10]. Its accuracy is equivalent to the accuracy of a variational calculation with the same basis [11]. This property is in particular verified for HH equations in Refs. [10,12].
The convergence of the HH method with respect to K is explored for the contact interaction, for which the exact solution is known, and for different ranges of a single-Gaussian interaction. The validity of extrapolation techniques is analysed.
In Sect. 2, the method of HH expansion is described for the one-dimensional case. The Lagrange-mesh method and an extension of this method are applied to the resulting hyperradial equations in Sect. 3. In Sect. 4, the special case of the lowest-order energy is considered and the convergence in the case of a contact interaction is analysed by comparison with the exact results. The same study is performed for Gaussian interactions in Sect. 5. Sect. 6 contains a summary while the Appendix provides information about some mathematical details relevant to the HH method.

Choice of Coordinates
An N -body system is described by a set of N − 1 normalised Jacobi coordinates ξ 1 , ξ 2 , . . . , ξ N −1 defined in one dimension, where r i is the individual coordinate of the i-th body. These coordinates form the (N − 1)-dimensional vector ρ. Its length is given by the hyperradius ρ, where R = ( N i=1 r i )/ √ N is the normalised coordinate of the centre of mass. The other hyperspherical coordinates, the hyperanglesρ ≡ {θ 1 , θ 2 , . . . , θ N −2 }, can be chosen in many different ways but they are not used in the present approach.

Schrödinger Equation in Hyperspherical Harmonics Basis
The kinetic-energy operatorT in hyperspherical coordinates is split into hyperradial and hyperangular parts, where n = N − 1 is the dimension of the space formed by the Jacobi coordinates. The eigenfunctions Y K γ (ρ) of the hyperangular part ρ , where K is the hypermomentum, form a complete orthonormal set. They can be highly degenerate and different harmonics belonging to the same K are labelled by γ , which represents all the angular momentum quantum numbers. Below, they are constructed symmetrical with respect to permutations of particles in order to be able to describe bosonic systems. The wave function of the N -body system in these coordinates, (ρ), can be expanded over complete sets of hyperspherical harmonics, The hyperradial parts χ K γ (ρ) satisfy the infinite set of coupled differential equations where is the generalized angular momentum, m is the mass of the particle and the hyperradial potentials V K γ,K γ (ρ) are the matrix elements of the interaction potentialV between HHs,

Construction of the Hyperspherical Harmonics
Following Refs. [5][6][7], we construct the hyperspherical harmonics Y K γ (ρ) from linear combinations of eigenfunctions of the harmonic-oscillator Hamiltonian, where ω is the oscillator frequency. This Hamiltonian can be rewritten in hyperspherical ρ,ρ and centre-ofmass R coordinates as The eigenfunctions of the Hamiltonian (10) corresponding to the lowest energy of the centre-of-mass motion and the internal energy eigenvalues (2κ + K + n/2)hω are factorised as where c.m.
The number of positive-parity hyperspherical harmonics for some values of K is shown in Table 1 for several many-boson systems. For N ≥ 4, the number of positive-parity HHs is non-zero for any K except K = 2. For N = 3, only those HHs where K is a multiple of six exist, one HH per K . This does not contradict the HH basis from [3], where the HHs are just the functions e i K θ of the angle θ = arctan(ξ 2 /ξ 1 ) with K values of 0, ±6, ±12, ±18, . . . Indeed, transforming the e i K θ basis into a cos K θ and sin K θ basis results in non-zero potential matrix elements only for cos K θ , leaving effectively only one HH channel per K , i.e.,

Hyperradial Potentials
To calculate the hyperradial potentials V K γ,K γ (ρ) that enter Eq. (6), we use the technique developed in Ref. [6] to obtain in which the integration path bypasses the origin in the counterclockwise direction. The integral in Eq. (18) is the inverse Laplace transform of the function s −(K +K +n)/2 V K γ,K γ (s), many of which are tabulated in Ref. [13]. The potentials contain the usual matrix elements K μ |V | K μ in the oscillator basis, which are calculated using standard many-body techniques for b = s −1/2 . The integration over s can be easily performed analytically before making the integration over {r i }. For some potentials, among which the contact and Gaussian ones, it is possible to make the integration over {r i } first and then integrate over s analytically using available tables of inverse Laplace transforms (see "Appendix C").

Solving the Hyperradial Equations
To solve the system (6) of coupled hyperradial equations, we use the Lagrange-mesh method [8,10]. This method is an approximate variational calculation making use of M Lagrange basis functions related to the M mesh points of an associated Gauss quadrature. As the matrix elements are evaluated with this quadrature, the resulting equations take the form of a calculation on a mesh. In spite of this approximation, numerical tests on a variety of problems show that a Lagrange-mesh calculation is essentially as accurate as the corresponding variational calculation involving numerically exact matrix elements [11] (see also Sect. 4.2). This striking property is obtained provided that no singularities appear in the matrix elements and that the behaviours of the wave function at the origin and at infinity are well simulated by the basis. The problem raised by the existence of a singularity can be solved with a regularization technique [9].

Lagrange-mesh Method
For a variable varying from zero to infinity, it is convenient to use the M regularized Lagrange-Laguerre functionsf where L α M (x) is a generalized Laguerre polynomial [14] depending on parameter α > −1 and h α M = (M + α + 1)/M! is the square of its norm. The associated zeros x j , j = 1, . . . , M, are given by They correspond to the mesh points of the Gauss-Laguerre quadrature where the λ k are the Gauss weights. The x k and λ k depend on α and M. This quadrature method is exact when g(x) is the product of a polynomial of degree at most 2M − 1 by the weight function of the Laguerre polynomials The Lagrange functions (20) are thus the product of the square root of the weight function w(x) by x and by a polynomial of degree M − 1. They vanish at the origin and verify the Lagrange propertŷ Without the factor x, the Lagrange functions would be exactly orthonormal. The usefulness of the regularizing factor x is explained below. The functions (20) are not orthogonal but are still orthonormal at the Gauss quadrature (22). As explained in Ref. [10], this basis will be treated as orthonormal in the following without loss of accuracy.
The hyperradial functions are expanded as where h is a scaling parameter and h −1/2f j (ρ/ h) is normed at the Gauss approximation. Their behaviour at the origin is given by For x → 0, the Lagrange functions behave asf We choose thanks to the Lagrange property (24) and are thus treated as orthonormal. The crucial simplification of the Lagrange-mesh method is that the matrix elements of the potential are approximated as The potential matrix is diagonal. Thanks to the regularization by x, the matrix elements of 1/x and 1/x 2 are exact with the Gauss quadrature. The hyperradial equations take the form of mesh equations, The kinetic-energy matrix elements of −d 2 /dx 2 read at the Gauss approximation [10,15] and They are simple expressions involving the Laguerre zeros x j . Although approximations (29) and (30) are rather poor for individual matrix elements, the resulting Lagrange-mesh calculation can be very accurate [10,11]. This striking property is illustrated here by the results in Sects. 4 and 5.

Hybrid Lagrange-mesh method
For large numbers of bosons, the increasing degeneracies of large K values (see Table 1) lead to strongly increasing sizes of the matrices. Computer times may then forbid reaching a sufficient convergence with respect to K . To partly compensate this effect, it would be efficient to use decreasing numbers of mesh points when K increases. If there is no significant loss of accuracy, this would reduce the size of the matrix and thus the computer time. However, the simplicity of (30) is then lost in the calculation of some matrix elements connecting different K values since the associated meshes and Gauss quadratures then differ. These matrix elements can nevertheless be accurately calculated in many cases with a third Gauss quadrature [16]. This leads to a hybrid calculation where some matrix elements are calculated like in the previous subsection and a different method of evaluation is used for the other ones.
Such a method opens the way to different strategies. One can keep a unique value for α as before and vary the numbers of mesh points and scale parameters. One can also vary the α values. For example, one can use basis functions with the correct behaviour (26) at the origin for each partial wave by choosing These hybrid calculations may reduce the matrix size at the cost of a more complicated search for optimal numbers of mesh points and scale parameters. As shown below, the best strategy we have found happens to be quite simple.
With respect to the Lagrange-mesh method, a new type of calculation is necessary. For α = α , one must calculate the matrix elements between Lagrange functions, In the case V (ρ) ∝ 1/ρ which will be encountered below, the integral (35) can be exactly computed with a Gauss-Laguerre quadrature with Laguerre parameter and scale parameter provided the number M of Gauss mesh points x k and weights λ k verifies The matrix element of 1/ρ is then exactly given by The choice α = (α + α )/2 + 1 is also possible. A similar formula should be accurate for more general hyperradial potentials. If these potentials are regular at the origin, the value of α may be increased by one unit. The value of M must be large enough to reach a sufficient accuracy.

Special Treatment of N = 3
For N = 3, one has L 0 = −1/2 and the K = 0 'centrifugal' term becomes attractive. This gives rise to difficulties in most numerical approaches. The previous methods are not valid since α = −1 appears. A special treatment of K = 0 is necessary. It is made possible by introducing non-regularized Lagrange functions for K = 0 in the hybrid method. As shown below, it leads to a fast and accurate computation.
The behaviour of the hyperradial functions at the origin is given by The function χ 0 thus behaves as ρ 1/2 . The M non-regularized Lagrange functions [8,10] are defined as With they have the required behaviour at the origin. With this basis, the matrix elements of 1/x 2 and of −d 2 /dx 2 do not exist but the matrix elements of the full kinetic operator converge [10,17] and are given for i = j by and for i = j by These expressions replace the parenthesis in the first term of (31) for K = 0. The Gauss quadrature is not accurate for matrix elements of 1/x but one can use the exact expression [10,17], This expression is useful in the case of contact interactions described below. In this case, the potential terms calculated with the Gauss quadrature in (31) must be replaced by exact expressions with the help of (45).
Other K values are treated like in the hybrid method. The expression of non-diagonal matrix elements between K = 0 and K = 0 must be adapted from (39) by replacingf α i by f i for K = 0. Finally, when V (ρ) ∝ 1/ρ, an exact variational calculation is also possible. All matrix elements for K = 0 and all matrix elements of the potential are exact. The overlaps of the regularized functionsf j (x) are given by Eq. (3.71) of Ref. [10] and the exact kinetic-energy matrix elements are given by Eq. (3.77). The variational results are then obtained from a generalized eigenvalue problem.

HH Expansion for Zero-Range Forces
For contact (zero-range) two-body forces (V 0 > 0), the exact energy is known analytically [18]. This allows precise tests of the convergence of the HH expansion. The Schrödinger equation verifies the scaling property The same property is verified by the system of coupled equations (6) and by its approximation (31), even when the system is truncated at any value of K . The matrix elements of the potential K μ |V | K μ in the oscillator basis are proportional to √ s, which gives with (19) the hyperradial potentials as The system (6) looks like a coupled Coulomb problem for a two-body system in three dimensions. The coefficients for the three-boson case in one dimension can be derived from Eq. (3.3) of Ref. [3]. In the basis (17) constructed in the present paper, the non-vanishing coefficients are where K and K are multiples of 6. For all other numbers N of bosons, the c K γ,K γ are obtained from numerical calculations.

K = 0 Limit
For K = 0, (19) reads with (16), With an inverse Laplace transform [14], the corresponding matrix element (18) reads If we keep in expansion (5) the K = 0 term only, then the hyperradial equation (6) is exactly the same as the one for a hydrogen-like atom with non-zero angular momentum L 0 = 1 2 (N − 4), with The solutions of this hydrogen-like equation are well known. Its lowest-energy eigenvalue for each L 0 value is and the corresponding wave function is where the inverse 'Bohr radius' reads In the N → ∞ limit, one has So, both E K =0 and E exact have the same cubic asymptotic behaviour. The N 2 behaviour of E K =0 in Ref. [2] comes from a missing factor ( where i is the solid angle in an (i + 1)-dimensional space, in their analytical expression for the K = 0 energy. Indeed, for N large, one has The evaluation of the solid angle d N −1 is explained in "Appendix D". Interestingly, the same N 3 behaviour can be obtained by minimizing the expectation value of the two-body interaction in the lowest-order oscillator basis. It gives from which it follows that E min HO is always higher than E K =0 but approaches the latter asymptotically. Both E min HO and E K =0 are shown in Table 2 in comparison with the exact solution E exact .

Convergence of HH Expansion for Contact Interaction in a Three-Boson Case
The convergence of the HH expansion for a three-body system in one dimension with the contact interaction was shown to be slow [3]. It was suggested there that this convergence has an exponential behaviour. We confirm the slow convergence but we disagree with its exponential form. Because of the small size of the matrices involved in the 3-body case and availability of exact expressions for the hyperradial potentials, it is possible to reach very high values of K max , such as K max = 6000, so that convergence can be studied precisely. Because of the scaling property (47), the convergence of the HH expansion for a contact interaction is independent of the choice of V 0 for any fixed N . Therefore, we choose the parameters As shown for hydrogen-like atoms in Ref. [10], a Lagrange-mesh or variational calculation with a single Lagrange function can give the exact lowest energy for K = 0 since the corresponding wave function can be simulated by this Lagrange function (compare (55) with (20) or (41) for M = 1) and the Gauss quadrature is exact. This is obtained with the scale parameter This should thus be a good scale parameter for the full calculation.
In calculating E(K max ), we had to use a modification of the Lagrange-mesh approach because its original version is not valid due to the attractive 'centrifugal' term at K = 0. This case is treated as explained in Sect The convergence with respect to the number of mesh points is fast as shown in Table 3. An accuracy of 10 −10 is reached with M = M = 4. The validity of this hybrid mesh calculation can be tested. Indeed, one can easily perform an exact variational calculation by exactly evaluating the remaining matrix elements, i.e., the overlaps and the kinetic-energy matrix elements of regularized functions. This leads to a generalized eigenvalue problem. The comparison in Table 3 shows that the hybrid Lagrange-mesh results converge faster than the variational results but that the difference is small and decreases with increasing M = M . The slower convergence of the variational method is probably be due to the non-orthogonality of the basis.
We have found that the behaviour of E(K max ) is well described by where E(∞) ≡ E exact , with a 0 = 0.183786 and a 1 = 3.45912. This is illustrated in Fig. 1, which shows that the ratio 1/(E(K max ) − E(∞)) is practically a linear function over the huge interval 30 ≤ K max ≤ 6000. This ratio should have increased exponentially if the convergence was exponential.

Convergence of HH Expansion for Contact Interaction for N > 3
Here we consider the typical values N = 5, 20 and 100 for which we discuss the convergence with respect to The convergence with respect to various strategies is analysed in Table 4 at K max = 20 for the three cases. With only two mesh points per partial wave, the relative error is of the order of 10    Table 1. In the following, we shall truncate the case with N = 5 bosons at K max = 40. The sum of the degeneracies over K is 145 and the matrix size is thus M × 145. The matrix sizes become much larger when the number of bosons increases. For N = 20 at K max = 32, the matrix size is M × 1451. For N = 100 at the same K max , the matrix size is M × 1507.
A hybrid treatment reducing these sizes would thus be helpful. We have compared various strategies and found that the choice (34) leads to complicated searches for optimal sets of h values, without significant improvement. The simplest hybrid strategy is to reduce the number of mesh points beyond some K value. The reduction from M = 4 to M = 3 beyond K = 10 is tested as H1 in Table 4 for K max = 20. There is a little but not significant loss of accuracy. It can be partly compensated by a slight increase of h as shown in row H2. For high K max values, this can reduce the diagonalisation time in a useful way. For N = 20, the size is reduced from 5804 to 4365 with a negligible loss of accuracy. For N = 100, the matrix size is reduced from 6028 to 4533.
The ratio E(K max )/E(∞) is shown in Fig. 1 as a function of the model-space K max for N = 5, 20 and 100. The K max value is 40 for N = 5 and 32 for N = 20 and 100. Figure 1 shows a very slow convergence for all selected N . However, for larger N , the energies for a fixed E(K max ) are closer to E(∞). We have also plotted the 1/(E(K max ) − E(∞)) values. Similar to the three-body case, their behaviour is very close to linear at large K max . However, a close look shows small-scale irregularities that prevent to unambiguously determine the parameters of such linear functions. If E(K max ) is indeed governed by Eq. (60) then its convergence within five significant digits would be achievable in a huge model space with K max as large as several hundreds. Working with such model spaces is not realistic for N ≥ 4.
We tried to parametrize E(K max ) at K max ≥ 20 by other functional forms, exponential and power functions. Exponential extrapolations result in asymptotic values of E which are slightly higher (of the order of 0.2%) than the exact solutions. The power extrapolations give better asymptotic values of the energies. However, Table 5 Converged and extrapolated energies E(∞) (in K) for an N -boson system with a Gaussian interaction of range a (in a.u.) and depth V g = V 0 /( √ πa) that gives the same volume integral as a zero-range force of strength V 0 N a = 0 a = 0.05 a = 0.1 a = 0.2 a = 0. The calculated values are shown for V 0 = 10 K. The underlined digits are uncertain by 1-2 units their accurate determination is problematic within the model space considered because of the small-scale irregularities.

HH expansion for Gaussian Two-Body Interaction
For a Gaussian two-body interaction with range a and depth V g , the convergence of the HH expansion is determined by V g a 2 , because of the scaling property r i → ar i , The energies E(K max ) converge faster with larger V g a 2 and larger N , which is exemplified below in Fig. 2. To compare this convergence with that for the contact interaction, V g has been chosen as V 0 /( √ πa) so that both the contact and Gaussian interactions have the same volume integral. For N = 3, the non-vanishing hyperradial potentials are obtained by a direct calculation of (8) using basis (17) as where K and K are multiples of 6 and I n is a modified Bessel function of the first kind [14]. The K = 0 kinetic energy is calculated with (43) and (44). For a → 0, one recovers (48) and (49). As for the contact case (a = 0), the Lagrange-mesh approximation (30) is not valid for a < 1. All potential matrix elements involving K = 0 are thus computed numerically with (39) and M = 40 − 60. The energies obtained for several two-body interactions that have the same volume integral, but different ranges a are presented in Table 5. The ranges a are given in a.u. while the depth V 0 = 10 is given in K and h 2 /m = 43.281307 is chosen to be the same as in Ref. [4]. The energy for a = 0 is exact. With K max = 120, convergence is reached for the displayed digits for the other a values. These results are obtained with h = 1.6 given by (59) and M = 20 or 30. The binding energy progressively decreases when a increases. The analysis of the E(K max ) behaviour suggests that it is determined by the E(∞) + a 1 /(a 2 + K max ) a 3 law. The coefficients a 3 are very sensitive to the number of energies E(K max ) used to make an extrapolation but, generally, they increase with the range a.
For N ≥ 4, the Lagrange-mesh method is used. The potential matrix elements are calculated as explained in "Appendix C". For the larger a values (a ≥ 0.5 a.u.), converged energies are shown. For a < 0.5 a.u., extrapolated values are given. The underlined digits in Table 5  One can see that, with increasing N , the energy of the N -boson system deviates more strongly from that obtained with the contact interaction as the interaction range gets larger. At a = 1, K max = 6 becomes a good approximation with an accuracy better than 10 −4 for N ≥ 20. For a > 3, the accuracy of the K = 0 approximation is better than 10 −3 .
The E(K max ) values show more irregularities with increasing N so that it was not possible to determine an exact functional form of the E(K max ) behaviour. The inverse power law, E(∞) + a 1 /(K max ) a 2 , however,  Table 2 for N = 3, 5, 20 and 100. For comparison, calculations with the corresponding zero-range interaction are also presented seems to be appropriate to describe E(K max ) at 20 ≤ K max ≤ 40, and it was used to determine the E(∞) values shown in Table 5. To illustrate the E(K max ) behaviour, the ratios of energies E(K max ) to the extrapolated values E from Table 5 are displayed as a function of K max in Fig. 2. The accuracy of these values depends on V g a 2 but in all cases their uncertainties are too small to be seen.

Summary
The hyperspherical-harmonics expansion for bosonic systems in one dimension is revisited and formulated in terms of an expansion on a single-particle oscillator basis. The hyperspherical harmonics are expressed via linear combinations of symmetrized products of single-particle states that contain only the lowest-energy centre-of-mass motion and do not contain hyperradial excitations. They are constructed by diagonalising the matrices of the square of the centre-of-mass coordinate and the hyperangular part of the multidimensional Laplacian. Then the hyperradial potentials are calculated with known methods of many-body physics that use an oscillator basis. The HH expansion reduces the many-body problem into a one-dimensional coupled-channel problem. We have solved this problem using the Lagrange-mesh method, often with surprisingly small bases. Its convergence has been tested by comparison with the exact analytical ground-state energy for a contact interaction. We found that, in some cases, an extension of the method mixing different meshes can reduce computational times. The N = 3 case is treated separately with an exact evaluation of the K = 0 component. It easily allows a comparison with an exact variational calculation using the same Lagrange basis. Both methods converge very fast with the number of basis functions and agree.
The HH expansion was applied to treat several bosonic systems interacting with contact and Gaussian forces. It is shown that, for contact interactions, the lowest-order energy in the HH expansion has the correct N 3dependence at large N and is asymptotically higher than the known exact solution by a factor 48/16π ≈ 0.955.
The convergence on K max for these interactions is slow and can be described by an E(∞) + a 0 /(a 1 + K max ) dependence. This implies that, to achieve numerically converged values by the HH expansion, very large K values, of an order of a few hundreds, would be needed.
Numerical studies of the HH expansion, carried out for three, five, twenty and hundred bosons with a Gaussian interaction, have shown that their HH convergence rate is strongly determined by the Gaussian range. For a sufficiently large range, the K max = 0 energy can already be close to the converged value. For smaller ranges the convergence deteriorates and extrapolations are needed to get a converged value. The functional behaviour of E with respect to K max has been found to be reasonably described by an inverse power law at 20 ≤ K max ≤ 40, which makes it very difficult to obtain more than a few converged digits by further basis expansion, except for N = 3. As in the three-dimensional cases from [4,5], the convergence gets better with larger boson numbers N .

A Removing the Centre-of-Mass
To guarantee that centre-of mass excitations are absent in the wave function The operator R 2 in individual particle coordinates reads The calculation of its matrix elements in the oscillator basis involves one-body matrix elements of r 2 i , and two-body matrix elements of r i r j which include matrix elements of r ,

B Constructing States with Well-Defined Hypermomentum
According to Ref. [19] the angular part of the kinetic energy operator associated with N arbitrary coordinates If the first N − 1 of these coordinates are the Jacobi coordinates (1) and the last one is the coordinate of the centre-of-mass, x N = R, then Eq. (B.1) can be split into two parts, where the first term is the angular part of the kinetic energy operator associated with Jacobi coordinates, which is exactly the operator ρ , and the second one depends on the centre-of-mass coordinate, In order the wave function (0) κ K γ (R, ρ,ρ) to be the eigenfunction of the operator ρ , the expansion coefficients C μ K γ in (15) should be chosen as the eigenvectors of the matrix of this operator, calculated in the oscillator basis K μ (r 1 , .r 2 , . . . , r N ), corresponding to the eigenvalue K (K + n − 2). However, it is inconvenient to calculate matrix elements of the operator defined in Jacobi coordinates by (B.3) in a basis of states constructed in individual coordinates. Instead, we calculate K μ | N −1 | K μ as taking advantage that N is represented by N individual coordinates and that the diagonalization of the operator N −1 takes place after diagonalization of the matrix R 2 so that the new basis functions 0 K ν = μ C (0)ν K μ K μ do not contain centre-of-mass excitations. In this case, one has The operator N can be rewritten as To calculate K μ | N | K μ , one needs matrix elements of r 2 given by Eq. (A.1) and the matrix elements (C.7) The contribution from these terms into the hyperradial potential is obtained by applying the Laplace transform from Eq. (18), which gives (K + n 2 ) (K + n 2 ) ρ K +K +n−2 1 2πi i∞ −i∞ ds e sρ 2 s − K +K +n where 1 F 1 is the confluent hypergeometric function. At very small values of the range a the hypergeometrical function will the dominated by the first term of its asymptotic expansion, proportional to 1/ρ, which is exactly the same as the hyperradial potential of the contact interaction. With (2), the same integral is also evaluated in hyperspherical coordinates by Hence, one has