Orthogonal polynomial approach to calculate the two-nucleon transition operator in three dimensions

We give a short report on the possibility to use orthogonal polynomials (OP) in calculations that involve the two-nucleon (2N) transition operator. The presented work adds another approach to the set of previously developed methods (described in Phys. Rev. C 81, 034006 (2010); Few-Body Syst. 53, 237 (2012); K. Topolnicki, PhD thesis, Jagiellonian University (2014)) and is applied to the transition operator calculated at laboratory kinetic energy 300MeV. The new results for neutron-neutron and neutron-proton scattering observables converge to the results presented in Few-Body Syst. 53, 237 (2012) and to results obtained using the Arnoldi algorithm (Y. Saad, Iterative methods for sparse linear systems (SIAM Philadelphia, PA, USA 2003)). The numerical cost of the calculations performed using the new scheme is large and the new method can serve only as a backup to cross-check the previously used calculation schemes.


Introduction
The transition operatorť satisfying the Lippmann-Schwinger equation (LSE), contains all the necessary information to calculate any observable in the nucleon-nucleon scattering process [1]. It is also an important dynamical ingredient in momentum space few-nucleon calculations [2][3][4]. This important object was a subject of our previous work [2,3,5] and in this paper we apply an alternative approach to solve the LSE that complements the previous methods. In eq. (1)V is the 2N potential,Ǧ 0 (E + i ) = (E + i − p 2 m ) −1 is the free propagator, E is the energy, m is the nucleon mass and p is the relative momentum between the nucleons.
The traditional approach to solve eq. (1) requires the partial wave decomposition of all relevant operators. This standard technique has been very successful in describing experimental data at lower energies, however at higher energies, sometimes, a very large number of partial waves needs to be taken into account in order to achieve convergence for certain observables [2,6,7]. In [2,3,5] we describe and utilize an alternative approach that uses directly the three-dimensional (3D) degrees of freedom of the nucleon. In this method, where we work in the isospin framework, a e-mail: kacper.topolnicki@uj.edu.pl a general parity, time reversal and rotation invariant form of the 2N potential [8] is used: with |γ being one of the four possible isospin states of the 2N system, | p being a relative momentum eigenstate, w i (p , p) are known 2N spin operators (they are listed for example in [2]) and v γ i (|p |, |p|,p ·p) are scalar functions that define the potential. The transition operator has a similar operator form: but with energy-dependent scalar functions t γ i (E; |p |, |p|, p ·p).
After inserting (2) and (3) into the LSE (1) the spin dependencies can be removed. What remains is a set of equations that involve the scalar functions v γ i (|p |, |p|,p ·p) and t γ i (E; |p |, |p|,p ·p). A further observation that the resulting equations constitute a linear system with operators acting on the scalar functions makes it possible to use Krylev subspace algorithms (for example the Arnoldi algorithm [9]) to calculate the solution. A careful consideration of the resulting relations leads to some additional simplifications. It can be easily verified that every total isospin state γ, initial relative momentum magnitude |p| and energy E case is governed by an independent equation and can be treated separately. It is therefore possible, for each independant case, to consider only two-dimensional functions t: and the resulting equation has a form that is very similar to the LSE, but this time with operators acting on the functions (vectors) t, The linear operatorB is defined for example in [5]. For negative energies E it transforms the scalar function t into a new function t according to the following: where the momentum p = |p |( √ 1−x 2 cos φ , √ 1−x 2 sin φ , x ) and the expressions for the functions B kjj (|p |, |p|, x , |p |, x , φ ) are listed, e.g., in [5]. For nucleonnucleon scattering and positive energies, eq. (6) has to be reformulated to treat the singularity of the free propagator. This can be performed with the use of the standard methods with the principle value of the integral [2,3,5].
The numerical solution of (5) for the scalar functions "t" that can be used to reproduce the operator form of the transition operator (3) can be worked out using various approaches. One approach requires the construction of the matrix representation of the operatorB directly from (6). This is a difficult task and the resulting matrices are large making this method unsuitable where the transition operator has to be calculated for a wide spectrum of energies and many different isospin and initial momentum magnitude cases. Another approach uses the Arnoldi algorithm [5,9] to calculate the Krylov subspace ofB. Given a starting vector (scalar function) "t 1 ", it is spanned by The result of using the Arnoldi algorithm is a set of N orthogonal (with respect to a given scalar product) vectors that span (7) and a N × N matrix representation ofB. In practice N rarely has to be greater then 40, making the final equation easy to solve using standard linear solvers. The cut-off value for the momentum is 6 [fm −1 ], which is sufficient for the chiral NNLO potential [10]. In actual calculations the functions are scaled by a constant factor so that the normalization with respect to (11) is equal to 1.
A third approach that is the subject of this paper uses a basis composed from products of orthogonal polynomials to recreate the scalar functions "t" and to create a matrix representation ofB.

Orthogonal polynomials
The starting point for the new calculation scheme is a choice of basis for the scalar functions t i (|p |,p ·p) from (4). We consider two-dimensional functions T j composed from products of orthogonal polynomials X, P : where c j is a complex number, o(j) is the operator number (seew i=o(j) (p , p) from eq. (3)) and m(j), n(j) are polynomial numbers assigned to basis function j = 1 . . . N. We will be considering two choices for the polynomials P , X from [11] that are depicted on figs. 1 and 2. The cut-off value for the momentum is 6 [fm −1 ], which is sufficient for the chiral NNLO potential [10]. In actual calculations the functions are scaled by a constant factor so that the normalization with respect to (11) is equal to 1.
In order to create a matrix representation ofB from eq. (5) we act with this operator on all the basis functions (8). In the next step it is necessary to calculate the scalar product of the resulting functions with all the basis functions. There is a natural choice of scalar product for every type of, one-dimensional, orthogonal polynomial. This is outlined, for example, in chapts. 18.2 and 18.3 of [11]. If the scalar product of two polynomials P , P in the p = |p | direction with the cut-off value for the momentump is (P , P ) P = p 0 dpP (p)P (p)w P (p) (9) and the product of two polynomials X , X in the x =p ·p direction is then, for the two-dimensional functions Q i (|p |,p ·p), Q i (|p |,p ·p), we use: In these relations, w P (p) and w X (x) are appropriate weighing functions listed for example in [11]. Note the additional complex conjugation in (11) since we will be dealing with complex functions. The result of applyingB and calculating the scalar products is the possibility to rewrite eq. (5) in terms of N × N matrices and N -dimensional vectors: where B N ×N jk = (T j ,BT k ) and [t] N j = (T j , t). Typically in order to get agreement with results obtained using Arnoldi iterations and PWD we needed to use around 16 polynomials in each direction. This corresponds to N = 6 × 16 × 16 = 1536 basis functions and thus 1536, numerically demanding, applications ofB. The next section contains a comparison of the results obtained using orthogonal polynomials with the Krylev method that uses only around 40 applications ofB.

Numerical results
The transition operator was calculated for the on-shell case with the laboratory (projectile) energy 300 MeV using the method outlined in sect. 2. The two-dimensional functions from (4), (8) were discretized over a lattice of 32 × 72 points in the |p |,p ·p directions in order to perform all integrations numerically. The Wolfenstein parameters and scattering observables were calculated from the resulting transition operator scalar functions using the expressions from [1]. Finally the results were compared with calculations that used Arnoldi iterations [5,9]. The comparisons for the neutron-neutron and neutron-proton cases is shown in figs. 3, 4, 5 and 6. Convergence is observed to the Arnoldi results but only after 16 polynomials are used in each direction. The figure captions additionally contain the relative difference between the 16×16 polynomial observables A and the referential Arnoldi observables B calculated using where θ is the center-of-mass scattering angle. These figures also agree with fig. 1 from [3] where other methods were applied to come up with the solution. All calculations presented here used the NNLO [10] nucleon-nucleon potential. In order to obtain the on-shell results the transition operator scalar functions were calculated for two values: Selected observables for neutron-proton scattering at laboratory kinetic energy 300 MeV using the NNLO [10] nucleon-nucleon potential. σ is the differential cross-section with respect to the center-of-mass scattering angle θ. For the definition of A, D, R see, e.g. [1]. The solid line is obtained using the Arnoldi algorithm [5,9]. The dashed, dotted and dash-dotted lines are obtained using 8, 12 and 16 orthogonal polynomials from fig. 1 in each argument, respectively. The differences (calculated using (13)) between the 16 × 16 polynomial and Arnoldi results for R, A, D and the cross-section are 0.0058, 0.0094, 0.011 and 0.0059, respectively. Selected observables for neutron-proton scattering at laboratory kinetic energy 300 MeV using the NNLO [10] nucleonnucleon potential. σ is the differential cross-section with respect to the center-of-mass scattering angle θ. For the definition of A, D, R see e.g. [1]. The solid line is obtained using the Arnoldi algorithm [5,9]. The dashed, dotted and dash-dotted lines are obtained using 8, 12 and 16 orthogonal polynomials from fig. 2 in each argument, respectively. The differences (calculated using (13)) between the 16 × 16 polynomial and Arnoldi results for R, A, D and the cross-section are 0.0053, 0.012, 0.012 and 0.0054, respectively. It has been shown that the method described in this paper can be applied to calculations that involve the twonucleon transition operator. However, the large numerical cost of the calculations does not make it suitable for most practical calculations. The orthogonal polynomial approach requires approximately 16 × 16/40 = 6.4 times more computing resources. We suggest using Krylov subspace algorithms which are shown to be much more efficient. It also does not appear that the change of polynomial type will lead to better numerical performance. The main benefit of the Arnoldi algorithm approach is a careful choice of the set of orthogonal functions that will be used to create a matrix representation of the equations. This choice is made in such a way that these functions span the most numerically relevant subspace. It is highly unlikely that a new choice of orthogonal polynomial type will be a good alternative for this approach since the polynomials, contrary to the Krylov subpace basis functions, are not directly related to the physical problem. This is also illustrated in the observation that the relative differences between the polynomial result and the referential Arnoldi result are similar for both types of polynomial choices. The slow convergence of observables calculated using the standard partial wave techniques at higher energies for two [2] and three nucleon scattering [7] remains, however, a good motivation for the development of "threedimensional" approaches. Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.