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

Golak, Jacek

Skibiński, Roman

Topolnicki, Kacper

Witała, Henryk

The European Physical Journal A

English

Jagiellonian Univeristy Repository

DOI 10.1140/epja/i2016-16022-5Regular Article – Theoretical PhysicsEur. Phys. J. A (2016) 52: 22 THE EUROPEANPHYSICAL JOURNAL AOrthogonal polynomial approach to calculate the two-nucleontransition operator in three dimensionsKacper Topolnickia, Jacek Golak, Roman Skibiński, and Henryk WitalaM. Smoluchowski Institute of Physics, Jagiellonian University, PL-30348, Kraków, PolandReceived: 8 October 2015 / Revised: 11 January 2016Published online: 11 February 2016c© The Author(s) 2016. This article is published with open access at Springerlink.comCommunicated by S. HandsAbstract. We give a short report on the possibility to use orthogonal polynomials (OP) in calculationsthat involve the two-nucleon (2N) transition operator. The presented work adds another approach to theset 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 operatorcalculated at laboratory kinetic energy 300MeV. The new results for neutron-neutron and neutron-protonscattering observables converge to the results presented in Few-Body Syst. 53, 237 (2012) and to resultsobtained using the Arnoldi algorithm (Y. Saad, Iterative methods for sparse linear systems (SIAM Philadel-phia, PA, USA 2003)). The numerical cost of the calculations performed using the new scheme is large andthe new method can serve only as a backup to cross-check the previously used calculation schemes.1 IntroductionThe transition operator ť satisfying the Lippmann-Schwinger equation (LSE),ť(E) = ť(E) + V̌ Ǧ0(E + iε)ť(E), (1)contains all the necessary information to calculate any ob-servable in the nucleon-nucleon scattering process [1]. Itis also an important dynamical ingredient in momentumspace few-nucleon calculations [2–4]. This important ob-ject was a subject of our previous work [2,3,5] and in thispaper we apply an alternative approach to solve the LSEthat complements the previous methods. In eq. (1) V̌ isthe 2N potential, Ǧ0(E + iε) = (E + iε− p2m )−1 is the freepropagator, E is the energy, m is the nucleon mass and pis the relative momentum between the nucleons.The traditional approach to solve eq. (1) requires thepartial wave decomposition of all relevant operators. Thisstandard technique has been very successful in describingexperimental data at lower energies, however at higherenergies, sometimes, a very large number of partial wavesneeds to be taken into account in order to achieve conver-gence for certain observables [2,6,7]. In [2,3,5] we describeand utilize an alternative approach that uses directly thethree-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.pla general parity, time reversal and rotation invariant formof the 2N potential [8] is used:〈p′ | V̌ | p〉 =6∑i=1∑γvγi (|p′|, |p|, p̂′ · p̂) | γ〉〈γ | ⊗w̌i(p′,p), (2)with |γ〉 being one of the four possible isospin states of the2N system, | p〉 being a relative momentum eigenstate,w̌i(p′,p) are known 2N spin operators (they are listed forexample in [2]) and vγi (|p′|, |p|, p̂′ · p̂) are scalar functionsthat define the potential. The transition operator has asimilar operator form:〈p′ | ť(E) | p〉 =6∑i=1∑γtγi (E; |p′|, |p|, p̂′ · p̂) | γ〉〈γ | ⊗w̌i(p′,p), (3)but with energy-dependent scalar functions tγi (E; |p′|, |p|,p̂′ · p̂).After inserting (2) and (3) into the LSE (1) the spin de-pendencies can be removed. What remains is a set of equa-tions that involve the scalar functions vγi (|p′|, |p|, p̂′ · p̂)and tγi (E; |p′|, |p|, p̂′ · p̂). A further observation that theresulting equations constitute a linear system with oper-ators acting on the scalar functions makes it possible touse Krylev subspace algorithms (for example the ArnoldiPage 2 of 6 Eur. Phys. J. A (2016) 52: 22algorithm [9]) to calculate the solution. A careful consid-eration of the resulting relations leads to some additionalsimplifications. It can be easily verified that every totalisospin state γ, initial relative momentum magnitude |p|and energy E case is governed by an independent equationand can be treated separately. It is therefore possible, foreach independant case, to consider only two-dimensionalfunctions t:ti(|p′|, p̂′ · p̂) = tγi (E; |p′|, |p|, p̂′ · p̂), (4)and the resulting equation has a form that is very similarto the LSE, but this time with operators acting on thefunctions (vectors) t,t = v + B̌t. (5)The linear operator B̌ is defined for example in [5]. Fornegative energies E it transforms the scalar function t intoa new function t′ according to the following:(B̌t)k(|p′|, x′) ≡ t′k(|p′|, x′) =∫ +∞0d|p′′|∫ 1−1dx′′∫ 2π0dφ′′6∑j=16∑j′=1|p′′|2E − |p′′|2m + iεvj(|p′|, |p′′|,√1 − x′2√1 − x′′2 cos φ′′ + x′x′′)Bkjj′(|p′|, |p|, x′, |p′′|, x′′, φ′′)tj′(|p′′|, x′′), (6)where the momentum p′′ = |p′′|(√1−x′′2 cos φ′′,√1−x′′2sin φ′′, x′′) and the expressions for the functions Bkjj′(|p′|,|p|, x′, |p′′|, x′′, φ′′) are listed, e.g., in [5]. For nucleon-nucleon scattering and positive energies, eq. (6) has tobe reformulated to treat the singularity of the free propa-gator. This can be performed with the use of the standardmethods 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 ofthe transition operator (3) can be worked out using var-ious approaches. One approach requires the constructionof the matrix representation of the operator B̌ directlyfrom (6). This is a difficult task and the resulting matri-ces are large making this method unsuitable where thetransition operator has to be calculated for a wide spec-trum of energies and many different isospin and initialmomentum magnitude cases. Another approach uses theArnoldi algorithm [5,9] to calculate the Krylov subspaceof B̌. Given a starting vector (scalar function) “t1”, it isspanned byK(t1, N) = span(t1, B̌t1, . . . , B̌N t1). (7)The result of using the Arnoldi algorithm is a set of Northogonal (with respect to a given scalar product) vectorsthat span (7) and a N ×N matrix representation of B̌. Inpractice N rarely has to be greater then 40, making thefinal equation easy to solve using standard linear solvers.−1.0 −0.5 0.0 0.5 1.0−4−2024xP(α,β) n(x)0 1 2 3 4 5 6−4−2024p [1/fm]Un(2p/6−1)Fig. 1. Jacobi polynomials in the x = p̂′ ·p̂ argument (top) andshifted Chebyshev polynomials in the p = |p| argument (bot-tom) with α = 1, β = 1. The polynomial numbers are n = 0, 1,2, 3 (dashed, dotted, dot-dashed and solid line, respectively).The cut-off value for the momentum is 6 [fm−1], which is suffi-cient for the chiral NNLO potential [10]. In actual calculationsthe functions are scaled by a constant factor so that the nor-malization with respect to (11) is equal to 1.A third approach that is the subject of this paper uses abasis composed from products of orthogonal polynomialsto recreate the scalar functions “t” and to create a matrixrepresentation of B̌.2 Orthogonal polynomialsThe starting point for the new calculation scheme isa choice of basis for the scalar functions ti(|p′|, p̂′ · p̂)from (4). We consider two-dimensional functions T j com-posed from products of orthogonal polynomials X, P :ti(|p′|, p̂′ · p̂) =N∑j=1cjT j(|p′|, p̂′ · p̂)T ji (|p′|, p̂′ · p̂) = δio(j)Pn(j)(|p′|)Xm(j)(p̂′ · p̂)), (8)where cj is a complex number, o(j) is the operator number(see w̌i=o(j)(p′,p) from eq. (3)) and m(j), n(j) are poly-nomial numbers assigned to basis function j = 1 . . . N . Wewill be considering two choices for the polynomials P , Xfrom [11] that are depicted on figs. 1 and 2.Eur. Phys. J. A (2016) 52: 22 Page 3 of 6−1.0 −0.5 0.0 0.5 1.0−4−2024xUn(x)0 1 2 3 4 5 6−6−4−2024p [1/ fm]P(α,β) n(2p/6−1)Fig. 2. Chebyshev polynomials in the x = p̂′ ·p̂ argument (top)and shifted Jacobi polynomials in the p = |p| argument (bot-tom) with α = 0, β = 1.5. The polynomial numbers are n = 0,1, 2, 3 (dashed, dotted, dot-dashed and solid line, respectively).The cut-off value for the momentum is 6 [fm−1], which is suffi-cient for the chiral NNLO potential [10]. In actual calculationsthe functions are scaled by a constant factor so that the nor-malization with respect to (11) is equal to 1.In order to create a matrix representation of B̌ fromeq. (5) we act with this operator on all the basis func-tions (8). In the next step it is necessary to calculate thescalar product of the resulting functions with all the ba-sis functions. There is a natural choice of scalar productfor every type of, one-dimensional, orthogonal polynomial.This is outlined, for example, in chapts. 18.2 and 18.3of [11]. If the scalar product of two polynomials P ′, P inthe p = |p′| direction with the cut-off value for the mo-mentum p̄ is(P ′, P )P =∫ p̄0dpP ′(p)P (p)wP (p) (9)and the product of two polynomials X ′, X in the x = p̂′ ·p̂direction is(X ′,X)X =∫ 1−1dxX ′(x)X(x)wX(x), (10)then, for the two-dimensional functions Q′i(|p′|, p̂′ · p̂),Qi(|p′|, p̂′ · p̂), we use:(Q′, Q) =6∑i=1∫ p̄0dp∫ 1−1dx(Q′i(p, x))∗Qi(p, x)wP (p)wX(x). (11)In these relations, wP (p) and wX(x) are appropriateweighing functions listed for example in [11]. Note theadditional complex conjugation in (11) since we will bedealing with complex functions.The result of applying B̌ and calculating the scalarproducts is the possibility to rewrite eq. (5) in terms ofN × N matrices and N -dimensional vectors:[t]N = [v]N +[B̌]N×N[t]N . (12)where[B̌]N×Njk= (T j , B̌T k) and [t]Nj = (T j , t). Typi-cally in order to get agreement with results obtained us-ing Arnoldi iterations and PWD we needed to use around16 polynomials in each direction. This corresponds toN = 6 × 16 × 16 = 1536 basis functions and thus 1536,numerically demanding, applications of B̌. The next sec-tion contains a comparison of the results obtained usingorthogonal polynomials with the Krylev method that usesonly around 40 applications of B̌.3 Numerical resultsThe transition operator was calculated for the on-shellcase with the laboratory (projectile) energy 300 MeV us-ing the method outlined in sect. 2. The two-dimensionalfunctions from (4), (8) were discretized over a lattice of32× 72 points in the |p′|, p̂′ · p̂ directions in order to per-form all integrations numerically. The Wolfenstein param-eters and scattering observables were calculated from theresulting transition operator scalar functions using the ex-pressions from [1]. Finally the results were compared withcalculations that used Arnoldi iterations [5,9]. The com-parisons for the neutron-neutron and neutron-proton casesis shown in figs. 3, 4, 5 and 6. Convergence is observed tothe Arnoldi results but only after 16 polynomials are usedin each direction. The figure captions additionally con-tain the relative difference between the 16×16 polynomialobservables A and the referential Arnoldi observables Bcalculated using〈A − B〉 =∫ 180◦0dθ|A(θ) − B(θ)|∫ 180◦0dθ|B(θ)|, (13)where θ is the center-of-mass scattering angle. These fig-ures also agree with fig. 1 from [3] where other methodswere applied to come up with the solution.All calculations presented here used the NNLO [10]nucleon-nucleon potential. In order to obtain the on-shellresults the transition operator scalar functions were cal-culated for two values:tγi (E; p0, p0 + δ, p̂′ · p̂)andtγi (E; p0, p0 − δ, p̂′ · p̂)with p0 =√2mE. The on-shell value is the average of thetwo with δ = 0.01[fm−1].Page 4 of 6 Eur. Phys. J. A (2016) 52: 220 50 100 150024681012θ [deg]σ[mb/sr]0 50 100 150−0.4−0.20.00.20.40.6θ [deg]A0 50 100 150−0.4−0.20.00.20.40.60.8θ [deg]D0 50 100 150−0.20.00.20.40.6θ [deg]RFig. 3. 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 thedefinition of A, D, R see, e.g. [1]. The solid line is obtained using the Arnoldi algorithm [5,9]. The dashed, dotted anddash-dotted lines are obtained using 8, 12 and 16 orthogonal polynomials from fig. 1 in each argument, respectively. Thedifferences (calculated using (13)) between the 16 × 16 polynomial and Arnoldi results for R, A, D and the cross-section are0.0058, 0.0094, 0.011 and 0.0059, respectively.0 50 100 150024681012θ [deg]σ[mb/sr]0 50 100 150−0.4−0.20.00.20.40.6θ [deg]A0 50 100 150−0.4−0.20.00.20.40.60.8θ [deg]D0 50 100 150−0.20.00.20.40.6θ [deg]RFig. 4. 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 ofA, D, R see e.g. [1]. The solid line is obtained using the Arnoldi algorithm [5,9]. The dashed, dotted and dash-dotted linesare obtained using 8, 12 and 16 orthogonal polynomials from fig. 2 in each argument, respectively. The differences (calculatedusing (13)) between the 16 × 16 polynomial and Arnoldi results for R, A, D and the cross-section are 0.0053, 0.012, 0.012 and0.0054, respectively.Eur. Phys. J. A (2016) 52: 22 Page 5 of 60 50 100 15002468101214θ [deg]σ[mb/sr]0 50 100 150−0.6−0.4−0.20.00.2θ [deg]A0 50 100 1500.00.20.40.6θ [deg]D0 50 100 150−0.2−0.10.00.10.20.30.4θ [deg]RFig. 5. Same as in fig. 3 but for selected neutron-neutron scattering observables. The differences (calculated using (13)) betweenthe 16 × 16 polynomial and Arnoldi results for R, A, D and the cross-section are 0.016, 0.0074, 0.0073 and 0.0045, respectively.0 50 100 15002468101214θ [deg]σ[mb/sr]0 50 100 150−0.6−0.4−0.20.00.2θ [deg]A0 50 100 1500.00.20.40.60.8θ [deg]D0 50 100 150−0.2−0.10.00.10.20.30.4θ [deg]RFig. 6. Same as in fig. 4 but for selected neutron-neutron scattering observables. The differences (calculated using (13)) betweenthe 16×16 polynomial and Arnoldi results for R, A, D and the cross-section are 0.011682, 0.0073, 0.0076 and 0.0057, respectively.Page 6 of 6 Eur. Phys. J. A (2016) 52: 22It has been shown that the method described in thispaper can be applied to calculations that involve the two-nucleon transition operator. However, the large numeri-cal cost of the calculations does not make it suitable formost practical calculations. The orthogonal polynomialapproach requires approximately 16 × 16/40 = 6.4 timesmore computing resources. We suggest using Krylov sub-space algorithms which are shown to be much more effi-cient. It also does not appear that the change of polyno-mial type will lead to better numerical performance. Themain benefit of the Arnoldi algorithm approach is a care-ful choice of the set of orthogonal functions that will beused to create a matrix representation of the equations.This choice is made in such a way that these functionsspan the most numerically relevant subspace. It is highlyunlikely that a new choice of orthogonal polynomial typewill be a good alternative for this approach since the poly-nomials, contrary to the Krylov subpace basis functions,are not directly related to the physical problem. This isalso illustrated in the observation that the relative differ-ences between the polynomial result and the referentialArnoldi result are similar for both types of polynomialchoices. The slow convergence of observables calculatedusing the standard partial wave techniques at higher en-ergies for two [2] and three nucleon scattering [7] remains,however, a good motivation for the development of “three-dimensional” approaches.The project was financed from the resources of the Na-tional Science Center (Poland) under grants No. DEC-2013/11/N/ST2/03733 and DEC-2013/10/ M/ST2/00420.Some numerical calculations have been performed on the su-percomputer clusters of the JSC, Jülich, Germany.Open Access This is an open access article distributedunder the terms of the Creative Commons AttributionLicense (http://creativecommons.org/licenses/by/4.0), whichpermits unrestricted use, distribution, and reproduction in anymedium, provided the original work is properly cited.References1. W. Glöckle, The Quantum Mechanical Few-Body Problem(Springer-Verlag, Berlin-Heidelberg, 1983).2. J. Golak, W. Glöckle, R. Skibiński, H. Witala, D.Rozpedzik, K. Topolnicki, I. Fachruddin, Ch. Elster, A.Nogga, Phys. Rev. C 81, 034006 (2010).3. J. Golak, R. Skibiński, H. Witala, K. Topolnicki, W.Glöckle, A. Nogga, H. Kamada, Few-Body Syst. 53, 237(2012).4. J. Golak, R. Skibiński, H. Witala, K. Topolnicki, A.E.Elmeshneb, H. Kamada, A. Nogga, L. Marcucci, Phys.Rev. C 90, 024001 (2014).5. K. Topolnicki, PhD thesis, Jagiellonian University (2014)unpublished, available online at the following address:http://www.fais.uj.edu.pl/dla-studentow/studia-doktoranckie/prace-doktorskie#2014.6. H. Witala, J. Golak, R. Skibiński, K. Topolnicki, J. Phys.G: Nucl. Part. Phys. 41, 094011 (2014).7. K. Topolnicki, J. Golak, R. Skibiński, H. Witala, C.A.Bertulani, Eur. Phys. J. A 51, 132 (2015).8. L. Wolfenstein, Phys. Rev. 96, 1654 (1954).9. Y. Saad, Iterative methods for sparse linear systems (SIAMPhiladelphia, PA, USA 2003).10. E. Epelbaum, Prog. Part. Nucl. Phys. 57, 654 (2006).11. NIST Digital Library of Mathematical Functions,http://dlmf.nist.gov, Release 1.0.9 of 2014-08-29,online companion to OLBC10.

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

Kacper Topolnicki

Jacek Golak

Roman Skibiński

Henryk Witała

Springer - Publisher Connector

Crossref

DOI 10.1140/epja/i2016-16022-5Regular Article – Theoretical PhysicsEur. Phys. J. A (2016) 52: 22 THE EUROPEANPHYSICAL JOURNAL AOrthogonal polynomial approach to calculate the two-nucleontransition operator in three dimensionsKacper Topolnickia, Jacek Golak, Roman Skibin´ski, and Henryk WitalaM. Smoluchowski Institute of Physics, Jagiellonian University, PL-30348, Krako´w, PolandReceived: 8 October 2015 / Revised: 11 January 2016Published online: 11 February 2016c© The Author(s) 2016. This article is published with open access at Springerlink.comCommunicated by S. HandsAbstract. We give a short report on the possibility to use orthogonal polynomials (OP) in calculationsthat involve the two-nucleon (2N) transition operator. The presented work adds another approach to theset 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 operatorcalculated at laboratory kinetic energy 300MeV. The new results for neutron-neutron and neutron-protonscattering observables converge to the results presented in Few-Body Syst. 53, 237 (2012) and to resultsobtained using the Arnoldi algorithm (Y. Saad, Iterative methods for sparse linear systems (SIAM Philadel-phia, PA, USA 2003)). The numerical cost of the calculations performed using the new scheme is large andthe new method can serve only as a backup to cross-check the previously used calculation schemes.1 IntroductionThe transition operator tˇ satisfying the Lippmann-Schwinger equation (LSE),tˇ(E) = tˇ(E) + Vˇ Gˇ0(E + i)tˇ(E), (1)contains all the necessary information to calculate any ob-servable in the nucleon-nucleon scattering process [1]. Itis also an important dynamical ingredient in momentumspace few-nucleon calculations [2–4]. This important ob-ject was a subject of our previous work [2,3,5] and in thispaper we apply an alternative approach to solve the LSEthat complements the previous methods. In eq. (1) Vˇ isthe 2N potential, Gˇ0(E + i) = (E + i− p2m )−1 is the freepropagator, E is the energy, m is the nucleon mass and pis the relative momentum between the nucleons.The traditional approach to solve eq. (1) requires thepartial wave decomposition of all relevant operators. Thisstandard technique has been very successful in describingexperimental data at lower energies, however at higherenergies, sometimes, a very large number of partial wavesneeds to be taken into account in order to achieve conver-gence for certain observables [2,6,7]. In [2,3,5] we describeand utilize an alternative approach that uses directly thethree-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.pla general parity, time reversal and rotation invariant formof the 2N potential [8] is used:〈p′ | Vˇ | p〉 =6∑i=1∑γvγi (|p′|, |p|, pˆ′ · pˆ) | γ〉〈γ | ⊗wˇi(p′,p), (2)with |γ〉 being one of the four possible isospin states of the2N system, | p〉 being a relative momentum eigenstate,wˇi(p′,p) are known 2N spin operators (they are listed forexample in [2]) and vγi (|p′|, |p|, pˆ′ · pˆ) are scalar functionsthat deﬁne the potential. The transition operator has asimilar operator form:〈p′ | tˇ(E) | p〉 =6∑i=1∑γtγi (E; |p′|, |p|, pˆ′ · pˆ) | γ〉〈γ | ⊗wˇi(p′,p), (3)but with energy-dependent scalar functions tγi (E; |p′|, |p|,pˆ′ · pˆ).After inserting (2) and (3) into the LSE (1) the spin de-pendencies can be removed. What remains is a set of equa-tions that involve the scalar functions vγi (|p′|, |p|, pˆ′ · pˆ)and tγi (E; |p′|, |p|, pˆ′ · pˆ). A further observation that theresulting equations constitute a linear system with oper-ators acting on the scalar functions makes it possible touse Krylev subspace algorithms (for example the ArnoldiPage 2 of 6 Eur. Phys. J. A (2016) 52: 22algorithm [9]) to calculate the solution. A careful consid-eration of the resulting relations leads to some additionalsimpliﬁcations. It can be easily veriﬁed that every totalisospin state γ, initial relative momentum magnitude |p|and energy E case is governed by an independent equationand can be treated separately. It is therefore possible, foreach independant case, to consider only two-dimensionalfunctions t:ti(|p′|, pˆ′ · pˆ) = tγi (E; |p′|, |p|, pˆ′ · pˆ), (4)and the resulting equation has a form that is very similarto the LSE, but this time with operators acting on thefunctions (vectors) t,t = v + Bˇt. (5)The linear operator Bˇ is deﬁned for example in [5]. Fornegative energies E it transforms the scalar function t intoa new function t′ according to the following:(Bˇt)k(|p′|, x′) ≡ t′k(|p′|, x′) =∫ +∞0d|p′′|∫ 1−1dx′′∫ 2π0dφ′′6∑j=16∑j′=1|p′′|2E − |p′′|2m + ivj(|p′|, |p′′|,√1− x′2√1− x′′2 cosφ′′ + x′x′′)Bkjj′(|p′|, |p|, x′, |p′′|, x′′, φ′′)tj′(|p′′|, x′′), (6)where the momentum p′′ = |p′′|(√1−x′′2 cosφ′′,√1−x′′2sinφ′′, x′′) and the expressions for the functions Bkjj′(|p′|,|p|, x′, |p′′|, x′′, φ′′) are listed, e.g., in [5]. For nucleon-nucleon scattering and positive energies, eq. (6) has tobe reformulated to treat the singularity of the free propa-gator. This can be performed with the use of the standardmethods 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 ofthe transition operator (3) can be worked out using var-ious approaches. One approach requires the constructionof the matrix representation of the operator Bˇ directlyfrom (6). This is a diﬃcult task and the resulting matri-ces are large making this method unsuitable where thetransition operator has to be calculated for a wide spec-trum of energies and many diﬀerent isospin and initialmomentum magnitude cases. Another approach uses theArnoldi algorithm [5,9] to calculate the Krylov subspaceof Bˇ. Given a starting vector (scalar function) “t1”, it isspanned byK(t1, N) = span(t1, Bˇt1, . . . , BˇN t1). (7)The result of using the Arnoldi algorithm is a set of Northogonal (with respect to a given scalar product) vectorsthat span (7) and a N ×N matrix representation of Bˇ. Inpractice N rarely has to be greater then 40, making theﬁnal equation easy to solve using standard linear solvers.−1.0 −0.5 0.0 0.5 1.0−4−2024xP(α,β) n(x)0 1 2 3 4 5 6−4−2024p [1/fm]Un(2p/6−1)Fig. 1. Jacobi polynomials in the x = pˆ′ ·pˆ argument (top) andshifted Chebyshev polynomials in the p = |p| argument (bot-tom) with α = 1, β = 1. The polynomial numbers are n = 0, 1,2, 3 (dashed, dotted, dot-dashed and solid line, respectively).The cut-oﬀ value for the momentum is 6 [fm−1], which is suﬃ-cient for the chiral NNLO potential [10]. In actual calculationsthe functions are scaled by a constant factor so that the nor-malization with respect to (11) is equal to 1.A third approach that is the subject of this paper uses abasis composed from products of orthogonal polynomialsto recreate the scalar functions “t” and to create a matrixrepresentation of Bˇ.2 Orthogonal polynomialsThe starting point for the new calculation scheme isa choice of basis for the scalar functions ti(|p′|, pˆ′ · pˆ)from (4). We consider two-dimensional functions T j com-posed from products of orthogonal polynomials X, P :ti(|p′|, pˆ′ · pˆ) =N∑j=1cjT j(|p′|, pˆ′ · pˆ)T ji (|p′|, pˆ′ · pˆ) = δio(j)Pn(j)(|p′|)Xm(j)(pˆ′ · pˆ)), (8)where cj is a complex number, o(j) is the operator number(see wˇi=o(j)(p′,p) from eq. (3)) and m(j), n(j) are poly-nomial numbers assigned to basis function j = 1 . . . N . Wewill be considering two choices for the polynomials P , Xfrom [11] that are depicted on ﬁgs. 1 and 2.Eur. Phys. J. A (2016) 52: 22 Page 3 of 6−1.0 −0.5 0.0 0.5 1.0−4−2024xUn(x)0 1 2 3 4 5 6−6−4−2024p [1/ fm]P(α,β) n(2p/6−1)Fig. 2. Chebyshev polynomials in the x = pˆ′ ·pˆ argument (top)and shifted Jacobi polynomials in the p = |p| argument (bot-tom) with α = 0, β = 1.5. The polynomial numbers are n = 0,1, 2, 3 (dashed, dotted, dot-dashed and solid line, respectively).The cut-oﬀ value for the momentum is 6 [fm−1], which is suﬃ-cient for the chiral NNLO potential [10]. In actual calculationsthe functions are scaled by a constant factor so that the nor-malization with respect to (11) is equal to 1.In order to create a matrix representation of Bˇ fromeq. (5) we act with this operator on all the basis func-tions (8). In the next step it is necessary to calculate thescalar product of the resulting functions with all the ba-sis functions. There is a natural choice of scalar productfor every type of, one-dimensional, orthogonal polynomial.This is outlined, for example, in chapts. 18.2 and 18.3of [11]. If the scalar product of two polynomials P ′, P inthe p = |p′| direction with the cut-oﬀ value for the mo-mentum p¯ is(P ′, P )P =∫ p¯0dpP ′(p)P (p)wP (p) (9)and the product of two polynomials X ′, X in the x = pˆ′ ·pˆdirection is(X ′,X)X =∫ 1−1dxX ′(x)X(x)wX(x), (10)then, for the two-dimensional functions Q′i(|p′|, pˆ′ · pˆ),Qi(|p′|, pˆ′ · pˆ), we use:(Q′, Q) =6∑i=1∫ p¯0dp∫ 1−1dx(Q′i(p, x))∗Qi(p, x)wP (p)wX(x). (11)In these relations, wP (p) and wX(x) are appropriateweighing functions listed for example in [11]. Note theadditional complex conjugation in (11) since we will bedealing with complex functions.The result of applying Bˇ and calculating the scalarproducts is the possibility to rewrite eq. (5) in terms ofN ×N matrices and N -dimensional vectors:[t]N = [v]N +[Bˇ]N×N [t]N . (12)where[Bˇ]N×Njk= (T j , BˇT k) and [t]Nj = (T j , t). Typi-cally in order to get agreement with results obtained us-ing Arnoldi iterations and PWD we needed to use around16 polynomials in each direction. This corresponds toN = 6 × 16 × 16 = 1536 basis functions and thus 1536,numerically demanding, applications of Bˇ. The next sec-tion contains a comparison of the results obtained usingorthogonal polynomials with the Krylev method that usesonly around 40 applications of Bˇ.3 Numerical resultsThe transition operator was calculated for the on-shellcase with the laboratory (projectile) energy 300 MeV us-ing the method outlined in sect. 2. The two-dimensionalfunctions from (4), (8) were discretized over a lattice of32× 72 points in the |p′|, pˆ′ · pˆ directions in order to per-form all integrations numerically. The Wolfenstein param-eters and scattering observables were calculated from theresulting transition operator scalar functions using the ex-pressions from [1]. Finally the results were compared withcalculations that used Arnoldi iterations [5,9]. The com-parisons for the neutron-neutron and neutron-proton casesis shown in ﬁgs. 3, 4, 5 and 6. Convergence is observed tothe Arnoldi results but only after 16 polynomials are usedin each direction. The ﬁgure captions additionally con-tain the relative diﬀerence between the 16×16 polynomialobservables A and the referential Arnoldi observables Bcalculated using〈A−B〉 =∫ 180◦0dθ|A(θ)−B(θ)|∫ 180◦0dθ|B(θ)|, (13)where θ is the center-of-mass scattering angle. These ﬁg-ures also agree with ﬁg. 1 from [3] where other methodswere applied to come up with the solution.All calculations presented here used the NNLO [10]nucleon-nucleon potential. In order to obtain the on-shellresults the transition operator scalar functions were cal-culated for two values:tγi (E; p0, p0 + δ, pˆ′ · pˆ)andtγi (E; p0, p0 − δ, pˆ′ · pˆ)with p0 =√2mE. The on-shell value is the average of thetwo with δ = 0.01[fm−1].Page 4 of 6 Eur. Phys. J. A (2016) 52: 220 50 100 150024681012θ [deg]σ[mb/sr]0 50 100 150−0.4−0.20.00.20.40.6θ [deg]A0 50 100 150−0.4−0.20.00.20.40.60.8θ [deg]D0 50 100 150−0.20.00.20.40.6θ [deg]RFig. 3. Selected observables for neutron-proton scattering at laboratory kinetic energy 300MeV using the NNLO [10]nucleon-nucleon potential. σ is the diﬀerential cross-section with respect to the center-of-mass scattering angle θ. For thedeﬁnition of A, D, R see, e.g. [1]. The solid line is obtained using the Arnoldi algorithm [5,9]. The dashed, dotted anddash-dotted lines are obtained using 8, 12 and 16 orthogonal polynomials from ﬁg. 1 in each argument, respectively. Thediﬀerences (calculated using (13)) between the 16 × 16 polynomial and Arnoldi results for R, A, D and the cross-section are0.0058, 0.0094, 0.011 and 0.0059, respectively.0 50 100 150024681012θ [deg]σ[mb/sr]0 50 100 150−0.4−0.20.00.20.40.6θ [deg]A0 50 100 150−0.4−0.20.00.20.40.60.8θ [deg]D0 50 100 150−0.20.00.20.40.6θ [deg]RFig. 4. Selected observables for neutron-proton scattering at laboratory kinetic energy 300MeV using the NNLO [10] nucleon-nucleon potential. σ is the diﬀerential cross-section with respect to the center-of-mass scattering angle θ. For the deﬁnition ofA, D, R see e.g. [1]. The solid line is obtained using the Arnoldi algorithm [5,9]. The dashed, dotted and dash-dotted linesare obtained using 8, 12 and 16 orthogonal polynomials from ﬁg. 2 in each argument, respectively. The diﬀerences (calculatedusing (13)) between the 16× 16 polynomial and Arnoldi results for R, A, D and the cross-section are 0.0053, 0.012, 0.012 and0.0054, respectively.Eur. Phys. J. A (2016) 52: 22 Page 5 of 60 50 100 15002468101214θ [deg]σ[mb/sr]0 50 100 150−0.6−0.4−0.20.00.2θ [deg]A0 50 100 1500.00.20.40.6θ [deg]D0 50 100 150−0.2−0.10.00.10.20.30.4θ [deg]RFig. 5. Same as in ﬁg. 3 but for selected neutron-neutron scattering observables. The diﬀerences (calculated using (13)) betweenthe 16× 16 polynomial and Arnoldi results for R, A, D and the cross-section are 0.016, 0.0074, 0.0073 and 0.0045, respectively.0 50 100 15002468101214θ [deg]σ[mb/sr]0 50 100 150−0.6−0.4−0.20.00.2θ [deg]A0 50 100 1500.00.20.40.60.8θ [deg]D0 50 100 150−0.2−0.10.00.10.20.30.4θ [deg]RFig. 6. Same as in ﬁg. 4 but for selected neutron-neutron scattering observables. The diﬀerences (calculated using (13)) betweenthe 16×16 polynomial and Arnoldi results for R, A, D and the cross-section are 0.011682, 0.0073, 0.0076 and 0.0057, respectively.Page 6 of 6 Eur. Phys. J. A (2016) 52: 22It has been shown that the method described in thispaper can be applied to calculations that involve the two-nucleon transition operator. However, the large numeri-cal cost of the calculations does not make it suitable formost practical calculations. The orthogonal polynomialapproach requires approximately 16 × 16/40 = 6.4 timesmore computing resources. We suggest using Krylov sub-space algorithms which are shown to be much more eﬃ-cient. It also does not appear that the change of polyno-mial type will lead to better numerical performance. Themain beneﬁt of the Arnoldi algorithm approach is a care-ful choice of the set of orthogonal functions that will beused to create a matrix representation of the equations.This choice is made in such a way that these functionsspan the most numerically relevant subspace. It is highlyunlikely that a new choice of orthogonal polynomial typewill be a good alternative for this approach since the poly-nomials, contrary to the Krylov subpace basis functions,are not directly related to the physical problem. This isalso illustrated in the observation that the relative diﬀer-ences between the polynomial result and the referentialArnoldi result are similar for both types of polynomialchoices. The slow convergence of observables calculatedusing the standard partial wave techniques at higher en-ergies for two [2] and three nucleon scattering [7] remains,however, a good motivation for the development of “three-dimensional” approaches.The project was ﬁnanced from the resources of the Na-tional Science Center (Poland) under grants No. DEC-2013/11/N/ST2/03733 and DEC-2013/10/ M/ST2/00420.Some numerical calculations have been performed on the su-percomputer clusters of the JSC, Ju¨lich, Germany.Open Access This is an open access article distributedunder the terms of the Creative Commons AttributionLicense (http://creativecommons.org/licenses/by/4.0), whichpermits unrestricted use, distribution, and reproduction in anymedium, provided the original work is properly cited.References1. W. Glo¨ckle, The Quantum Mechanical Few-Body Problem(Springer-Verlag, Berlin-Heidelberg, 1983).2. J. Golak, W. Glo¨ckle, R. Skibin´ski, H. Witala, D.Rozpedzik, K. Topolnicki, I. Fachruddin, Ch. Elster, A.Nogga, Phys. Rev. C 81, 034006 (2010).3. J. Golak, R. Skibin´ski, H. Witala, K. Topolnicki, W.Glo¨ckle, A. Nogga, H. Kamada, Few-Body Syst. 53, 237(2012).4. J. Golak, R. Skibin´ski, H. Witala, K. Topolnicki, A.E.Elmeshneb, H. Kamada, A. Nogga, L. Marcucci, Phys.Rev. C 90, 024001 (2014).5. K. Topolnicki, PhD thesis, Jagiellonian University (2014)unpublished, available online at the following address:http://www.fais.uj.edu.pl/dla-studentow/studia-doktoranckie/prace-doktorskie#2014.6. H. Witala, J. Golak, R. Skibin´ski, K. Topolnicki, J. Phys.G: Nucl. Part. Phys. 41, 094011 (2014).7. K. Topolnicki, J. Golak, R. Skibin´ski, H. Witala, C.A.Bertulani, Eur. Phys. J. A 51, 132 (2015).8. L. Wolfenstein, Phys. Rev. 96, 1654 (1954).9. Y. Saad, Iterative methods for sparse linear systems (SIAMPhiladelphia, PA, USA 2003).10. E. Epelbaum, Prog. Part. Nucl. Phys. 57, 654 (2006).11. NIST Digital Library of Mathematical Functions,http://dlmf.nist.gov, Release 1.0.9 of 2014-08-29,online companion to OLBC10.

https://link.springer.com/content/pdf/10.1140%2Fepja%2Fi2016-16022-5.pdf

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

Abstract

Similar works

Full text

Available Versions

Jagiellonian Univeristy Repository

Springer - Publisher Connector

Crossref