The Tucson-Melbourne Three-Nucleon Force in the automatized Partial Wave Decomposition

A recently developed procedure for a partial wave decomposition of a three-nucleon force is applied to the pi-pi, pi-rho and rho-rho components of the Tucson-Melbourne three-nucleon potential. The resulting matrix elements for the pi-pi and pi-rho components are compared with the values obtained using the standard approach to the partial wave decomposition, in which the pi-rho expressions for the matrix elements are also derived and presented. Several numerical tests and results for the triton binding energy and the correlation function prove the reliability and efficiency of the new method.


I. INTRODUCTION
The Tucson-Melbourne (TM) three nucleon force (3NF) [1][2][3][4] is an important model of the three-nucleon (3N) interaction. It consists of three parts stemming from exchanges of π −π , π − ρ and ρ − ρ mesons. The main ingredient of the TM force, the meson-nucleon scattering amplitude with the off-shell mesons, was derived using the current algebra techniques. This was done in [1] and improved in [3] for the π − π part. The π − ρ and ρ − ρ contributions were derived in [2,4]. In [5] the structure of the π − π part of the TM 3NF was revisited to achieve a consistency with the chiral symmetry and the modified force is known as the TM' model.
The effects of all terms on the triton binding energy were studied in [6]. It turned out that the π − ρ force acts repulsively for the 3 H contrarily to the π − π interaction and combining them leads to the 3 H binding energy close to the experimental value. The ρ − ρ force has only a small influence on the triton binding energy. A similar behaviour was observed for scattering observables in the three-nucleon system [7]: the largest effects came from the dominant π − π part and the influence of the π − ρ part was smaller and in the opposite direction. The ρ − ρ contribution proved to be much smaller and practically negligible. However, the results of Refs. [6,7] were based on partial waves restricted to the total angular momenta in the two-nucleon subsystem j ≤ 2. Thus conclusions of [7] are valid only in a low energy domain of the three-nucleon continuum. For higher energies, where more partial waves are required to achieve convergence, only the π − π part of the TM force was used, (see e.g. [8,9]). While the inclusion of this main component of the TM 3NF improves the description of many scattering observables, some serious discrepancies with data remain and they become larger at higher energies. One of the possible explanations for this disagreement is a lack of shorter-range parts of the 3NF in those calculations, what calls for a reliable and fast method to obtain matrix elements for all components of the TM force in higher partial waves.
Recently we have proposed a novel, automatized way to perform a partial wave decomposition of any two-and three-nucleon potential [10]. This approach makes use of a software for symbolic calculations to generate the part of the code which is specific for a considered force model. More precisely, in this way we calculate exactly the isospin and spin-momentum parts of the nuclear interactions and generate a corresponding FORTRAN (or C) code. That momentum dependent output forms an integrand for further five-dimensional numerical integrations.
In this paper we present results of applying that new scheme to the original TM 3N force. They confirm the feasibility and efficiency of our method and its numerical implementation. The existence of such a reliable procedure is especially important in view of available and forthcoming results from the chiral perturbation theory (χPT) [11] for 3N forces at higher orders of the chiral expansion. A big number of different momentum-spin-isospin structures contained in those interactions requires a safe and automatized method to perform partial wave decompositions, which is guaranteed by our method. Furthermore our scheme avoids the application of partial wave decomposed permutation operators when dealing with products of 3NF's and permutation operators as they are often required, e.g., in 3N Faddeev equations. Such an application is numerically demanding because it requires a huge number of partial waves. Thus, again an efficient, fast and precise method is needed.
Our novel scheme of an automatized partial wave decomposition (aPWD) is described in Sec. II. Results and additional tests for our numerical realization are presented in Sec. III and conclusions are given in Section IV. The standard PWD of the π − ρ component of the TM 3NF is given in the Appendix.

II. AUTOMATIZED PARTIAL WAVE DECOMPOSITION
The 3NF, V 123 , is an indispensable ingredient in a theoretical description of the few-body systems. It can be always written as a sum of three terms where each V (i) is symmetrical under the exchange of nucleons j and k (i, j, k = 1, 2, 3, i = j = k). Such a splitting in the case of the π − π exchange TM 3NF corresponds to the possible choices of the nucleon undergoing off-shell πN scattering. The 3NF typically enters the dynamical equations via its part V (1) . In the case of the three-nucleon bound state, the Faddeev component ψ fulfils the following equation [12] where G 0 is the free 3N propagator and t is the two-body t-operator generated from a given nucleon-nucleon (NN) potential through the Lippmann-Schwinger equation. The permutation operator P ≡ P 12 P 23 + P 13 P 23 is given in terms of the transpositions P ij , which interchange particles i and j. The full bound state wave function Ψ is then obtained as Ψ = (1 + P )ψ. Transition amplitudes for the elastic nucleon-deuteron scattering, U, and for the breakup reaction, U 0 , are given as [13] where the auxiliary state T fulfils the 3N Faddeev equation with Φ being the initial state composed of the deuteron wave function and a momentum eigenstate of the projectile nucleon. Equations (2.2) and (2.4) are solved [13,14] in the momentum space using 3N partial-wave states | p, q, α in the jJ-coupling [14,15] | p, q, α ≡| pq(ls)j(λ where p and q are magnitudes of the standard Jacobi momenta and α denotes a set of discrete quantum numbers arising in the following way: the spin s of the subsystem composed from nucleons 2 and 3 is coupled with their orbital angular momentum l to the total angular momentum j. The spin 1 2 of the spectator particle 1 couples with its relative orbital angular momentum λ to the total angular momentum of nucleon 1, I. Finally, j and I are coupled to the total 3N angular momentum J with the projection M J . For the isospin part, the total isospin t of the (23) subsystem is coupled with the isospin 1 2 of the spectator nucleon to the total 3N isospin T with the projection M T .
Any three-nucleon force enters Eqs. (2.2)- (2.4) in the form of V (1) (1 + P ). Therefore a partial wave decomposition of V (1) as well as V (1) P has to be performed. The standard approach to perform a partial wave decomposition of V (1) [16] is very tedious, even with improvements suggested in [17], since each momentum-spin-isospin structure, which occurs in a 3NF, has to be treated separately. In the case when a 3NF consists of a big number of such structures, like chiral 3NF's at higher orders of the chiral expansion, the traditional approach to a partial wave decomposition is very inefficient and extremely time consuming. In addition, the application of the permutation operator, when calculating V (1) P , causes an additional numerical problem, which originates from a slow convergence of the V (1) P matrix elements with respect to the number of intermediate states |α ′′ >: In order to calculate precisely these matrix elements, a big number of intermediate states | α ′′ is required, and, thus, one is forced to calculate matrix elements of V (1) operator for a much bigger set of α ′′ states than actually needed in order to get converged solutions of the Faddeev equations. In our new approach, called in the following automatized partial wave decomposition (aPWD), to get matrix elements of V (1) and V (1) P , that drawback is removed because matrix elements of V (1 ) and V (1) P are calculated directly.
The starting point of our method is the observation, that any 3N interaction and thus also its V (1) component in momentum space can be written as a sum of terms in the form whereÔ spin andÔ isospin are operators acting on spin and isospin degrees of freedom, respectively, which are built from the spin ( σ i ) and isospin ( τ i ) operators of individual nucleons. Scalar factors f ( q 1 , q 2 , q 3 ) and spin operatorsÔ spin ( q 1 , q 2 , q 3 , σ 1 , σ 2 , σ 3 ) depend on momentum transfers q i to nucleon i which are expressed in terms of the initial and final Jacobi momenta p, q and p ′ , q ′ , respectively, as: For example, in the π − π part of the TM 3N force one meets the following spin-isospin structures: Note, that not all combinations ofÔ spin andÔ isospin actually appear in the above example.
In the first step of aPWD we calculate 3NF matrix elements using partial wave states | p, q, β [15] in the so-called LS-coupling | p, q, β 1 ≡| pq(lλ)L(s where the relative orbital angular momentum l (within the pair (23)) and λ (between the pair (23) and nucleon 1) are coupled to the total orbital angular momentum L. In the spin space, the spin of the (23) pair is coupled with the spin 1 2 of the nucleon 1 to the total spin S. Finally L and S are coupled to the total 3N angular momentum J with the projection M J . The index 1 emphasizes, that the spectator particle is nucleon 1. β describes the set of discrete quantum numbers discussed above. The isospin state is the same as in the basis state |p, q, α .
In this basis, it is easy to decouple the isospin and spin parts from the momentum part, what leads to the following form of a 3NF matrix element with the standard Clebsch-Gordan coefficients and the spherical harmonics. For abbreviation we skip in (2.11) and in the following the spin σ i and the isospin τ i operators in the arguments ofÔ spin andÔ isospin operators. The matrix element in the spin space appearing in (2.11), (s ′ 1 2 )S ′ M J − m L ′ |Ô spin ( p ′ , q ′ , p, q) | (s 1 2 )S M J − m L , depends on the momenta q i and spin quantum numbers. Using a software for symbolic calculations (such as Mathematica c [18] in our case) it is very easy to calculate this matrix element for all combinations of spin quantum numbers as a function of the momentum vectors q i . To this aim we use the Kronecker product built in Mathematica, which allows us to express the spin matrix element in terms of simple matrix operations. This is even more straightforward in the case of the isospin matrix element, which does not depend on any additional parameters. Another advantage of using software for symbolic calculations is the possibility to generate a Fortran (or C) code in an automatized way. This eliminates possible errors which can be introduced during programing of very lengthy formulas for the spin matrix element. The calculation of the 3NF matrix elements requires finally an eight-dimensional integration shown in (2.11). In a typical case the total isospin and its projection is conserved. We also assume that the considered 3N force is rotationally invariant. Then the matrix elements in (2.11) vanish unless J = J ′ and M J = M J ′ , and, additionally, do not depend on M J . Thus we can calculate which is equal to the original matrix element of V (1) given in Eq. (2.11). The integrand in G(l ′ , λ ′ , L ′ , s ′ , S ′ , t ′ , l, λ, L, s, S, J, t, T, M T ): is a scalar and thus does not depend on all directions of the Jacobi momenta [19]. Therefore we are free to choose for example p along z-axis ( p = (0, 0, p)) and φ q = 0 and are left with five-fold integrations only The reduction of the number of integrations for a simple example of 3NF is numerically exemplified in Ref. [10]. The remaining summations over m L ′ , m L and M J and 5-fold integrations can be performed for a small number of (p,q,p',q') combinations even on a personal computer. However, a large number of five-dimensional integrations, as they are needed to obtain all matrix elements needed for the solution of the 3N Faddeev equations, has to be carried out on a powerful parallel computer. Once the matrix elements p ′ , q ′ , β ′ | V (1) | p, q, β are calculated, recoupling to the jI-representation, p ′ , q ′ , α ′ | V (1) | p, q, α , can be easily performed [15] Now let us turn to the V (1) (1 + P ) operator and discuss its V (1) P 12 P 23 matrix element where P spin ij (P isospin ij ) is the part of P ij operator acting in the spin (isospin) space, one gets Similarly, for V (1) P 13 P 23 one gets That means that the calculation of these two contributions proceeds in the same way as calculation of the V (1) matrix element. Only the arguments of the termÔ spin have to be changed and additional factors originating from the recoupling of the spin and isospin quantum numbers have to be taken into account. As for the V (1) operator also here the eight-fold integrations can be reduced to the five-fold ones and recalculation to |p, q, α states can be performed.
It is important to note that, since our basis states |p, q, α are antisymmetric with respect to the exchange of nucleons 2 and 3, (2.19) and (2.20) yield the same values for the matrix elements. This allows one to reduce significantly the size of the codes and the required computation time.

III. RESULTS
A. The TM 3NF and its π − π , π − ρ , and ρ − ρ components Since the aim of this work is not to study the dependence of the matrix elements of the TM force on its parameters, in the following we use their values given in Tab.I of Ref. [4]: In the numerical implementation of (2.11) we use the same number of gaussian points for each of the five angular domains. It might be more efficient to relax this constraint in future applications and to optimize the grids further. Thus our integration method leaves room for improvement, even if we will later demonstrate in subsection III G that it leads to fully converged results.
The TM 3NF matrix elements calculated in the basis (2.5) are functions of four momentum magnitudes and two sets of discrete quantum numbers. In Figs. 1-2, examples of the TM force V (1) matrix elements are shown together with its π − π, π − ρ and ρ − ρ components in one-dimensional plots. In Fig. 1, the matrix elements < p ′ , q ′ , α ′ |V (1) |p, q, α > for p ′ = q ′ = q = 0.132 fm −1 and for different channel pairs (α ′ , α) (see Tab. I) are shown as a function of the momentum p. The same matrix elements but for the momenta p ′ = 0.711 fm −1 , q ′ = 0.132 fm −1 , and q = 2.84 fm −1 are shown in Fig. 2 again as a function of p. The π − π part dominates in all cases but the π − ρ part is also important (see Fig. 1b,2a-c). The ρ − ρ part is of less importance for all the considered matrix elements.
B. The aPWD for V (1) (1+P) operator As was described in Sec. II, aPWD can be applied not only to the V (1) alone but also to the V (1) (1 + P ) operator. Using aPWD for V (1) (1 + P ) has the same advantages as for the V (1) operator: the automatized procedure can be easy tuned to any kind of 3NF and reduces the possibility of errors. In the current implementation of aPWD the calculation of  Fig. 3 and (1,1) and (6,8) in Fig. 4, where π − π force dominates, the picture is similar to the corresponding ones in Figs. 1 and 2. For the remaining channel combinations the differences are more visible, for example the inclusion of the permutation operator for the π − π component for the (6,8) pair in Fig. 3 leads to the change of the sign and strength of this force. In that case also the π − ρ part becomes bigger after the permutation operator is applied. Also for the (6,3) case in Figs. 3 and 4 the action of the permutation operator changes the strength of the matrix element and increases the momentum range, where both π − π and π − ρ components play a significant role. For all the here presented cases the ρ − ρ force is much smaller than the remaining interactions. The aPWD method allows us to study the role played by different isospin structures entering the TM force. An example is given in Fig. 5 where, for the π − ρ force, the contribution from the so-called "Kroll-Ruderman" and two "∆" terms [4] (see also Appendix A 1) are shown. For the presented matrix elements (< p ′ = 0.132f m −1 , q ′ = 0.132f m −1 , α ′ = 1|V (1) π−ρ (1 + P )|p, q = 0.132f m −1 , α = 1 >) the "Kroll-Ruderman" term dominates for small momenta p, while the two "∆" terms are bigger for p > 2 fm −1 . However, they have opposite signs, so their combined effect is weak and leads to a reduction of the strength of the dominant "Kroll-Ruderman" term.
C. The comparison of the standard and automatized PWD schemes for π − π and π − ρ forces.
For the π − π force the partial wave decomposition has been presented in [3] and in an alternative way in [16]. The comparison of results obtained by aPWD and the ones obtained in Ref. [16] is presented in Fig. 6. Again the channel pairs and momenta are chosen as in Fig. 1. A very good agreement between the both methods is clearly seen.
In Appendix A we present expressions for the partial wave decomposition of the π − ρ force. This decomposition is in the spirit of the decomposition of the π−π interaction given in Ref. [16]. In Fig. 7 we compare results obtained in the aPWD scheme with those based on PWD given in Appendix A. Because of the internal construction of the PWD from Appendix A, we compare matrix elements of V (1) P 13 P 23 instead of V (1) . The matrix elements of the standard PWD are obtained using partial waves up to j max = 5 in intermediate states. For this truncation, the matrix elements considered here are converged (see Sec. III E). Again, for all given examples, the agreement between both methods is excellent.
Though the CPU time is smaller for the scheme presented in Appendix A, the long time which is needed for the derivation of the partial wave decomposition of complicated spinmomentum structures and its programming in the standard way is incomparable with the relatively short time demanded by aPWD. Another advantage of aPWD lies in its flexibility which allows one to use it easily for different operators. In the case of the standard PWD each spin-momentum structure has to be treated separately.
D. The equality of V (1) P 12 P 23 and V (1) P 13 P 23 The equality of V (1) P 13 P 23 and V (1) P 12 P 23 matrix elements between the states antisymmetrized in the (23) subsystem forms another nontrivial test of numerics. To check this, we compare some matrix elements for V (1) P 12 P 23 obtained via Eq.(2.21) with the corresponding ones for V (1) P 13 P 23 from Eq.(2.22). Results are displayed in Fig.8 again for four combinations of channel pairs and selected values of p ′ , q ′ and q momenta (the same as in Fig. 1). The numerical confirmation of the equality of the V (1) P 13 P 23 and V (1) P 12 P 23 matrix -5.0×10 The aPWD result for the V (1) (1 + P ) operator, which corresponds to the infinite number of the intermediate partial waves taken into account during the action of the permutation operator, gives the limit to which results of the traditional scheme should converge. This convergence is confirmed in Figs. 9 and 10 for the π − ρ part of the TM and the full TM 3NF, respectively. The channels and momenta are the same as in Fig. 2. While for the channel combination (1,1) already the smallest number of partial waves gives the aPWD limit, for the other combinations much more partial waves have to be taken into account. For one of the shown here cases (Fig. 10c) taking all partial waves up to j max = 5 is still insufficient to achieve the limit of aPWD. Note, however, that the magnitude of this matrix element is relatively small. In general, the convergence of the traditional PWD scheme is fully confirmed. <p',q',α'| V (1) π−ρ (1+P)| p,q,α> [fm 5 ] FIG. 5: (Color online) The contributions from the different parts of the π − ρ force for matrix elements < p ′ = 0.132 fm −1 , q ′ = 0.132 fm −1 , α ′ = 1|V 1 π−ρ (1 + P )|p, q = 0.132 fm −1 , α = 1 >. The black solid line represents the total π − ρ TM 3NF while the red dotted, green dashed and blue dot-dashed lines represent the "Kroll-Ruderman", the isospin even ∆ and the isospin odd ∆ terms, respectively.
F. The binding energy and correlation function for 3 H.
As a first application, we would like to calculate in the following the binding energy of 3 H, some energy expectation values and the correlation function. The obtained binding energies and expectation values of the kinetic energy H 0 , the NN potential energy V N N and the 3N potential energy V 3N are given in Tab. II for several realistic NN interactions alone and together with the TM force. The TM force was included for all states with subsystem total angular momentum j ≤ 2. The inclusion of the TM force leads to a stronger binding of 3 H. The binding energy changes, after the inclusion of the TM force, by approximately -1.093 MeV for the CDBonn potential and from -1.122 to -1.334 MeV for Nijmegen potentials. These results are in a reasonable agreement with the ones given in Tab. 2 of Ref. [6] for the Bonn OBEPQ (-9.596 MeV) and the Nijmegen (-8.689 MeV) potentials. Note, that in Ref. [6] slightly different values of the a, b and c parameters were used. In our calculations, we include partial waves up to j max = 5 for the two-body interaction. This is also different from Ref. [6] where only partial waves up to j max = 2 were included. Of course, for the given set of parameters, the binding energies do not accurately reproduce the experimental value of -8.482 MeV. Because of the well-known scaling behavior of many N-d scattering observables with the triton binding energy (see for example [20]), it will be necessary to finetune the TM model such that the triton binding energy is more accurately reproduced, e.g. along the lines of Ref. [12] The inclusion of the TM 3NF leads for all the NN potentials to higher expectation values of the kinetic energy and lower expectation values of the NN potential energy (about 3-6 MeV). The expectation values of the 3N potential energy amounts from 3.3% to 4.8% of the expectation values of the NN potential, depending on the particular NN potential. This observations are in line with the general expectations for the strength of 3NF's and the more compact state of 3 H when the binding energy is increased. The correlation function is defined in the configuration space as [12] C(r) ≡ 1 3 where r ij is the relative distance operator conjugate to the operator of the Jacobi momentum p. It is shown in Fig. 11 for the different NN potentials alone and combined with the TM 3NF. For the smaller distances shown in Fig. 11, the probability to find two nucleons increases when the TM 3NF is included. At least in part, this can be understood because the correlation functions drop more quickly for larger r due to the increased binding energy. Note that at short distances, the effect of the 3NF's is much larger than the dependence on the NN interaction model. The here presented correlation functions are in good agreement with the ones presented in [12] for the same NN potentials combined with the π − π part of the TM force.
G. The quality of the five-dimensional integration Finally, we would like to give an example of the stability of aPWD against the number of points used in the numerical integrations. In Tab. III the V (1) (1 + P ) matrix elements are NN potential  given for the same channels and momenta as in Fig. 1 for two values of momentum p=0.711 fm −1 and p=5.959 fm −1 . Results were obtained using N=12 or N=15 gaussian points in each of the five integrations in (2.11). The agreement seen in Tab. III between both predictions clearly demonstrates that the numerical integration is well under control and leads to fully converged numbers.   (1) (1 + P ) matrix elements (in fm 5 ) for channels combinations and momenta p ′ , q ′ , q as in Fig. 1, depending on the number of gaussian points N used in the five-fold integration of Eq.(2.14). The value of momentum p is 0.711 fm −1 (left) and 5.959 fm −1 (right).

IV. SUMMARY
We apply an automatized method of partial wave decomposition to the Tucson-Melbourne three-nucleon force. The obtained results agree very well with the traditional way of a partial wave decomposition for π − π and π − ρ contributions to the TM 3NF. For the latter one, we also give formulas of the partial wave decomposition in the traditional approach. Matrix elements obtained in the new way are used in the calculations of the triton wave function with different underlying nucleon-nucleon potentials. We performed also different numerical tests, which confirm the reliability of our method and computer codes.
Among many advantages of aPWD, we would like to emphasize its generality, efficiency, the semi-automatized process of preparing a code and the possibility of a calculation of the higher partial waves. The latter point gives hope for the future use of the full Tucson-Melbourne force in a description of 3N scattering at higher energies. The expected strong effects on observables coming from a 3NF should be tested also for short-range parts of three-body interactions. Such parts are included in the Tucson-Melbourne force.
The automatized partial wave decomposition is especially important in view of future applications of 3NFs arising from the χPT. In this approach, consistent two-and threebody forces are derived [11]. The numerous spin-momentum and isospin structures, which occur at higher orders of the chiral expansion require an efficient and automatized method for the PWD. The here presented results for the TM force prove, that such a method already exists. and The R πρ KR , R πρ ∆ + , and R πρ ∆ − form factors are given in terms of regularization form factors at the meson-baryon-baryon vertices F i as [6] R πρ KR (q 2 , q ′ 2 ) = and The F i are taken in monopole form with i = {πNN, πN∆, ρNN D , ρNN P , ρN∆}. The mass of the boson at the corresponding vertex, m b , is either m π or m ρ with exception of the case i = ρN∆ when m b = 0. We would like to have matrix elements of the three-body force in a partial wave basis |pqα 1 , where p and q are magnitudes of Jacobi momenta ( p is the relative momentum between particles 2 and 3 and q is the momentum of spectator particle 1 relative to the 2-3 pair) and α denotes discrete quantum numbers which we separate in spin, α J , and isospin, α T , parts The partial wave states corresponding to different spectator nucleon i (i = 1, 2, 3) can be obtained from |pqα 1 acting with proper permutation operator, for instance P 13 P 23 |pqα 1 = |p ′ q ′ α ′ 3 . According to the scheme presented in Ref. [16] these matrix elements can be calculated as where spin operators F and form factors R are taken among F KR , F I ∆ + , F II ∆ + , F I ∆ − , F II ∆ − , and R πρ KR , R πρ ∆ + and R πρ ∆ − , respectively, for different contributing terms. The matrix elements of the isospin parts appearing in Eq. (A3) are given [16] by where we use abbreviationâ ≡ 2a + 1.

The term F
The B Λ πN∆ and B ΛρN ∆ are given by The summation over µ in Eq. (A21) can be carried through resulting in