Non-Symmetrized Hyperspherical Harmonics Method for Non-Equal Mass Three-Body Systems

The non-symmetrized hyperspherical harmonics method for a three-body system, composed by two particles having equal masses, but different from the mass of the third particle, is reviewed and applied to the $^3$H, $^3$He nuclei and $^3_{\Lambda}$H hyper-nucleus, seen respectively as $nnp$, $ppn$ and $NN\Lambda$ three-body systems. The convergence of the method is first tested in order to estimate its accuracy. Then, the difference of binding energy between $^3$H and $^3$He due to the difference of the proton and the neutron masses is studied using several central spin-independent and spin-dependent potentials. Finally, the $^3_{\Lambda}$H hypernucleus binding energy is calculated using different $NN$ and $\Lambda N$ potential models. The results have been compared with those present in the literature, finding a very nice agreement.


I. INTRODUCTION
The hyperspherical harmonics (HH) method has been widely applied in the study of the bound states of few-body systems, starting from A = 3 nuclei [1,2]. Usually, the use of the HH basis is preceded by a symmetrization procedure that takes into account the fact that protons and neutrons are fermions, and the wave function has to be antisymmetric under exchange of any pair of these particles. For instance, for A = 3, antisymmetry is guaranteed by writing the wave function as p = 1, 2, 3 corresponding to the three different particle permutations [1]. However, it was shown in Refs. [3][4][5][6][7] that this preliminary step is in fact not strictly necessary, since, after the diagonalization of the Hamiltonian, the eigenvectors turn out to have a well-defined symmetry under particle permutation. In this second version, the method is known as non-symmetrized hyperspherical harmonics (NSHH) method. As we will also show below, the prize to pay for the non-antisymmetrization is that a quite larger number of the expansion elements are necessary with respect to the "standard" HH method. However, the NSHH method has the advantage to reduce the computational effort due to the symmetrization procedure, and, moreover, the same expansion can be easily re-arranged for systems of different particles with different masses. In fact, the steps to be done within the NSHH method from the case of equal-mass to the case of non-equal mass particles are quite straightforward and will be illustrated below. In this work, we apply the NSHH method to study the 3 H, 3 He and 3 Λ H systems, seen as nnp, ppn, N N Λ respectively (we used the standard notation of N for nucleon and Y for hyperon). In order to test our method, we study the first two systems listed above with five different potential models, and the hypernucleus with three potential models. We start with simple central spin-independent N N and Y N interactions, and then we move to central spin-dependent potentials. To be noticed that none of the interactions considered is realistic. Furthermore, we do not include three-body forces. Therefore, the comparison of our results with the experimental data is meaningless. However, the considered interactions are useful to test step by step our method and to compare with results obtained in the literature.
The paper is organized as follows: in Section II we describe the NSHH method, in Section III we discuss the results obtained for the considered nuclear systems. Some concluding remarks and an outlook are presented in Section IV.

II. THEORETICAL FORMALISM
We briefly review the formalism of the present calculation. We start by introducing the Jacobi coordinates for a system of A = 3 particles, with mass m i , position r i , and momentum p i . By defining x i = √ m i r i [8], they are taken as a linear combination of x i , i.e.
where the coefficients c ij need to satisfy the following conditions [8] 3 Here M is a reference mass. The advantage of using Eqs. (2)-(4) is that the kinetic energy operator can be cast in the form where y 3 is the center-of-mass coordinate. For a three-body system, there are three possible permutations of the particles. Therefore, the Jacobi coordinates depend on this permutations. For p = 3, i.e. i, j, k = 1, 2, 3, the Jacobi coordinates are explicitly given by They reduce to the familiar expressions for equal-mass particles when m 1 = m 2 = m 3 = M (see for instance Ref. [1]). We then introduce the hyperspherical coordinates, by replacing, in a standard way, the moduli of y 1,2 by the hyperradius and one hyperangle, given by tan φ (p) = y To be noticed that the hyperangle φ (p) depends on the permutation p, while the hyperradius ρ does not. The wellknown advantage of using the hyperspherical coordinates is that the Laplace operator can be cast in the form [8] where Λ 2 (Ω (p) ) is called the grand-angular momentum operator, and is explicitly written as Herel 2 1 andl 2 2 are the (ordinary) angular momentum operators associated with the Jacobi vectors y respectively, and Ω (p) ≡ (ŷ 2 , φ (p) ). The HH functions are the eigenfunctions of the grand-angular momentum operator Λ 2 (Ω (p) ), with eigenvalue −G(G + 4), i.e.
Here the HH function Y G (Ω (p) ) is defined as with N ℓ1,ℓ2 n a normalization factor [8] and ch Λ Σ G min 1 0 1/2 0 2 1 1/2 2 3 1 3/2 2 4 2 3/2 2 is the so-called grand-angular momentum. We remark that the HH functions depend on the considered permutation via Ω (p) . It is useful to combine the HH functions in order to assign them a well defined total orbital angular momentum Λ. Using the Clebsch-Gordan coefficients, we introduce the functions H G (Ω (p) ) as where [G] stands for [ℓ 1 , ℓ 2 , Λ, n], and We now consider our system made of three particles, two with equal masses, different from the mass of the third particle. We choose to fix the two equal mass particles in position 1 and 2, and we set the third particle with different mass as particle 3. Therefore, we will work with the Jacobi and hyperspherical coordinates with fixed permutation p = 3. The wave function that describes our system can now be cast in the form where u {G} (ρ) is a function of only the hyperradius ρ, and BH J {G} (Ω (3) ) is given by Eq. (14) multiplied by the spin part, i.e.
Here S is the spin of the first couple with third component s, Σ is the total spin of the system and Σ z its third component, and {G} now stands for {l 1 , l 2 , n, Λ, S, Σ}. To be noticed that the LS-coupling scheme is used, so that the total spin of the system is combined, using the Clebsh-Gordan coefficient (ΛΛ z , ΣΣ z |JJ z ), with the total orbital angular momentum to give the total spin J. Furthermore, (i) ℓ 1 , ℓ 2 and n are taken such that Eq. (13) is satisfied for G that runs from G min = ℓ 1 + ℓ 2 to a given G max , to be chosen in order to reach the desired accuracy, and (ii) we have imposed ℓ 1 + ℓ 2 = even, since the systems under consideration have positive parity. The possible values for Λ, Σ, and G min , which together with G max identify a channel, are listed in Table I for a system with J π = 1/2 + . Note that, since we are using central potentials, only the first channel of Table I will be in fact necessary.
In the present work, the hyperradial function is itself expanded on a suitable basis, i.e. a set of generalized Laguerre polynomials [1]. Therefore, we can write where c {G},l are unknown coefficients, and Here (5) L l (γρ) are generalized Laguerre polynomials, and the numerical factor in front of them is chosen so that f l (ρ) are normalized to unit. Furthermore, γ is a non-linear parameter, whose typical values are in the range (2 − 5) fm −1 .
The results have to be stable against γ, as we will show in Section III. With these assumptions, the functions f l (ρ) go to zero for ρ → ∞, and constitute an orthonormal basis. By using Eq. (18), the wave function can now be cast in the form where we have dropped the superscript (3) in Ω (3) to simplify the notation, and we have indicated with N max the maximum number of Laguerre polynomials in Eq. (18).
In an even more compact notation, we can write where Ψ ξ is a complete set of states, and ξ is the index that labels all the quantum numbers defining the basis elements. The expansion coefficients c ξ can be determined using the Rayleigh-Ritz variational principle [1], which states that where δ c Ψ denotes the variation of the wave function with respect to the coefficients c ξ . By doing the differentiation, the problem is then reduced to a generalized eigenvalue-eigenvector problem of the form that is solved using the Lanczos diagonalization algorithm [9]. The use of the Lanczos algorithm is dictated by the large size (∼ 50000 × 50000) of the involved matrices (see below). All the computational problem is now shifted in having to calculate the norm, kinetic energy and potential energy matrix elements. One of the advantage of using a fixed permutation is that the norm and kinetic energy matrix elements are or analytical, or involve just a one-dimensional integration. In fact, they are written as where J is the total Jacobian of the transformation, given by and f ′ l (ρ) and f ′′ l (ρ) are, respectively, the first and the second derivatives of the functions f l (ρ) defined in Eq. (19). The potential matrix elements in Eq. (23) can be written as with V ij indicating the two-body interaction between particle i and particle j. Note that in the present work we do not consider three-body forces. Since it is easier to evaluate the matrix elements of V ij when the Jacobi coordinate y 2 is proportional to r i − r j , we proceed as follows. We make use of the fact that the hyperradius is permutationindependent, and we use the fact that the HH function written in terms of Ω (p) can be expressed as function of the HH written using Ω (p ′ ) , with p ′ = p. Basically it can be shown that [1] H [G] (Ω (p) ) = where the grandangular momentum G and the total angular momentum Λ remain constant, i.e. G = G ′ and Λ = Λ ′ , but we have [G] = [G ′ ], since all possible combinations of ℓ 1 , ℓ 2 , n are allowed. The spin-part written in terms of permutation p can be easily expressed in terms of permutation p ′ via the standard 6j Wigner coefficients [10]. The transformation coefficients a can be calculated, for A = 3, through the Raynal-Revai recurrence relations [11]. Alternately we can use the orthonormality of the HH basis [1] Their explicit expression can be found for instance in Ref. [1] as is reported in the Appendix for completeness. The final expression for the potential matrix elements is given by It is then clear the advantage of using the NSHH method also for the calculation of the potential matrix elements, as in fact all what is needed is the calculation of one integral of the type with p the permutation corresponding to the order i, j, k.

III. RESULTS
We present in this section the results obtained with the NSHH method described above. In particular, we present in Section III A the study of the convergence of the method, in the case of the triton binding energy, calculated with m p = m n . We then present in Section III B the results for the triton and 3 He binding energy, when m p = m n . In Section III C we present the results of the hypertriton.
The potential models used in our study are central spin-independent and spin-dependent. In particular, the 3 H and 3 He systems have been investigated using the spin-independent Volkov [12], Afnan-Tang [13] and Malfliet-Tjon [14] potential models, and the two spin-dependent Minnesota [15] and Argonne AV4 ′ [16] potential models. Note that the AV4 ′ potential is a reprojection of the much more realistic Argonne AV18 [17] potential model. In the case of the hypernucleus 3 Λ H, we have used the Gaussian spin-independent central potential of Ref. [18], and two spin-dependent potentials: the first one, labeled MN9 [19], combines a Minnesota [15] potential for the N N interaction with the S = 1 component of the same Minnesota potential multiplied by a factor 0.9 for the ΛN interaction. The second one, labeled AU, uses the Argonne AV4 ′ of Ref. [16] for the N N interaction, and the Usmani potential of Ref. [20] for the ΛN interaction (see also Ref. [21]).

A. Convergence study
We recall that the wave function is written as in Eq. (20), and that, since we are using central potentials, only the first channel of Table I is considered, as for instance in Ref. [1]. Therefore we need to study the convergence of our results on G max and N max . Furthermore, we introduce the value of j as j = ℓ 2 + S, ℓ 2 and S being the orbital angular momentum and the spin of the pair ij on which the potential acts. This allows to set up the theoretical framework also in the case of projecting potentials. Therefore, we will study the convergence of our results also on the maximum value of j, called j max . Finally, the radial function written as in Eq. (19), presents a non-linear parameter γ, for which   we need to find a range of values such that the binding energy is stable. Note that in these convergence studies we have used m n = m p . We start by considering the parameter γ. The behaviour of the binding energy as a function of γ is shown for the Volkov potential in the top panel of Fig. 1. We mention here that for all the other potential models we have considered, the results are similar. The other parameters were kept constant, i.e. G max = 20, N max = 16 and j max = 6. The particular dependence on γ of the binding energy, that increases for low values of γ, is constant for some central values, and decreases again for large values of γ, allows to determine a so-called plateau, and the optimal value for γ has to be chosen on this plateau. Alternatively, we can chose γ such that for a given N max the binding energy is maximum. A choice of γ outside the plateau would require just a larger value of N max . To be noticed that this particular choice of γ is not universal. As an example, in the "standard" HH method, γ = 2.5 − 4.5 fm −1 for the AV18 potential, but much larger (≃ 7 fm −1 ) for the chiral non-local potentials [1]. In our case, different values of γ for different potentials might improve the convergence on N max , but not that on j max and G max , determined by the structure of the HH functions. Since, as shown below, the convergence on N max is not difficult to be achieved, we have chosen to keep γ at a fixed value, i.e. γ = 4 fm −1 for all the potentials.
In the bottom panel of Fig. 1 we fix j max = 8, γ = 4 fm −1 and G max = 20, and we show the pattern of convergence for the binding energy B with respect to N max , in the case of the Argonne AV4 ′ model. Here convergence is reached for N max = 24, i.e. we have verified that, for higher N max value, B changes by less than 1 keV. To be noticed that for the other potentials, convergence is already reached for N max = 16 − 20.
The variation of the binding energy as a function of j max and G max depends significantly on the adopted potential model. Therefore, we need to analyze every single case. As we can see from the data of Tables II and III, the convergence on j max and G max for the Volkov and the Minnesota potentials is really quick, and we can reach an accuracy better than 2 keV for j max = 10 and G max = 40. On the other hand, in the case of the Afnan-Tang potential, we need to go up to j max = 14 and G max = 50, in order to get a total accuracy of our results of about 2 keV (1 keV is due to the dependence on N max ). This can be seen by inspection of Table IV. The Malfliet-Tjon potential model implies a convergence even slower of the expansion, and we have to go up to G max = 90 and j max = 22, to get an uncertainty of about 3 keV, as shown in Table V. In fact, being a sum of Yukawa functions, the Malfliet-Tjon potential model is quite difficult to be treated also with the "standard" symmetrized HH method [1].
In Tables VI and VII, we show the convergence study for the AV4 ′ , which is the most realistic potential model used   for the Volkov potential model [12], with G max = 20, jmax = 6 and Nmax = 16, and using mn = mp. Bottom panel: The binding energy B (in MeV) as function of the parameter Nmax for the AV4 ′ potential model [16], with G max = 20, jmax = 8 and γ = 4 fm −1 , and using mn = mp.
here for the A = 3 nuclear systems. As we can see by inspection of the tables, in order to reach an accuracy of about 3 keV, we have to push the calculation up to G max = 80, N max = 24 and j max = 20. Our final result of B = 8.991 MeV, though, agrees well with the one of Ref. [22], obtained with the "standard" symmetrized HH method, for which B = 8.992 MeV.
The results for the binding energy of 3 H and 3 He with the different potentials will be summarized in the next Subsection.  [16] as function of Nmax, jmax and G max , using mn = mp and γ = 4 fm −1 .

B. The 3 H and 3 He systems
Having verified that our method can be pushed up to convergence, we present in the third column of Table VIII the results for the 3 H binding energy with all the different potential models, obtained still keeping m p = m n . The results are compared with those present in the literature, finding an overall nice agreement.
We now turn our attention to the 3 H and 3 He nuclei, considering them as made of different mass particles. Therefore, we impose m p = m n and we calculate the 3 H and 3 He binding energy and the difference of these binding energies, i.e.
To be noticed that we have not yet included the effect of the (point) Coulomb interaction. The results are listed in   verified that ∆B is equally distributed, i.e. we have verified that as can be seen from Table VIII. We would like to remark that in the NSHH method, the inclusion of the difference of masses is quite straightforward, and ∆B can be calculated "exactly". This is not so trivial within the symmetrized HH method. Furthermore, we compare our results with those of Ref. [23], where ∆B was calculated within the Faddeev equation method using realistic Argonne AV18 [17] potential, and it was found ∆B = 14 keV, in perfect agreement with our AV4 ′ result. In order to test our results for ∆B, we try to get a perturbative rough estimate of ∆B, proceeding as follows: since the neutron-proton difference of mass ∆m = m n − m p = 1.2934 MeV is about three orders of magnitude smaller than their average mass m = (m n + m p )/2 = 938.9187 MeV, we can assume also ∆B to be small. Furthermore, we suppose the potential to be insensitive to ∆m, and we consider only the kinetic energy. In the center of mass frame, the kinetic energy operator can be cast in the form where m e stands for the mass of the two equal particles, i.e. m n for 3 H and m p for 3 He, and m d is the mass of the third particle, different from the previous ones. By defining E = H = T + V , where H is the average value of the Hamiltonian H, we obtain where we have indicated T e/d ≈ p 2 i /(2m e/d ). Moreover, we define the proton and neutron mass difference ∆m p/n as and the 3 He and 3 H binding energy difference ∆B3 He/ 3 H as (40) In conclusion where the last equality holds assuming that T e = T d = T /3, since the 3 He and 3 H have a large S-wave component (about 90 %). The results of T and ∆B P T are listed in Table IX, and are compared with the values for ∆B calculated within the NSHH and already listed in Table VIII. By inspection of the table we can see an overall nice agreement between this rough estimate and the exact calculation for all the potential models. Only in the case of the Minnesota and AV4 ′ potentials, ∆B P T is 4 and 3 keV larger than ∆B, respectively. This can be understood by noticing that these potentials are spin-dependent, giving rise to mixed-symmetry components in the wave functions. These components are responsible for a reduction in ∆B P T [24], related to the fact that the nuclear force for the 3 S 1 np pair is stronger than for the 1 S 0 nn (or pp) pair. Therefore, the kinetic energy for equal particles T e is less than the kinetic energy for different particles T d .

C. The 3 Λ H hypernucleus
The hypernucleus 3 Λ H is a bound system composed by a neutron, a proton and the Λ hyperon. In order to study this system, we have considered the proton and the neutron as reference-pair, with equal mass m n = m p = m, while the Λ particle has been taken as the third particle with different mass. The Λ hyperon mass has been chosen depending on the considered potential. We remind that we have used three different potential models: a central spin-independent Gaussian model [18], and two spin-dependent central potentials, labelled MN9 [19] and AU [21] potentials. Therefore, when the 3 Λ H hypernucleus has been studied using the Gaussian potential of Ref. [18], we have set M Λ = 6/5 m N , accordingly. In the other two cases, we have used M Λ = 1115.683 MeV. We first study the convergence pattern of our method, which in the case of the Gaussian potential of Ref. [18] is really fast, with a reached accuracy of 1 keV on the binding energy already with N max = 20, j max = 10 and G max = 50. This can be seen directly by inspection of Tables X and XI.
The convergence pattern in the case of the spin-dependent central MN9 and AU potentials has been found quite slower. This is shown in Tables XII and XIII, respectively. By inspection of Table XII, we can conclude that B = 2.280 MeV, with an accuracy of about 3 keV, obtained with G max = 100, N max = 34, and j max = 14. By inspection of Table XIII The results obtained with our method for the three potential models considered in this work are compared with those present in the literature [18,19,21] in Table XIV, finding a very nice agreement, within the reached accuracy.    [18], using G max = 20, jmax = 8 and γ = 4 fm −1 .

IV. CONCLUSIONS AND OUTLOOK
In this work we present a study of the bound state of a three-body system, composed of different particles, by means of the NSHH method. The method has been reviewed in Section II. In order to verify its validity, we have started by considering a system of three equal-mass nucleons interacting via different central potential models, three spin-independent and two spin-dependent. We have studied the convergence pattern, and we have compared our results at convergence with those present in the literature, finding an overall nice agreement. Then, we have switched on the difference of mass between protons and neutrons and we have calculated the difference of binding energy ∆B due to the difference between the neutron and proton masses. We have found that ∆B depends on the considered potential model, but is always symmetrically distributed (see Eq. (33)).
Finally we have implemented our method for the 3 Λ H hypernucleus, studied with three different potentials, i.e. the Gaussian potential of Ref. [18], for which we have found a fast convergence of the NSHH method, the MN9 and the AU potentials of Ref. [19], for which the convergence is much slower. In these last two cases, in particular, we had found necessary to include a large number of the HH basis (46104 for the MN9 and 52704 for the AU potentials), but the agreement with the results in the literature has been found quite nice. To be noticed that we have included only two-body interactions, and therefore a comparison with the experimental data is meaningless.