A momentum-space Argonne V18 interaction

S. Veerasamy and W. N. Polyzou Department of Physics and Astronomy, The University of Iowa, Iowa City, IA 52242

(Dated: June 7, 2011)

Abstract

This paper gives a momentum-space representation of the Argonne V18 potential as an expansion in products of spin-isospin operators with scalar coefficient functions of the momentum transfer. Two representations of the scalar coefficient functions for the strong part of the interaction are given. One is as an expansion in an orthonormal basis of rational functions and the other as an expansion in Chebyshev polynomials on different intervals. Both provide practical and efficient representations for computing the momentum-space potential that do not require integration or interpolation. Programs based on both expansions are available as supplementary material. Analytic expressions are given for the scalar coefficient functions of the Fourier transform of the electromagnetic part of the Argonne V18. A simple method for computing the partial-wave projections of these interactions from the operator expressions is also given.

PACS numbers: 21.45.Bc ,21.30.Cb

I. INTRODUCTION

The Argonne V18 potential [1] is one of a number of nucleon-nucleon interactions [2][1][3] that provide a quantitative description of experimental two-body observables below the pionproduction threshold. It is distinguished from the other realistic interactions because it is expressed as an operator expansion with local configuration-space coefficient functions. This representation has advantages when used in variational Monte Carlo calculations. On the other hand, there are a number of calculations that require a realistic interaction that are more naturally performed in momentum space. These include some Faddeev calculations, relativistic few-body calculations, and calculations involving electromagnetic probes. In the momentum representation the variable conjugate to the relative coordinate is the momentum transfer. In calculations, both momenta appear, which requires either an interpolation or a separate Fourier transform for each pair of momenta. Fourier transforms of the V18 potential have been used in some applications [4]. The purpose of this paper is to provide useful, tested and reproducible analytic approximations of the Fourier transform of the Argonne V18 potential for use in momentum-space calculations. The analytic forms allow for a direct calculation of the momentum-space interaction for any pair of initial and final momenta. In keeping with the traditional Argonne form, the momentum-space potential is given as a linear combination of products of spin-isospin operators with scalar functions of the momentum transfer. The resulting momentum-space potential has 24 terms. The additional six operators appear because the Fourier transform of the terms involving the operators L^{2}V_{i}(r) and (L · S)^{2}V_{i}(r) each become a sum of two different momentum-space operators with different coefficient functions. In this work the Fourier transform is given for the strong part of the Argonne V18 potential, without the electromagnetic terms. This part of the potential must be treated numerically. The electromagnetic terms have analytic Fourier transforms, which are discussed in Appendix 3. The partial-wave projection of the momentum space potential is discussed in Appendix 2. It is constructed from the operator expressions by integrating over the angle between the initial and final momentum vectors, however unlike the configuration-space partial-wave projection, the integrals involve both the operator and the scalar coefficient functions.

The Argonne V18 potential has the form

V = ∑_{n=1}^{18} V_{n}(r)O_{n}  (1.1)

where V_{n}(r) are rotationally-invariant coefficient functions of the relative coordinate of the nucleons and the O_{n} are the eighteen spin-isospin operators given in Table 1.,

Table 1: Argonne V18 spin-isospin operators

in coordinate-space

Term spin-isospin Operator in r-space
O_{1}  I
O_{2}  (τ_{1} · τ_{2})
O_{3}  (σ σ σ_{1} · σ σ σ_{2}),
O_{4}  (σ σ σ_{1} · σ σ σ_{2})(τ_{1} · τ_{2})
O_{5} S_{12} = 3(σ σ σ_{1} · r̂)(σ σ σ_{2} · r̂) − σ σ σ_{1} · σ σ σ_{2}
O_{6}  S_{12}(τ_{1} · τ_{2}),
O_{7}  (L · S)
O_{8}  (L · S)(τ_{1} · τ_{2})
O_{9}  (L · L)
O_{10}  (L · L)(τ_{1} · τ_{2})
O_{11}  (L · L)(σ σ σ_{1} · σ σ σ_{2})
O_{12}  (L · L)(σ σ σ_{1} · σ σ σ_{2})(τ_{1} · τ_{2})
O_{13}  (L · S)^{2}
O_{14}  (L · S)^{2}(τ_{1} · τ_{2})
O_{15}  T_{12} = (3τ_{1z}τ_{2z} − τ · τ)
O_{16}  (σ σ σ_{1} · σ σ σ_{2})T_{12}
O_{17}  S_{12}T_{12}
O_{18}  (τ_{1z} + τ_{2z})

In this table T_{12} is the isotensor operator T_{12} := 3τ_{1z}τ_{2z} − τττ_{1} · τττ_{2}. While the isospin operators, τττ_{i}, factor out of the Fourier transforms, the operators L^{2}, L · S, (L · S)^{2} and the tensor operator S_{12} contribute to the Fourier transform.

The Fourier transform of this potential can be expressed as a linear combination of 24 momentum-space operators with scalar coefficient functions of the momentum transfer.

There are 24 operators because the L·L and (L·S)^{2} operators have two distinct contributions in momentum space. In appendix 1 it is shown that the potential matrix element ⟨k′|V |k⟩, with q := k′ − k, has the following five types of contributions:

1. I

\frac{1}{(2π)^{3}} ∫ e^{−i(k′−k)·r}V_{j}(r)Idr = I\frac{1}{2π^{2}} ∫_{0}^{∞} j_{0}(qr)V_{j}(r)r^{2}dr.  (1.2)

2. L · S

\frac{1}{(2π)^{3}} ∫ e^{−ik′·r}V_{j}(r)L · Se^{ik·r}dr = i(k × k′) · S\frac{1}{2π^{2}q} ∫_{0}^{∞} j_{1}(qr)V_{j}(r)r^{3}dr.  (1.3)

3. L · L

\frac{1}{(2π)^{3}} ∫ e^{−ik′·r}V_{j}(r)L · Le^{ik·r}dr =

−(k′×k)·(k′×k)\frac{1}{2π^{2}q^{2}} ∫_{0}^{∞} j_{2}(qr)V_{j}(r)r^{4}dr+2(k′·k)\frac{1}{2π^{2}q} ∫_{0}^{∞} j_{1}(qr)V_{j}(r)r^{3}dr. (1.4)

4. (L · S)^{2}

\frac{1}{(2π)^{3}} ∫ e^{−ik′·r}V_{j}(r)(L · S)^{2}e^{ik·r}dr =

−(S · (k × k′))^{2} \frac{1}{2π^{2}q^{2}} ∫_{0}^{∞} j_{2}(qr)V_{j}(r)r^{4}dr + (k′ × S) · (k × S)\frac{1}{2π^{2}q} ∫_{0}^{∞} j_{1}(qr)V_{j}(r)r^{3}dr.

(1.5)

5. S_{12} = 3(r̂ · σ σ σ_{1})(r̂ · σ σ σ_{2}) − σ σ σ_{1} · σ σ σ_{2}

\frac{1}{(2π)^{3}} ∫ e^{−ik′·r}V (r) (3(r̂ · σ σ σ_{1})(r̂ · σ σ σ_{2}) − σ σ σ_{1} · σ σ σ_{2}) e^{ik·r}dr =

− (3(q · σ σ σ_{1})(q · σ σ σ_{2}) − q^{2}σ σ σ_{1} · σ σ σ_{2}) \frac{1}{2π^{2}q^{2}} ∫_{0}^{∞} j_{2}(qr)V (r)r^{2}dr.  (1.6)

These expressions are used to represent the momentum-space interaction as a sum of scalar functions of q := |q| multiplied by spin-isospin operators. These scalar coefficient functions of the momentum transfer that multiply the spinisospin operators have the form of one of the integrals listed in Table 2:

Table 2: Momentum-space scalar coefficient functions

Scalar coefficient function dim indices
Ṽ_{m}(q) := \frac{1}{2π^{2}} ∫_{0}^{∞} j_{0}(qr)V_{m}(r)r^{2}dr MeV fm^{3} m ∈ \{1, 2, 3, 4, 15, 16, 18\}
Ṽ_{m}(q) := \frac{1}{2π^{2}q} ∫_{0}^{∞} j_{1}(qr)V_{m}(r)r^{3}dr MeV fm^{5} m ∈ \{7, 8, 9b, 10b, 11b, 12b, 13b, 14b\}
Ṽ_{m}(q) := \frac{1}{2π^{2}q^{2}} ∫_{0}^{∞} j_{2}(qr)V_{m}(r)r^{4}dr MeV fm^{7} m ∈ \{9a, 10a, 11a, 12a, 13a, 14a\}
Ṽ_{m}(q) := \frac{1}{2π^{2}q^{2}} ∫_{0}^{∞} j_{2}(qr)V_{m}(r)r^{2}dr MeV fm^{5} m ∈ \{5, 6, 17\}

where V_{m}(r) is the m^{th} potential in the expansion (1.1) and Ṽ_{ma}(q) and Ṽ_{mb}(q) are the two different functions that appear in (1.4) and (1.5). These functions have finite limits as q → 0 in spite of the 1/q^{l} coefficients since the Bessel function j_{l}(qr) vanishes like q^{l} as q → 0. The strong interaction contribution to the 24 scalar coefficients listed in Table 2 are numerically computed. The computational methods are discussed in section 3. Programs that compute these scalar coefficients are available as supplementary material to the electronic version of this paper. Quantities, like the binding energies in the test calculations, exhibit small sensitivities (in the sixth significant figure) to the precision of input constants. In the supplementary programs these constants are taken from the original V18 potential.

The electromagnetic contribution to each of these operators can be represented in terms of known special functions. These contributions are important for precise low-energy calculations and can be added to the strong interaction coefficient functions when they are needed. The analytic expressions for the electromagnetic terms are given in Appendix 2.

The resulting momentum-space potential has an operator expansion of the form

⟨k′|V |k⟩ = ∑_{m∈S} Ṽ_{m}(q)Õ_{m}  (1.7)

where S = {1, 2, 3, 4, 5, 6, 7, 8, 9a, 9b, 10a, 10b, 11a, 11b, 12a, 12b, 13a, 13b, 14a, 14b, 15, 16, 17,18} and the 24 operators Õ_{m} are given in Table 3.

Table 3: Argonne V18 momentum-space spin-isospin operators

term spin-isospin operator
Õ_{1} I
Õ_{2} (τ_{1} · τ_{2})
Õ_{3} (σ σ σ_{1} · σ σ σ_{2})
Õ_{4} (σ σ σ_{1} · σ σ σ_{2})(τ_{1} · τ_{2})
Õ_{5} − (3(q · σ σ σ_{1})(q · σ σ σ_{2}) − q^{2}σ σ σ_{1} · σ σ σ_{2})
Õ_{6} − (3(q · σ σ σ_{1})(q · σ σ σ_{2}) − q^{2}σ σ σ_{1} · σ σ σ_{2}) (τ_{1} · τ_{2})
Õ_{7} i(k × k′) · S
Õ_{8} i(k × k′) · S(τ_{1} · τ_{2})
Õ_{9a} −(k′ × k) · (k′ × k)
Õ_{9b} 2(k′ · k)
Õ_{10a} −(k′ × k) · (k′ × k)(τ_{1} · τ_{2})
Õ_{10b} 2(k′ · k)(τ_{1} · τ_{2})
Õ_{11a} −(k′ × k) · (k′ × k)(σ σ σ_{1} · σ σ σ_{2})
Õ_{11b} 2(k′ · k)(σ σ σ_{1} · σ σ σ_{2})
Õ_{12a} −(k′ × k) · (k′ × k)(σ σ σ_{1} · σ σ σ_{2})(τ_{1} · τ_{2})
Õ_{12b} 2(k′ · k)(σ σ σ_{1} · σ σ σ_{2})(τ_{1} · τ_{2})
Õ_{13a} −(S · (k × k′))^{2}
Õ_{13b} (k′ × S) · (k × S)
Õ_{14a} −(S · (k × k′))^{2}(τ_{1} · τ_{2})
Õ_{14b} (k′ × S) · (k × S)(τ_{1} · τ_{2})
Õ_{15} T_{12}
Õ_{16} (σ σ σ_{1} · σ σ σ_{2})T_{12}
Õ_{17} − (3(q · σ σ σ_{1})(q · σ σ σ_{2}) − q^{2}σ σ σ_{1} · σ σ σ_{2}) T_{12}
Õ_{18} (τ_{1z} + τ_{2z}).

The Argonne V18 potential in momentum-space has the dimension MeV fm^{3}. Dividing by ℏc in Mev-fermi can be used to convert the momentum-space potential to a consistent set of units, (fm)^{2}.

II. NUMERICAL FOURIER BESSEL TRANSFORMS

This section summarizes an accurate numerical computation of the integrals in Table 2. These computations are used to test the accuracy of the approximations discussed in the next section.

Because the configuration-space potential falls off asymptotically like e^{−m_{π}r}, the radial integrals are evaluated with a finite cutoff at 20 fm. The Fourier-Bessel transforms are evaluated for momentum transfers q < 100 fm^{−1}. With these cutoffs the maximum value of x := qr that can appear in the argument of the spherical Bessel functions in the integrals in Table 2. is x_{max} = 2000. To evaluate these integrals the zeros of the spherical Bessel functions j_{0}(x), j_{1}(x), and j_{2}(x) for 0 ≤ x ≤ 2000 are computed for each fixed value of q. For each value of q the integrals are expressed as sums of integrals between successive zeros of the spherical Bessel function that appear in the integral. If q is such that qr is never a zero of j_{l}(qr) for 0 < r < 20fm then the integral over r is performed using a 100 point Gauss-Legendre quadrature on the interval [0, 20fm]. If q is such that qr has zeros of j_{l}(qr) for 0 < r < 20fm, then the integrals between zeros [qr_{i}, qr_{i+1}] are computed using 20 Gauss-Legendre points when r_{i+1} ≤ 5 fm, 40 Gauss-Legendre points when 5 fm < r_{i+1} ≤ 10 fm and 80 Gauss-Legendre points when 10 fm < r_{i+1} ≤ 20 fm. For further details see [5]

III. APPROXIMATIONS

This section discusses two approximations of the potential functions Ṽ_{m}(q) in Table 2 by expansions in known elementary functions. The first method approximates these potential functions by linear combinations of Chebyshev polynomials on three distinct intervals of momenta, for momenta up to 100fm^{−1}. The second approach approximates these potential functions by a finite linear combination of orthonormal functions of the momentum transfer that have analytic Fourier-Bessel transforms. The configuration-space basis functions are associated Laguerre polynomials multiplied by decaying exponentials. These functions have analytic Fourier transforms that are rational functions of the momentum transfer [7]. In both approaches the coefficients of the expansion function are stored. The basis functions at any point can be generated efficiently by recursion and the potentials can be expressed as a finite linear combination of the basis functions. Both methods lead to efficient and accurate approximations to the momentum-space potential.

Figures 1 and 2 show the potential functions for the central and tensor parts (V_{1}(q) and V_{5}(q)) of the interaction to illustrate the structure of typical potentials.

A. Chebyshev expansions

This section discusses the Chebyshev basis. The functions Ṽ_{m}(q) are replaced by a Chebyshev polynomial approximation on the interval q ∈ [a, b] using [6]

Ṽ_{m}(q) ≈ c_{0}/2 + ∑_{n=1}^{100} c_{n}T_{n}(−\frac{a + b}{b − a} + \frac{2}{b − a}q)  (3.1)

where

T_{n}(x) = cos(n cos^{−1}(x)))  (3.2)

are Chebyshev polynomials and the coefficients c_{n} are computed using a Clenshaw-Curtiss quadrature [6] :

c_{n} = \frac{2}{N} [\frac{1}{2}Ṽ_{m}(b) + ∑_{j=1}^{N−1} Ṽ_{m}(\frac{a + b}{2} + \frac{b − a}{2} cos(πj/N)) cos(njπ/N) + (−)^{n} \frac{1}{2}Ṽ_{m}(a)] (3.3)

with N = 101. The functions Ṽ_{m}(q) are evaluated at the quadrature points q_{j} := \frac{a+b}{2} + \frac{b−a}{2} cos(πj/N) using the methods discussed above. This is repeated for q in each of three intervals, [a, b] = [0, 10], [10, 50], [50, 100] and the 101 expansion coefficients associated with each of these three intervals are stored. The Chebyshev polynomials are computed using the recurrence relations

T_{n+1}(x) = 2xT_{n}(x) − T_{n−1}(x),  T_{0}(x) = 1, T_{1}(x) = x.  (3.4)

For q larger than 100 fm^{−1} Ṽ_{m}(q) is approximated by 0.

For the potentials Ṽ_{4}(q), Ṽ_{6}(q) and Ṽ_{17}(q) it was necessary to add additional Chebyshev expansions intervals between zero and ten fm^{−1}. For Ṽ_{4}(q) 21 polynomials were used on [0, .2]fm^{−1}, 31 polynomials were used on [.2, .5]fm^{−1}, 41 polynomials were used on [.5, 2.0]fm^{−1} and 71 polynomials were used on [2.0, 10.0]fm^{−1}. For Ṽ_{6}(q) 31 polynomials were used on [0, .5]fm^{−1}, 41 polynomials were used on [.5, 2.0]fm^{−1} and 41 polynomials were used on [2.0, 10.0]fm^{−1}. Similarly for Ṽ_{17}(q) 31 polynomials were used on [0, 1.0]fm^{−1}, 51 polynomials were used on [1.0, 5.0]fm^{−1} and 51 polynomials were used on [5.0, 10.0]fm^{−1}.

This method provides an accurate and efficient representation for computing a momentum space V 18 interaction. One of the supplementary programs (chebyshev-argonne.c) uses this method to compute the 24 coefficient functions in Table 2.

B. Rational basis functions

While the method of the previous section gives accurate results, a more straightforward approach is to represent the potential directly as an expansion in basis functions that have analytic Fourier transforms. In order to represent the potential, each of the scalar potentials Ṽ_{m}(q), is approximated by an expansion in known basis functions. A method to compute both the expansion coefficients and a recursion formula to compute basis functions are given below.

The functions V_{m}(r), rV_{m}(r), and r^{2}V_{m}(r) that appear in the integrands of the integrals in Table 2 are expanded using an orthonormal set of radial functions that have analytic Fourier-Bessel transforms [7]. These functions are associated Laguerre polynomials multiplied by decaying exponentials in configuration space. Their Fourier-Bessel transforms have powerlaw fall of in momentum space. In addition, they vanish at the origin in a manner that can be used to explicitly cancel the factors 1/q and 1/q^{2} that appear in the definitions of Ṽ_{m} in Table 2. Both sets of basis functions can be generated efficiently using recursion relations. The cancellation of the factors 1/q and 1/q^{2} can be directly incorporated into the recursion that generates the momentum-space basis functions so the final expression for the potential does not require a special treatment for q near 0.

The radial basis functions for different values of l are given below. The dimensionless parameter x := Λr is used in the basis functions, where Λ is a scale parameter that can be chosen to improve efficiency. The parameterization of the Argonne V18 interaction uses the value Λ = 7(fm)^{−1}. The configuration-space basis functions are

ϕ_{nl}(r) = \frac{1}{\sqrt{N_{nl}}} x^{l}L_{n}^{2l+2}(2x)e^{−x}  (3.5)

where

L_{n}^{α} = ∑_{m=0}^{n}(−)^{m} ( \begin{matrix}n + α \\ n − m \end{matrix} ) \frac{x^{m}}{m!}  (3.6)

and the normalization coefficient is

N_{nl} = Λ^{−3}(\frac{1}{2})^{2l+3} \frac{Γ(n + α + 1)}{n!}.  (3.7)

These functions satisfy the orthogonality relations

∫_{0}^{∞} ϕ_{nl}(r)ϕ_{ml}(r)r^{2}dr = δ_{mn}.  (3.8)

They have analytic Fourier-Bessel transforms given by

ϕ̃_{nl}(q) = \sqrt{\frac{2}{π}} ∫_{0}^{∞} j_{l}(qr)ϕ_{nl}(r)r^{2}dr.  (3.9)

For y = q/Λ the ϕ̃_{nl}(q) can be expressed in terms of Jacobi polynomials:

ϕ̃_{nl}(q) = \frac{1}{\sqrt{Ñ_{nl}}} \frac{y^{l}}{(y^{2} + 1)^{l+2}} P_{n}^{l+ \frac{3}{2} ,l+ \frac{1}{2}} (\frac{y^{2} − 1}{y^{2} + 1})  (3.10)

with normalization coefficient

Ñ_{nl} = \frac{Λ^{3}}{2(2n + 2l + 3)} \frac{Γ(n + l + \frac{5}{2})Γ(n + l + \frac{3}{2})}{n!Γ(n + 2l + 3)}  (3.11)

and

P_{n}^{(α,β)}(x) = \frac{Γ(α + n + 1)}{n!Γ(α + β + n + 1)} ∑_{m=0}^{n} ( \begin{matrix}n \\ m \end{matrix} ) \frac{Γ(α + β + n + m + 1)}{2^{m}Γ(α + m + 1)}(x − 1)^{m}.  (3.12)

These functions satisfy the orthogonality relations

∫_{0}^{∞} ϕ̃_{nl}(q)ϕ̃_{ml}(q)q^{2}dq = δ_{mn}.  (3.13)

These basis functions can be generated by using the recursion formulas for the associated Laguerre functions and Jacobi polynomials

(n + 1)L_{n+1}^{α}(x) = (2n + α + 1 − x)L_{n+}^{α}(x) − (n + α)L_{n−1}^{α}(x)  (3.14)

and

2(n + 1)(n + α + β + 1)(2n + α + β)P_{n+1}^{(α,β)}(x) =

[(2n + α + β + 1)(α^{2} − β^{2}) + x((2n + α + β)(2n + α + β + 1)(2n + α + β + 2)]P_{n}^{(α,β)}(x)

−2(n + α)(n + β)(2n + α + β + 2)P_{n−1}^{(α,β)}(x).  (3.15)

These recursion relations can be modified to incorporate the normalization constants (3.7) and (3.11) directly into the recursion. The recursion for the normalized radial basis functions with (x = Λr) is given by:

ϕ_{0l}(r) = \frac{1}{\sqrt{(2l + 1)!}} \frac{1}{\sqrt{2^{2l+3}}} Λ^{3/2}x^{l}e^{−x}  (3.16)
ϕ_{1l}(r) = \frac{2l + 3 − 2x}{\sqrt{2l + 3}}ϕ_{0l}(r)  (3.17)
ϕ_{nl}(r) = \frac{2n + 1 + 2l − 2x}{\sqrt{n + 1 + 2l}\sqrt{n}}ϕ_{n−1,l}(r) − \sqrt{\frac{(n − 1)(n + 1 + 2l)}{n(n + 2 + 2l)}}ϕ_{n−2,l}(r).  (3.18)

Similarly, the normalized momentum-space basis functions with (y = q/Λ) are generated by the recursion:

ϕ̃_{0l}(q) =
\frac{1}{\sqrt{(2l + 3)!}} \frac{1}{\sqrt{Λ^{3}}} \frac{1}{\sqrt{\frac{1}{2} ⋯ \frac{2l+3}{2}}} \frac{1}{\sqrt{\frac{1}{2} ⋯ \frac{2l+1}{2}}} y^{l} \frac{1}{(y^{2} + 1)^{2l+2}}  (3.19)
ϕ̃_{1l}(q) = (\frac{1}{2} + (l + 2)\frac{y^{2} − 1}{y^{2} + 1})\sqrt{\frac{2l + 5}{(l + 2 + \frac{1}{2})(l + 1 + \frac{1}{2})}}ϕ̃_{0l}(q)  (3.20)
ϕ̃_{nl}(q) =
\sqrt{\frac{(2n + 2l + 3)n(n + 2l + 2)}{(2n + 2l + 1)(n + l + \frac{3}{2})(n + l + \frac{1}{2})}}

×\frac{(2n + 2l + 1)(2l + 2)(y^{2} + 1) + (2n + 2l + 1)(2n + 2l)(2n + 2l + 2)(y^{2} − 1)}{2n(n + 2l + 2)(2n + 2l)(y^{2} − 1)}ϕ̃_{n−1l}(q)

−\sqrt{\frac{(2n + 2l + 3)(n − 1)n(n + 1 + 2l)(n + 2l + 2)}{(2n + 2l − 1)(n + l + \frac{3}{2})(n + l + \frac{1}{2})(n + l + \frac{1}{2})(n + l − \frac{1}{2})}}

×\frac{(n + l + \frac{1}{2})(n + l − \frac{1}{2})(2n + 2l + 2)}{(n)(n + 2l + 2)(2n + 2l)}ϕ̃_{n−1l}(q).  (3.21)

Replacing ϕ̃_{0l}(q) in (3.19) by ϕ̂_{0l}(q) := ϕ̃_{0l}(q)/q^{l} given by

ϕ̂_{0l}(q) =

\frac{1}{\sqrt{(2l + 3)!}} \frac{1}{\sqrt{Λ^{3}}} \frac{1}{\sqrt{\frac{1}{2} ⋯ \frac{2l+3}{2}}} \frac{1}{\sqrt{\frac{1}{2} ⋯ \frac{2l+1}{2}}} Λ^{−l} \frac{1}{(y^{2} + 1)^{2l+2}}  (3.22)

to start the recursion in equations (3.20)-(3.21) generates ϕ̂_{nl}(q) := ϕ̃_{nl}(q)/q^{l}, which are wellbehaved as q → 0. Seventy expansion coefficients are used to construct the momentum-space potential for each value of m

c_{nm} = \frac{1}{2π^{2}} ∫_{0}^{∞} ϕ_{n0}(r)V_{m}(r)r^{2}dr m ∈ \{1, 2, 3, 4, 15, 16, 18\}  (3.23)
c_{nm} = \frac{1}{2π^{2}} ∫_{0}^{∞} ϕ_{n1}(r)V_{m}(r)r^{3}dr m ∈ \{7, 8, 9b, 10b, 11b, 12b, 13b, 14b\}  (3.24)
c_{nm} = \frac{1}{2π^{2}} ∫_{0}^{∞} ϕ_{n2}(r)V_{m}(r)r^{4}dr m ∈ \{9a, 10a, 11a, 12a, 13a, 14a\}  (3.25)
c_{nm} = \frac{1}{2π^{2}} ∫_{0}^{∞} ϕ_{n2}(r)V_{m}(r)r^{2}dr m ∈ \{5, 6, 17\}.  (3.26)

The integrals are approximated using an 80 point Gauss-Legendre quadrature between 0 and 10fm. The basis functions ϕ_{nl}(r) are generated using (3.16-3.18). The scale parameter in the recursion for ϕ_{nl}(r) is taken as Λ = 7fm^{−1}.

The 70x24 expansion coefficients c_{nm} are stored. The momentum-space potential functions are then given by

Ṽ_{m}(q) = ∑_{n=1}^{70} c_{nm}ϕ̂_{nl}(q)  (3.27)

where the reduced expansion functions ϕ̂_{nl}(q) := ϕ̃_{nl}(q)/q^{l} are generated recursively using (3.20-3.22).

The full momentum-space potential in operator form is given by

V = ∑_{m∈S} Ṽ_{m}(q)Õ_{m}  (3.28)

where Õ_{m} are the 24 operators in Table 3 and q = \sqrt{k^{2} + k^{′2} − 2k′ · k}.

One of the supplementary programs (rational-argonne.c) uses this method to compute the 24 coefficient functions in Table 2.

IV. TESTS

Two tests are performed on the potentials. First, the momentum-space coefficient functions, Ṽ_{m}(q), computed using the accurate numerical Fourier Bessel transforms, the Chebyshev expansion and the rational basis function expansion are compared. For the second test both representations of potential are used to compute the deuteron binding energy and wave functions. These results are compared to a direct calculation of these quantities using the partial-wave expansion of the original configuration space potential.

The results of the first test are shown in Tables 4-7, which list values of the Fourier-Bessel transforms of the 24 radial functions computed using these three different methods for momentum transfers of 1,5,15 and 25 fm^{−1}.

These results are shown in Tables 4,5,6 and 7 for all 24 operators and a representative range of the momentum transfers. The columns labeled RFExp show the scalar potential functions using the rational function expansion, the columns labeled CExp show the same quantities using the Chebyshev expansion, while the columns labeled NFT show the results of the direct numerical Fourier transform. Figures 1-24 plot the difference of the approximate Fourier transforms with an accurate Fourier Bessel transform divided by half of the sum of these quantities. The solid curves are for the rational function expansion and the dotted curved are for the Chebyshev expansion.

Table 4: Values of scalar coefficients at 1 fm^{−1}

n  RFExp  CExp  NFT
1 6.789973e-01 6.789977e-01 6.789977e-01
2 -4.019392e-01 -4.019392e-01 -4.019392e-01
3 -1.692090e-01 -1.692090e-01 -1.692090e-01
4 2.358519e-01 2.356704e-01 2.356705e-01
5 7.216739e-03 7.218217e-03 7.218217e-03
6 2.857732e-01 2.860471e-01 2.860467e-01
7 -5.511547e-01 -5.511547e-01 -5.511547e-01
8 -1.678888e-01 -1.678888e-01 -1.678888e-01
9 1.741415e-01 1.741415e-01 1.741415e-01
10 -3.272988e-02 -3.272987e-02 -3.272987e-02
11 1.999136e-02 1.999136e-02 1.999136e-02
12 -7.414060e-03 -7.414060e-03 -7.414060e-03
13 9.084422e-02 9.084424e-02 9.084424e-02
14 1.245017e-01 1.245017e-01 1.245017e-01
15 1.122388e-02 1.122389e-02 1.122389e-02
16 -1.214926e-02 -1.216021e-02 -1.216031e-02
17 2.403290e-03 2.420818e-03 2.420818e-03
18 6.124964e-03 6.124964e-03 6.124964e-03
19 1.304278e-02 1.304274e-02 1.304274e-02
20 -1.702409e-02 -1.702401e-02 -1.702401e-02
21 -7.227244e-03 -7.227256e-03 -7.227256e-03
22 -7.849686e-03 -7.849707e-03 -7.849707e-03
23 4.518193e-02 4.518262e-02 4.518262e-02
24 3.980251e-02 3.980269e-02 3.980269e-02

Table 5: Value of scalar coefficients at 5 fm^{−1}

n  RFExp  CExp  NFT
1 1.160699e+00 1.160699e-00 1.160699e-00
2 -1.360382e-02 -1.360382e-02 -1.360382e-02
3 -1.148807e-01 -1.148807e-01 -1.148807e-01
4 -1.065288e-01 -1.065203e-01 -1.065203e-01
5 4.489757e-03 4.489763e-03 4.489763e-03
6 4.405849e-03 4.405371e-03 4.405370e-03
7 -4.623736e-02 -4.623736e-02 -4.623736e-02
8 -1.871380e-02 -1.871380e-02 -1.871380e-02
9 2.471311e-02 2.471311e-02 2.471311e-02
10 1.480758e-03 1.480758e-03 1.480758e-03
11 6.027203e-03 6.027203e-03 6.027203e-03
12 1.465070e-03 1.465070e-03 1.465070e-03
13 5.222260e-03 5.222260e-03 5.222260e-03
14 8.233502e-03 8.233502e-03 8.233502e-03
15 4.828280e-03 4.828280e-03 4.828280e-03
16 -4.815794e-03 -4.815305e-03 -4.815305e-03
17 1.656921e-06 1.627533e-06 1.627533e-06
18 4.274306e-04 4.274306e-04 4.274306e-04
19 4.273833e-03 4.273832e-03 4.273832e-03
20 1.791462e-04 1.791461e-04 1.791461e-04
21 9.672551e-04 9.672550e-04 9.672550e-04
22 1.814761e-04 1.814761e-04 1.814761e-04
23 1.620319e-03 1.620318e-03 1.620318e-03
24 1.790086e-03 1.790085e-03 1.790085e-03

Table 6: Value of scalar coefficients at 15 fm^{−1}

n  RFExp  CExp  NFT
1 9.321365e-04 9.321031e-04 9.321031e-04
2 4.123439e-05 4.123387e-05 4.123387e-05
3 -1.924812e-05 -1.924669e-05 -1.924669e-05
4 -6.648375e-05 -6.643904e-05 -6.643770e-05
5 -9.010902e-06 -9.010512e-06 -9.010512e-06
6 1.026393e-05 1.026324e-05 1.026323e-05
7 5.541260e-06 5.540856e-06 5.540856e-06
8 2.632043e-06 2.631901e-06 2.631901e-06
9 -1.962835e-06 -1.962585e-06 -1.962585e-06
10 -9.304609e-07 -9.304846e-07 -9.304846e-07
11 -6.015901e-07 -6.015363e-07 -6.015363e-07
12 -1.047669e-07 -1.047529e-07 -1.047529e-07
13 -4.725022e-06 -4.725152e-06 -4.725152e-06
14 -1.527634e-06 -1.527584e-06 -1.527584e-06
15 2.942747e-06 2.942623e-06 2.942623e-06
16 -2.895027e-06 -2.892432e-06 -2.892244e-06
17 -2.865458e-10 -3.049671e-10 -3.061292e-10
18 9.986465e-08 9.985107e-08 9.985107e-08
19 -2.604487e-07 -2.604660e-07 -2.604660e-07
20 -6.335039e-08 -6.334951e-08 -6.334951e-08
21 -7.055132e-08 -7.055521e-08 -7.055521e-08
22 -1.454468e-08 -1.454569e-08 -1.454569e-08
23 -3.115148e-07 -3.115089e-07 -3.115089e-07
24 -1.394089e-07 -1.394129e-07 -1.394129e-07

Table 7: Value of scalar coefficients at 25 fm^{−1}

n RFExp CExp NFT
1 -1.386301e-05 -1.383431e-05 -1.383431e-05
2 -6.108349e-08 -6.010007e-08 -6.010007e-08
3 8.598072e-07 8.595154e-07 8.595154e-07
4 1.014189e-06 1.003839e-06 1.003915e-06
5 -4.600082e-07 -4.599210e-07 -4.599210e-07
6 4.739733e-07 4.738710e-07 4.738711e-07
7 2.443040e-08 2.442088e-08 2.442088e-08
8 9.428095e-09 9.412965e-09 9.412965e-09
9 -1.534834e-08 -1.533919e-08 -1.533919e-08
10 3.457372e-10 3.607579e-10 3.607579e-10
11 -3.619628e-09 -3.613201e-09 -3.613201e-09
12 -1.005784e-09 -1.003137e-09 -1.003137e-09
13 4.666338e-09 4.709390e-09 4.709390e-09
14 -3.274714e-09 -3.270324e-09 -3.270324e-09
15 -5.425469e-08 -5.415686e-08 -5.415686e-08
16 5.452722e-08 5.398602e-08 5.399643e-08
17 -2.888773e-12 -4.263050e-12 -4.213031e-12
18 -5.852151e-09 -5.841440e-09 -5.841440e-09
19 -2.512190e-10 -2.555613e-10 -2.555613e-10
20 7.827015e-12 7.411785e-12 7.411785e-12
21 -5.864134e-11 -5.980761e-11 -5.980761e-11
22 -1.617550e-11 -1.652711e-11 -1.652711e-11
23 8.297311e-11 8.271942e-11 8.271943e-11
24 -5.322500e-11 -5.424546e-11 -5.424546e-11

These tables show generally good agreement among the three methods of computation. At 1fm^{−1} and 5fm^{−1} the Chebyshev expansion agrees with the direct numerical Fourier transform to between 5-7 significant figures for all 24 potentials. There is similar agreement at 15fm^{−1} and 25fm^{−1} except in potentials 17. The agreement between the potentials calculated using the rational function expansion do not agree with the direct numerical Fourier transforms as well as the Chebyshev expansion. The accuracy depends on the particular potential and gets worse as the momentum transfer increases. Thus for precision calculations the Chebyshev expansion is preferred.

Figures 1-24 provide a more complete picture of the nature of the errors in both approximations. Spikes in the errors occur near points where the potentials change sign. Some of the errors near zero are enhanced because the some of the plotted potential are divided by powers of the momentum transfer. For these terms the operators include compensating powers of the momentum transfer that vanish near the origin, so the contribution of the error in the full potential near the origin is reduced. The rational function expansions have larger relative errors near higher and lower values of the momentum transfer. This is not surprising because the basis functions are not local. The Chebyshev expansion is uniformly good, in part because it is a local expansion, so more intervals can be added as needed. The largest errors are in potential 17. At 10fm^{−1} its value is about −1.2 × 10^{−8}, which is several orders of magnitude smaller than any of the other potentials at that momentum transfer.

Table 8: deuteron s and d wave functions using Chebyshev expansion, rational function expansion and r-space partial waves

kfm^{−1} u_{s}(k)-CExp. u_{s}(k)-RFExp u_{s}(k)-pw u_{d}(k)-CExp, u_{d}(k)-RFExp u_{d}(k)-pw
0.0 1.2695e+01 1.2695e+01 1.2693e+01 0.00000e+00 0.00000e+00 0.00000e+00
0.5 1.9609e+00 1.9609e+00 1.9609e+00 -2.19827e-01 -2.19811e-01 -2.19826e-01
1.0 3.7684e-01 3.7685e-01 3.7684e-01 -1.72164e-01 -1.72159e-01 -1.72164e-01
1.5 8.2472e-02 8.2472e-02 8.2471e-02 -1.12429e-01 -1.12429e-01 -1.12429e-01
2.0 6.0809e-03 6.0808e-03 6.0806e-03 -7.10857e-02 -7.10863e-02 -7.10857e-02
2.5 -1.3615e-02 -1.3616e-02 -1.3615e-02 -4.45428e-02 -4.45432e-02 -4.45428e-02
3.0 -1.6153e-02 -1.6153e-02 -1.6153e-02 -2.76853e-02 -2.76854e-01 -2.76853e-02
3.5 -1.3648e-02 -1.3648e-02 -1.3648e-02 -1.69880e-02 -1.69881e-02 -1.69881e-02
4.0 -1.0153e-02 -1.0153e-02 -1.0153e-02 -1.02233e-02 -1.02234e-02 -1.02234e-02
4.5 -6.9954e-03 -6.9954e-03 -6.9953e-03 -5.98472e-03 -5.98475e-03 -5.98472e-03
5.0 -4.5270e-03 -4.5270e-03 -4.5270e-03 -3.37040e-03 -3.37043e-03 -3.37041e-03

As a second test the deuteron binding energy and wave functions are computed using the two different momentum-space potentials are compared to each other and to the direct Fourier transform of wave functions computed using a configuration-space partial-wave calculation.

The deuteron binding energy and the s and d wave functions are computed using the operator form of the Fourier transformed potential, by direct integration of the vector variables. The method of solution, which is discussed in [8], uses the expansion (1.7) directly without using partial waves. Calculations are performed for both the Chebyshev and rationalfunction representations of the momentum space potentials.

These calculations are compared to a configuration-space partial-wave calculation. In that calculation, labeled pw in table 9, the wave functions are represented by an expansion in 70 configuration-space basis functions using the configuration-space basis functions (3.5). Matrix elements of the partial wave projection of the Hamiltonian, with the configuration space Argonne V18 potential, are computed in this basis and the eigenvalue problem is solved numerically. The Fourier transform is computed by analytically Fourier transforming the basis functions. The solution of the eigenvalue problem gives an independent evaluation of both the binding energy and wave functions constructed directly from the configuration space potential.

The deuteron binding energy obtained using the Chebyshev representation of the Fourier transform gives a deuteron binding energy of E_{d} = −2.242233 MeV. The rational function representation gives a d deuteron binding energy of E_{d} = −2.242193 MeV compared with E_{d} = −2.242211 MeV using the configuration space partial-wave calculation. The binding energies based on all three calculation agree to within 22 eV. The computation used in the configuration-space partial-wave calculation is a Galerkin calculation and thus gives a variational bound on the binding energy.

These eigenvalues differ from the experimental deuteron binding energy. This is because the momentum-space potentials used in these computations do not include electromagnetic corrections that appear in the Argonne V18 codes. The electromagnetic corrections contribute an additional +17.6 keV [9] to the binding energy which is consistent with the experimental binding energy of −2.2246 MeV.

The s and d wave functions for all three calculations are compared in Table 10. The wave functions differ in the fifth or sixth significant figure, while binding energies of all three calculations differ in the sixth significant figure.

The electromagnetic contributions to the Argonne V18 potential are important for lowenergy high-precision calculations. The Fourier transform of the electromagnetic contributions of the Argonne V18 potential can be computed analytically, and can be added to the strong interaction contribution discussed in this paper when necessary. The analytic Fourier transform of the electromagnetic contribution is discussed in appendix 3.

While this paper gives the momentum-space version of the operator expansion of the Argonne V18 potential, it is often useful to have partial-wave contributions of the potential. These can be computed from the operator matrix elements using a one-dimensional integration over the cosine of the angle between the initial and final momenta. A simple method to compute the partial-wave projections from the operator expansion is given in the appendix 2.

The programs to compute the potentials Ṽ_{m}(q) using both methods are available as supplementary material to the electronic version of this article.

This work supported by the U.S. Department of Energy, contract # DE-FG02-86ER40286 and National Science Foundation grants NSF-PHY-1005578 and NSF-PHY-1005501. The authors would like to acknowledge useful comments from Prof. C. Elster in preparing this manuscript.

V. APPENDIX 1

In this appendix we compute the Fourier transform of the parts of the potential containing the five types of operators, I, L · S, L · L,(L · S)^{2}, and S_{12} that appear in equations(1.3-1.6).

L · S:

Let q := k′ − k.

The Fourier transform is

\frac{1}{(2π)^{3}} ∫ e^{−ik′·r}V_{j}(r)L · Se^{ik·r}dr = \frac{1}{(2π)^{3}} ∫ e^{−ik′·r}V_{j}(r)S · (r × p)e^{ik·r}dr =

\frac{1}{(2π)^{3}} ∫ e^{−ik′·r}V_{j}(r)S · (r × k)e^{ik·r}dr = \frac{1}{(2π)^{3}} ∫ e^{−iq·r}V_{j}(r)S · (r × k)dr =

\frac{4π}{(2π)^{3}} ∫ ∑_{l=0}^{∞} ∑_{m=−l}^{l}(−i)^{l}j_{l}(qr)Y_{lm}(q̂)Y_{lm}^{∗}(r̂)V_{j}(r)S · (r × k)dr.  (5.1)

Since r can be expanded as a linear combination of spherical harmonics, Y_{1m}(r̂), the only terms that survive are the l = 1 terms. The integral over angles and the spherical harmonics replace r̂ by q̂, giving

−\frac{4πi}{(2π)^{3}} ∫_{0}^{∞} j_{1}(qr)V_{j}(r)S · (q × k)r^{3}dr =

iS · (k × k′) × [\frac{1}{2π^{2}q} ∫_{0}^{∞} j_{1}(qr)V_{j}(r)r^{3}dr].  (5.2)

Thus the Fourier transform has the form

\frac{1}{(2π)^{3}} ∫ e^{−ik′·r}V_{j}(r)L · Se^{ik·r}dr = iS · (k × k′)I_{1}(q)  (5.3)

where

I_{1}(q) := \frac{1}{2π^{2}q} ∫_{0}^{∞} j_{1}(qr)V_{j}(r)r^{3}dr.  (5.4)

The following relations are used to compute Fourier transforms of the remaining three operators:

∇ ∇ ∇_{q}f(q) = f′(q)\frac{q}{q}  (5.5)
∇_{q}^{2}f(q) = f′′(q) + \frac{2}{q}f′(q)  (5.6)
(a · ∇ ∇ ∇_{q})(b · ∇ ∇ ∇_{q})f(q) =
\frac{a · q}{q} \frac{b · q}{q}(f′′(q) − \frac{f′(q)}{q}) + \frac{a · b}{q}f′(q).  (5.7)

L · L:

The Fourier transform of this operator is

\frac{1}{(2π)^{3}} ∫ e^{−ik′·r}V_{j}(r)(r × p) · (r × p)e^{ik·r}dr = \frac{1}{(2π)^{3}} ∫ e^{−ik′·r}V_{j}(r)(r × k′) · (r × k)e^{ik·r}dr =

(i∇ ∇ ∇_{q} ×k′)·(i∇ ∇ ∇_{q} ×k)\frac{1}{(2π)^{3}} ∫ V_{j}(r)e^{−iq·r}dr = −(∇ ∇ ∇_{q} ×k′)·(∇ ∇ ∇_{q} ×k)\frac{4π}{(2π)^{3}} ∫_{0}^{∞} V_{j}(r)j_{0}(qr)r^{2}dr.

(5.8)

To compute the derivatives use

(∇ ∇ ∇_{q} × k′) · (∇ ∇ ∇_{q} × k) = (k · k′)(∇ ∇ ∇_{q} · ∇ ∇ ∇_{q}) − (k · ∇ ∇ ∇_{q})(k · ∇ ∇ ∇_{q})

in the above to get

−(∇ ∇ ∇_{q} × k′) · (∇ ∇ ∇_{q} × k)\frac{4π}{(2π)^{3}} ∫_{0}^{∞} V_{j}(r)j_{0}(qr)r^{2}dr = − ((k′ · k)∇_{q}^{2} − (k′ · ∇ ∇ ∇_{q})(k · ∇ ∇ ∇_{q})) I_{0}(q)

(5.9)

where

I_{0}(q) = \frac{4π}{(2π)^{3}} ∫_{0}^{∞} V_{j}(r)j_{0}(qr)r^{2}dr = \frac{1}{2π^{2}} ∫_{0}^{∞} V_{j}(r)j_{0}(qr)r^{2}dr.  (5.10)

Evaluating this gives

− ((k′ · k)∇_{q}^{2} − (k′ · ∇ ∇ ∇_{q})(k · ∇ ∇ ∇_{q})) I_{0}(q) =
−(k′ · k)(I′′_{0}(q) + \frac{1}{q}I′_{0}(q)) + \frac{(k′ · q)(k · q)}{q^{2}}(I′′_{0}(q) − \frac{1}{q}I′_{0}(q)).  (5.11)

To eliminate the derivatives use

I′′_{0}(q) − \frac{1}{q}I′_{0}(q) =

\frac{1}{2π^{2}} ∫_{0}^{∞} V (r)(j′′_{0}(qr) − j′_{0}(qr)\frac{1}{qr})r^{4}dr = \frac{1}{2π^{2}} ∫_{0}^{∞} V (r)j_{2}(qr)r^{4}dr = I_{2}(q)

and

I′′_{0}(q) + \frac{1}{q}I′_{0}(q) = \frac{1}{2π^{2}} ∫_{0}^{∞} V (r)(j′′_{0}(qr) − j′_{0}(qr)\frac{1}{qr} + 2j′_{0}(qr)\frac{1}{qr})r^{4}dr =

\frac{1}{2π^{2}} ∫_{0}^{∞} V (r)j′′_{2}(qr)r^{4}dr − \frac{1}{2π^{2}} \frac{2}{q} ∫_{0}^{∞} V (r)j_{1}(qr)r^{3}dr = I_{2}(q) − \frac{2}{q}I_{1}(q).

This gives

\frac{1}{(2π)^{3}} ∫ e^{−ik′·r}V_{j}(r)(r × p) · (r × p)e^{ik·r}dr = −(k′ · k)(I_{2}(q) − \frac{2}{q}I_{1}(q)) + \frac{(k′ · q)(k · q)}{q^{2}}I_{2}(q)

(5.12)

which can be reexpressed in terms of cross products:

\frac{1}{(2π)^{3}} ∫ e^{−ik′·r}V_{j}(r)(r×p)·(r×p)e^{ik·r}dr = −I_{2}(q)\frac{(k′ × k) · (k′ × k)}{q^{2}} + \frac{2}{q}(k′·k)I_{1}(q). (5.13)

(L · S)^{2}:

The Fourier transform is

\frac{1}{(2π)^{3}} ∫ e^{−i(k′−k)·r}V_{j}(r)(L · S)^{2}dr = \frac{1}{(2π)^{3}} ∫ e^{−i(k′−k)·r}V_{j}(r)(S · (r × p))^{2}dr =

−\frac{4π}{(2π)^{3}} ((k′ × S) · ∇ ∇ ∇_{q})((k × S) · ∇ ∇ ∇_{q}) ∫ j_{0}(qr)V_{j}(r)r^{2}dr =

−((k′ × S) · ∇ ∇ ∇_{q})((k × S) · ∇ ∇ ∇_{q})I_{0}(q) =

−((k′ × S) · q((k × S) · q))\frac{1}{q^{2}} (I′′_{0}(q) − \frac{1}{q}I′_{0}(q)) − (k′ × S) · (k × S)\frac{1}{q}I′_{0}(q) =

−((k′ × S) · q((k × S) · q)) = \frac{1}{q^{2}} I_{2}(q) + (k′ × S) · (k × S)\frac{1}{q}I_{1}(q)  (5.14)

which gives

−((S · (k × k′))^{2} \frac{1}{q^{2}} I_{2}(q) + (k′ × S) · (k × S)\frac{1}{q}I_{1}(q)  (5.15)

or

\frac{1}{(2π)^{2}} ∫ e^{−i(k′−k)}V_{j}(r)(L · S)^{2}dr = −((S · (k × k′))^{2} \frac{1}{q^{2}} I_{2}(q) + (k′ × S) · (k × S)\frac{1}{q}I_{1}(q). (5.16)
\frac{1}{3}S_{12} = (r̂ · σ σ σ_{1})(r̂ · σ σ σ_{2}) − \frac{1}{3}σ σ σ_{1} · σ σ σ_{2}):

The Fourier transform is

\frac{1}{(2π)^{3}} ∫ e^{−ik′·r}V (r) ((r̂ · σ σ σ_{1})(r̂ · σ σ σ_{2}) − \frac{1}{3}σ σ σ_{1} · σ σ σ_{2}) e^{ik·r}dr =

− ((∇ ∇ ∇_{q} · σ σ σ_{1})(∇ ∇ ∇_{q} · σ σ σ_{2}) − \frac{1}{3}σ σ σ_{1} · σ σ σ_{2}∇ ∇ ∇_{q}^{2}) \frac{4π}{(2π)^{3}} ∫ V (r)\frac{r^{2}}{r^{2}} j_{0}(qr)dr =

−\frac{σ σ σ_{1} · q}{q} \frac{σ σ σ_{2} · q}{q}((I′′_{0−}(q) − \frac{I′_{0−}(q)}{q}) − \frac{s_{1} · s_{2}}{q} I′_{0}(q) =

+\frac{1}{3}σ σ σ_{1} · σ σ σ_{2}(I′′_{0−}(q) + \frac{2}{q}I′_{0−}(q))  (5.17)

Thus,

\frac{1}{(2π)^{3}} ∫ e^{−ik′·r}V (r) (3(r̂ · σ σ σ_{1})(r̂ · σ σ σ_{2}) − σ σ σ_{1} · σ σ σ_{2}) e^{ik·r}dr = − (3\frac{σ σ σ_{1} · q}{q} \frac{σ σ σ_{2} · q}{q} − σ σ σ_{1} · σ σ σ_{2}) I_{2−}(q)

(5.18)

where

I_{0−}(q) = \frac{4π}{(2π)^{3}} ∫ V (r)j_{2}(qr)r^{2}dr.  (5.19)

VI. APPENDIX 2

In this appendix we compute the partial-wave projection of potentials from the vector representation of the momentum-space potential. Using rotational invariance the partialwave potentials can be expressed using a one-dimensional integral over the cosine of the angle between the initial and final momentum vectors. The method below is similar to the method first proposed in [10].

Rotational invariance of the potential implies

⟨j, μ, k, l, s|V |j′, μ′, k′, l′, s′⟩
= ⟨j, μ, k, l, s|U^{†}(R)V U(R)|j′, μ′, k′, l′, s′⟩

= ∑_{νν′} D_{νμ}^{j∗}(R)D_{ν′μ′}^{j}(R)⟨j, ν, k, l, s|V |j′, ν′, k′, l′, s′⟩

= ∫ dR ∑_{νν′} D_{νν′}^{j∗}(R)D_{ν′μ′}^{j}(R)⟨j, ν, k, l, s|V |j′, ν′, k′, l′, s′⟩

= δ_{jj′} ∑_{μ′ν′} δ_{μμ′} \frac{δ_{νν′}}{2j + 1}⟨j, ν, k, l, s|V |j′, ν′, k′, l′, s′⟩

= δ_{μμ′}δ_{jj′}⟨k, l, s∥V^{j}∥k′, l′, s′⟩  (6.1)

where we have integrated over the SU(2) Haar measure with normalization ∫ dR = 1 and defined the partial-wave potentials by

⟨k, l, s∥V^{j}∥k′, l′, s′⟩ := \frac{1}{2j + 1} ∑_{μ=−j}^{j} ∑_{m_{l}m′_{l}} ∑_{m_{s}m′_{s}}⟨j, μ, k, l, s|V |j, μ, k′, l′, s′⟩.  (6.2)

This kernel is rotationally invariant. Formally the partial-wave potential is

⟨k, l, s∥V^{j}∥k′, l′, s′⟩

= ∫ dk̂dk̂′ \frac{1}{2j + 1} ∑_{μ=−j}^{j}⟨\frac{1}{2}, μ_{1}, \frac{1}{2}, μ_{2}|s, m_{s}⟩⟨l, m_{l}, s, m_{s}|j, μ⟩Y_{lm_{l}}^{∗}(k̂)

×⟨k, μ_{1}, μ_{2}∥V ∥k′, μ′_{1}, μ′_{2}⟩
×⟨\frac{1}{2}, μ′_{1}, \frac{1}{2}, μ′_{2}|s′, m′_{s}⟩⟨l′, m′_{l}, s′, m′_{s}|j, μ⟩Y_{l′m′_{l}}(k̂′).  (6.3)

For any fixed rotation, R, rotational invariance of V gives

⟨k, μ_{1}, μ_{2}∥V ∥k′, μ′_{1}, μ′_{2}⟩

= ⟨k, μ_{1}, μ_{2}∥U^{†}(R)V U(R)∥k′, μ′_{1}, μ′_{2}⟩
= ∑_{μ′′_{1}μ′′_{2}} ∑_{μ′′′_{1}μ′′′_{2}} D_{μ′′_{1}μ_{1}}^{∗1/2}(R)D_{μ′′_{2}μ_{2}}^{∗1/2}(R)⟨Rk, μ′′_{1}, μ′′_{2}∥V ∥Rk′, μ′′′_{1}, μ′′′_{2}⟩D_{μ′′′_{1}μ′_{1}}^{1/2}(R)D_{μ′′′_{2}μ′_{2}}^{1/2}(R). (6.4)

Using this expression in eq. (6.3) gives

⟨k, l, s∥V^{j}∥k′, l′, s′⟩

= ∫ dk̂dk̂′ \frac{1}{2j + 1} ∑_{μ=−j}^{j}⟨\frac{1}{2}, μ_{1}, \frac{1}{2}, μ_{2}|s, m_{s}⟩⟨l, m_{l}, s, m_{s}|jμ⟩⟨Y_{lm_{l}}^{∗}(k̂)×

D_{μ′′_{1}μ_{1}}^{∗1/2}(R)D_{μ′′_{2}μ_{2}}^{∗1/2}(R)⟨Rk, μ′′_{1}, μ′′_{2}∥V ∥Rk′, μ′′′_{1}, μ′′′_{2}⟩D_{μ′′′_{1}μ′_{1}}^{1/2}(R)D_{μ′′′_{2}μ′_{2}}^{1/2}(R)×

⟨\frac{1}{2}, μ′_{1}, \frac{1}{2}, μ′_{2}|s′, m′_{s}⟩⟨l′, m′_{l}, s′, m′_{s}|jμ⟩⟨Y_{l′m′_{l}}(k̂′).  (6.5)

Next we eliminate three of the integrals in the potential matrix. For any fixed k′ there is an R such that R(k′)k̂′ = ẑ. Obviously R^{−1}(k′) = R_{z}(ϕ′)R_{y}(θ′)R_{z}(ξ) where (θ′, ϕ′) are the polar angles of k′ and ξ is arbitrary has this property where

R_{y}(θ) = ( \begin{matrix}cos(θ) & 0 & sin(θ) \\ 0 & 1 \\ − sin(θ) & 0 & cos(θ) \end{matrix} )  (6.6)

R_{z}(ϕ) = ( \begin{matrix}cos(ϕ) & − sin(ϕ) & 0 \\ sin(ϕ) & cos(ϕ) & 0 \\ 0 & 0 & 1 \end{matrix} ) .  (6.7)

In this case

Rk′ = k′ẑ  (6.8)

and

Rk = R_{z}^{−1}(ξ)R_{y}^{−1}(θ′)R_{z}^{−1}(ϕ′)k.  (6.9)

For fixed k′ define k′′ by

k′′ = R_{y}^{−1}(θ′)R_{z}^{−1}(ϕ′)k.  (6.10)

Given these fixed (primed) angles we change the unprimed integration variables k → k′′. We also have

Rk = R_{z}^{−1}(ξ)k′′  (6.11)

We are also free to choose the angle ξ in R_{z}^{−1}(ξ). We choose it so it transforms k′′ to the x − z plane. This is achieved by letting ξ be the azimuthal angle of k′′′

R_{z}^{−1}(ϕ′′) = ( \begin{matrix}cos(ϕ′′) & sin(ϕ′′) & 0 \\ − sin(ϕ′′) & cos(ϕ′′) & 0 \\ 0 & 0 & 1 \end{matrix} )  (6.12)

R_{z}^{−1}(ϕ′′)k′′ = ( \begin{matrix}cos(ϕ′′) & sin(ϕ′′) & 0 \\ − sin(ϕ′′) & cos(ϕ′′) & 0 \\ 0 & 0 & 1 \end{matrix} ) ( \begin{matrix}k′′ sin(θ′′) cos(ϕ′′) \\ k′′ sin(θ′′) sin(ϕ′′) \\ k′′ cos(θ′′) \end{matrix} ) = ( \begin{matrix}k′′ sin(θ′′) \\ 0 \\ k′′ cos(θ′′) \end{matrix} ) . (6.13)

With these substitutions the partial wave integral becomes

⟨k, l, s∥V^{j}∥k′, l′, s′⟩

= ∫ dk̂′′dk̂′ \frac{1}{2j + 1} ∑_{μ=−j}^{j}⟨\frac{1}{2}, μ_{1}, \frac{1}{2}, μ_{2}|s, m_{s}⟩⟨l, m_{l}, s, m_{s}|j, μ⟩Y_{lm}^{∗}(R^{−1}R_{z}^{−1}(ϕ′′)k̂′′)

×D_{μ′′_{1}μ_{1}}^{∗1/2}(R)D_{μ′′_{2}μ_{2}}^{∗1/2}(R)⟨k′′(x̂ sin(θ′′) + ẑ cos(θ′′)), μ′′_{1}, μ′′_{2}∥V ∥k′ẑ, μ′′′_{1}, μ′′′_{2}⟩D_{μ′′′_{1}μ′_{1}}^{1/2}(R)D_{μ′′′_{2}μ′_{2}}^{1/2}(R)

×⟨\frac{1}{2}, μ′_{1}, \frac{1}{2}, μ′_{2}|s′, m′_{s}⟩⟨l′, m′_{l}, s′, m′_{s}|j, μ⟩⟨Y_{l′m′_{l}}(R^{−1}ẑ′)

= ∫ dk̂′′dk̂′ \frac{1}{2j + 1} ∑_{μ=−j}^{j}⟨\frac{1}{2}, μ_{1}, \frac{1}{2}, μ_{2}|s, m_{s}⟩⟨l, m_{l}, s, m_{s}|j, μ⟩Y_{lm}^{∗}(x̂ sin(θ′′)+ẑ cos(θ′′))D_{m′′_{l}m_{l}}^{l∗}(R)

×D_{μ′′_{1}μ_{1}}^{∗1/2}(R)D_{μ′′_{2}μ_{2}}^{∗1/2}(R)⟨k′′(x̂ sin(θ′′) + ẑ cos(θ′′), μ′′_{1}, μ′′_{2}∥V ∥k′ẑ, μ′′′_{1}, μ′′′_{2}⟩D_{μ′′′_{1}μ′_{1}}^{1/2}(R)D_{μ′′′_{2}μ′_{2}}^{1/2}(R)

×⟨\frac{1}{2}, μ′_{1}, \frac{1}{2}, μ′_{2}|s′, m′_{s}⟩⟨l′, m′_{l}, s′, m′_{s}|jμ⟩Y_{l′m′_{l}}(ẑ)D_{m′′′_{l}m′_{l}}^{l′}(R).  (6.14)

Using properties of Clebsch-Gordan coefficients (i.e. D(R) < | >=< | > D(R) ⊗ D(R)) we get

⟨k, l, s∥V^{j}∥k′, l′, s′⟩

= ∫ dk̂′′dk̂′ \frac{1}{2j + 1} ∑_{μ=−j}^{j} D_{μ′′μ}^{j∗}(R)⟨\frac{1}{2}, μ_{1}, \frac{1}{2}, μ_{2}|s, m_{s}⟩⟨l, m_{l}, s, m_{s}|j, μ⟩Y_{lm_{l}}^{∗}(x̂ sin(θ′′)+ẑ cos(θ′′)|l, m_{l})

×⟨k′′(x̂ sin(θ′′) + ẑ cos(θ′′), μ_{1}, μ_{2}∥V ∥k′ẑ, μ′_{1}, μ′_{2}⟩

×⟨\frac{1}{2}, μ′_{1}, \frac{1}{2}, μ′_{2}|s′, m′_{s}⟩⟨l′, m′_{l}, s′, m′_{s}|j, μ′′′⟩Y_{l′m′_{L}}(ẑ)D_{μ′′′μ}^{j}(R)

= ∫ dk̂′′dk̂′ \frac{1}{2j + 1} ∑_{μ′′=−j}^{j}⟨\frac{1}{2}, μ_{1}, \frac{1}{2}, μ_{2}|s, m_{s}⟩⟨l, m_{l}, s, m_{s}|j, μ′′⟩Y_{lm_{l}}^{∗}(x̂ sin(θ′′) + ẑ cos(θ′′))

×⟨k′′(x̂ sin(θ′′) + ẑ cos(θ′′), μ_{1}, μ_{2}∥V ∥k′ẑ, μ′_{1}, μ′_{2}⟩
×⟨\frac{1}{2}, μ′_{1}, \frac{1}{2}, μ′_{2}|s′, m′_{s}⟩⟨l′, m′_{l}, s′, m′_{s}|j, μ′′⟩Y_{l′m′_{l}}(ẑ).  (6.15)

Since all of the dependence on ϕ′, θ′, ϕ′′ is in R, and the integrand is independent of R, these three angular integrals can be computed giving the multiplicative phase space factor of 8π^{2}. What remains is the following integral over the cosine of the angle between the final and initial momentum:

⟨k, l, s∥V^{j}∥k′, l′, s′⟩

= ∫_{−1}^{1} du′′ \frac{8π^{2}}{2j + 1} ∑_{μ′′=−j}^{j}⟨\frac{1}{2}, μ_{1}, \frac{1}{2}, μ_{2}|s, m_{s}⟩⟨l, m_{l}, s, m_{s}|j, μ′′⟩Y_{lm_{l}}^{∗}(x̂\sqrt{1 − u^{′′2}}) + ẑu′′)

×⟨k′′(x̂\sqrt{1 − u^{′′2}}) + ẑu′′), μ_{1}, μ_{2}∥V ∥k′ẑ), μ′_{1}, μ′_{2}⟩
×⟨\frac{1}{2}, μ′_{1}, \frac{1}{2}, μ′_{2}|s′, m′_{s}⟩⟨l′, m′_{l}, s′, m′_{s}|j, μ′′⟩Y_{l′m′_{l}}(ẑ).  (6.16)

The last thing that needs for be addressed for an explicit formula is the spherical harmonics

Y_{lm}^{∗}(x̂\sqrt{1 − u^{′′2}})+ẑu′′) = (−)^{m_{l}} \sqrt{\frac{2l + 1}{4π}}\sqrt{\frac{(l − m′_{l})!}{(l + m_{l})!}}P_{l}^{m}(u′′) = \sqrt{\frac{2l + 1}{4π}}D_{m0}^{l}(R_{y}(θ′′)) (6.17)
Y_{l′m′_{l}}(ẑ) = δ_{m′_{l}0}\sqrt{\frac{2l′ + 1}{4π}}.  (6.18)

Using these in the above expression we are left with a 1 dimensional integral

⟨k, l, s∥V^{j}∥k′, l′, s′⟩

= ∫_{−1}^{1} du′′ \frac{8π^{2}}{2j + 1} ∑_{μ′′=−j}^{j}⟨\frac{1}{2}, μ_{1}, \frac{1}{2}, μ_{2}|s, m_{s}⟩⟨l, m_{l}, s, m_{s}|j, μ′′⟩\sqrt{\frac{2l + 1}{4π}}D_{m0}^{l}(R_{y}(θ′′))

×⟨k′′(x̂\sqrt{1 − u^{′′2}}) + ẑu′′), μ_{1}, μ_{2}∥V ∥k′ẑ, μ′_{1}, μ′_{2}⟩
×⟨\frac{1}{2}, μ′_{1}, \frac{1}{2}, μ′_{2}|s′, m′_{s}⟩⟨l′, 0, s′, m′_{s}|j, μ′′⟩\sqrt{\frac{2l′ + 1}{4π}}.  (6.19)

Cleaning this up gives the following expression for the partial wave amplitude:

⟨j, μ, k, l, s|V |j′, μ′, k′, l′, s′⟩ =
δ_{jj′}δ_{μμ′}⟨k, l, s∥V^{j}∥k′, l′, s′⟩

with

⟨k, l, s∥V^{j}∥k′, l′, s′⟩

= \frac{8π^{2}}{2j + 1}\sqrt{\frac{2l + 1}{4π}}\sqrt{\frac{2l′ + 1}{4π}} ∑_{μ′′=−j}^{j} ∑⟨\frac{1}{2}, μ_{1}, \frac{1}{2}, μ_{2}|s, m_{s}⟩⟨l, m_{l}, s, m_{s}|jμ′′⟩

× ∫_{−1}^{1} du′′D_{m0}^{l}(R_{y}(cos^{−1}(u′′))⟨k′′(x̂\sqrt{1 − u^{′′2}}) + ẑu′′), μ_{1}, μ_{2}∥V ∥k′ẑ, μ′_{1}, μ′_{2}⟩

×⟨\frac{1}{2}, μ′_{1}, \frac{1}{2}, μ′_{2}|s′, m′_{s}⟩⟨l′, 0, s′, m′_{s}|jμ′′⟩  (6.20)

where all repeated spin indicies are summed. This reduces the partial-wave integral to a one-dimensional integral. In this case there are no traces, and no momentum-dependent spin bases, but there are a number of spin sums.

Explicit computation requires

D_{m0}^{l}(R_{y}(cos^{−1}(u′′)) = ∑_{s=0}^{2j} \frac{\sqrt{j + m)!(j − m)!j!j!}}{(j + m − s)!s!(s − m)!(j − s)!}R_{11}^{j+m−s}R_{12}^{s}R_{21}^{s−m}R_{22}^{j+s} (6.21)

where all negative factorials are ^{∞} and

R_{ij}(cos^{−1}(u′′)) = ( −\sqrt{\frac{1−u^{′′2}}{2}} \sqrt{\frac{1+u^{′′2}}{2}}^{\sqrt{\frac{1+u^{′′2}}{2}} \sqrt{\frac{1−u^{′′2}}{2}}} )  (6.22)

or

D_{m0}^{l}(R_{y}(cos^{−1}(u′′)) =
= ∑_{s=0}^{2j} \frac{\sqrt{j + m)!(j − m)!j!j!}}{(j + m − s)!s!(s − m)!(j − s)!}(−1)^{s−m+0}(\sqrt{\frac{1 + u^{′′2}}{2}} )^{2j+m}(\sqrt{\frac{1 − u^{′′2}}{2}} )^{2s−m}. (6.23)

VII. APPENDIX 3

The electromagnetic corrections to the nucleon-nucleon interaction have the following forms for the pp, np, and nn systems:

V_{em pp}(r) =

(v_{em,1}(r)+v_{em,2}(r)+v_{em,3}(r)+v_{em,4}(r))I +v_{em,6}(r)σ σ σ_{1} ·σ σ σ_{2} +v_{em,9}(r)S_{12} +v_{em,12}(r)L·S (7.1)
V_{em np}(r) = v_{em,5}(r)I + v_{em,8}(r)σ σ σ_{1} · σ σ σ_{2} + v_{em,11}(r)S_{12} + v_{em,14}(r)L · S  (7.2)
V_{em nn}(r) = v_{em,7}(r)σ σ σ_{1} · σ σ σ_{2} + v_{em,10}(r)S_{12} + v_{em,13}(r)L · S  (7.3)

The Fourier transforms have the same structure as they do for the strong interactions.

Ṽ_{em pp}(q) = (v_{em,1}(q) + v_{em,2}(q) + v_{em,3}(q) + v_{em,4}(q))I+
v_{em,6}(q)σ σ σ_{1} · σ σ σ_{2} + v_{em,9}(q)S̃_{12} + iv_{em,12}(q)(k × k′) · S.  (7.4)
Ṽ_{em np}(q) = v_{em,5}(q)I + v_{em,8}(q)σ σ σ_{1} · σ σ σ_{2} + v_{em,11}(q)S̃_{12} + iv_{em,14}(q)(k × k′) · S.  (7.5)

Ṽ_{em nn}(q) = v_{em,7}(q)σ σ σ_{1} · σ σ σ_{2} + v_{em,10}(q)S̃_{12} + iv_{em,13}(q)(k × k′) · S.  (7.6)

where

S̃_{12} := −3((q · σ σ σ_{1})(q · σ σ σ_{2}) − σ σ σ_{1} · σ σ σ_{2}).  (7.7)

The coefficients of I,σ σ σ_{1} · σ σ σ_{2} are

ṽ_{em,n}(q) = \frac{1}{2π^{2}} ∫_{0}^{∞} j_{0}(qr)v_{em,n}(r)r^{2}dr n ∈ \{1, 2, 3, 4, 5, 6, 7, 8\};  (7.8)

the coefficients of i(k × k′) · S are

ṽ_{em,n}(q) = \frac{1}{2π^{2}q} ∫_{0}^{∞} j_{1}(qr)v_{em,n}(r)r^{3}dr n ∈ \{12, 13, 14\}  (7.9)

and the coefficients of S̃_{12}, are

ṽ_{em,n}(q) = \frac{1}{2π^{2}q^{2}} ∫_{0}^{∞} j_{2}(qr)v_{em,n}(r)r^{2}dr n ∈ \{9, 10, 11\}  (7.10)

where the individual terms are

v_{em,1}(r) = \frac{αℏc}{r}(1 − (1 + \frac{11}{16}br + \frac{3}{16}(br)^{2} + \frac{1}{48}(br)^{3})e^{−br}),  (7.11)

v_{em,2}(r) = −\frac{α(ℏc)^{3}}{4m_{p}^{2}} \frac{b^{3}}{16}(1 + br + \frac{1}{3}(br)^{2})e^{−br}),  (7.12)

v_{em,3}(r) = −\frac{1}{m_{p}} v_{em}(1)^{2},  (7.13)

v_{em,4}(r) = \frac{2α}{3π}(−γ − \frac{5}{6} + | ln(ar)|)v_{em,1}(r),  (7.14)

v_{em,5}(r) = αℏcβ \frac{b^{3}}{384}(15 + 15br + 6(br)^{2} + (br)^{3})e^{−br}),  (7.15)

v_{em,6}(r) = −\frac{α(ℏc)^{3}μ_{p}^{2}}{6m_{p}^{2}} \frac{b^{3}}{16}(1 + br + \frac{1}{3}(br)^{2})e^{−br}),  (7.16)
v_{em,7}(r) = −\frac{α(ℏc)^{3}μ_{n}^{2}}{6m_{n}^{2}} \frac{b^{3}}{16}(1 + br + \frac{1}{3}(br)^{2})e^{−br}),  (7.17)
v_{em,8}(r) = −\frac{α(ℏc)^{3}μ_{n}μ_{p}}{6m_{n}m_{p}} \frac{b^{3}}{16}(1 + br + \frac{1}{3}(br)^{2})e^{−br}),  (7.18)

v_{em,9}(r) = −\frac{α(ℏc)^{3}μ_{p}^{2}}{4m_{p}^{2}} \frac{1}{r^{3}} (1 − (1 + br + \frac{1}{2}(br)^{2} + \frac{1}{6}(br)^{3} + \frac{1}{24}(br)^{4} + \frac{1}{144}(br)^{5})e^{−br}), (7.19)

v_{em,10}(r) = −\frac{α(ℏc)^{3}μ_{n}^{2}}{4m_{n}^{2}} \frac{1}{r^{3}} (1 − (1 + br + \frac{1}{2}(br)^{2} + \frac{1}{6}(br)^{3} + \frac{1}{24}(br)^{4} + \frac{1}{144}(br)^{5})e^{−br}), (7.20)

v_{em,11}(r) = −\frac{α(ℏc)^{3}μ_{p}μ_{n}}{4m_{p}m_{n}} \frac{1}{r^{3}} (1−(1+br + \frac{1}{2}(br)^{2} + \frac{1}{6}(br)^{3} + \frac{1}{24}(br)^{4} + \frac{1}{144}(br)^{5})e^{−br}), (7.21)

v_{em,12}(r) = −\frac{α(ℏc)^{3}(4μ_{p} − 1)μ_{p}μ_{n}}{2m_{p}^{2}} \frac{1}{r^{3}} (1 − (1 + br + \frac{1}{2}(br)^{2} + \frac{7}{48}(br)^{3})e^{−br}),  (7.22)

v_{em,13}(r) = 0,  (7.23)

and

v_{em,14}(r) = −\frac{α(ℏc)^{3}μ_{n}}{2m_{n}m_{r}} \frac{1}{r^{3}} (1 − (1 + br + \frac{1}{2}(br)^{2} + \frac{7}{48}(br)^{3})e^{−br})  (7.24)

where b = 4.27 and a = m_{e}/(ℏc), μ_{p} = 2.7928474, and μ_{n} = −1.9130427, β = .0189.

The Fourier transform of most of the terms in the potential can be computed using direct integration, the identities

j_{1}(x) = −\frac{d}{dx}j_{0}(x)  \frac{1}{x}j_{2}(x) = −\frac{d}{dx} \frac{j_{1}(x)}{x},  (7.25)

and the following relation [11], with ν = l + \frac{1}{2} and μ = n + 1/2, gives the relation

∫_{0}^{∞} e^{−br}j_{l}(qr)r^{n}dx

= \frac{q^{l}}{b^{n+l+1}} \frac{(n + l)!}{(2l + 1)!!}F (\frac{n + l + 1}{2}, \frac{n + l + 2}{2}, \frac{2l + 3}{2}, −q^{2}/b^{2})  (7.26)

which is valid for n + l > −1.

The only integral that can not be computed using these formulas involves the | ln(kr)| term that appears in v_{em 4}(r), which is an approximation to the vacuum polarization correction to the pp interaction. The required integrals, which are evaluated in this appendix are:

∫_{0}^{∞} j_{0}(qr)| ln(ar)|rdr = \frac{a^{2}}{q} [γ + 3 ln(q/a) − 2ci(q/a)],  (7.27)
∫_{0}^{∞} j_{0}(qr)| ln(ar)|re^{−br}dr =

\frac{a^{3}}{b^{2} + q^{2}} [b tan^{−1}(\frac{q}{b}) − qγ + \frac{q}{2} ln(q^{2}/a^{2} + b^{2}/a^{2})]+

\frac{a^{3}}{q^{2} + b^{2}} [q \frac{(E_{1}(\frac{b−iq}{a}) + E_{1}(\frac{b+iq}{a}))}{2} + b\frac{(E_{1}(\frac{b−iq}{a}) − E_{1}(\frac{b+iq}{a}))}{2i}] ,  (7.28)
∫_{0}^{∞} j_{0}(qr)| ln(ar)|r^{2}e^{−br}dr = −\frac{∂}{∂b} ∫_{0}^{∞} j_{0}(qr)| ln(ar)|re^{−br}dr,  (7.29)
∫_{0}^{∞} j_{0}(qr)| ln(ar)|r^{3}e^{−br}dr = \frac{∂^{2}}{∂b^{2}} ∫_{0}^{∞} j_{0}(qr)| ln(ar)|re^{−br}dr,  (7.30)
∫_{0}^{∞} j_{0}(qr)| ln(ar)|r^{4}e^{−br}dr = −\frac{∂^{3}}{∂b^{3}} ∫_{0}^{∞} j_{0}(qr)| ln(ar)|re^{−br}dr  (7.31)

and

∫_{0}^{∞} j_{2}(qr)r^{−1}dr = lim_{x→0} j_{1}(x)/x = 1/3.  (7.32)

The vacuum polarization integral only appears in the Coulomb potential which has the approximate form given in [12] - this approximation is adequate for binding energy calculations. It appears in the following contribution to the proton-proton interaction

v_{em(4)}(r) = \frac{2α}{3π}(−γ − \frac{5}{6} +| ln(ar)|+ \frac{6πkr}{8})\frac{αℏc}{r}(1−e^{−br}(1+ \frac{11}{16}br+ \frac{3}{16}(br)^{2}+ \frac{1}{48}(br)^{3}). (7.33)

The Fourier Bessel transform of this interaction,

\frac{1}{2π^{2}q} ∫_{0}^{∞} j_{0}(qr)V_{em(4)}(r)r^{2}dr =

\frac{1}{2π^{2}q^{2}} ∫_{0}^{∞} sin(qr)\frac{2α}{3π}(−γ − \frac{5}{6} +| ln(ar)|+ \frac{6πkr}{8})αℏc(1−e^{−br}(1+ \frac{11}{16}br+ \frac{3}{16}(br)^{2}+ \frac{1}{48}(br)^{3}), dr

(7.34)

can be computed analytically. There are three types of contributions

(I) = \frac{1}{2π^{2}q^{2}} ∫_{0}^{∞} sin(qr)\frac{2α}{3π}(−γ − \frac{5}{6})αℏc(1 − e^{−br}(1 + \frac{11}{16}br + \frac{3}{16}(br)^{2} + \frac{1}{48}(br)^{3})dr, (7.35)
(II) = \frac{1}{2π^{2}q^{2}} ∫_{0}^{∞} sin(qr)\frac{2α}{3π}| ln(ar)|αℏc(1 − e^{−br}(1 + \frac{11}{16}br + \frac{3}{16}(br)^{2} + \frac{1}{48}(br)^{3}), dr, (7.36)

and

(III) = \frac{1}{2π^{2}q^{2}} ∫_{0}^{∞} sin(qr)\frac{2α}{3π} \frac{6πkr}{8}αℏc(1 − e^{−br}(1 + \frac{11}{16}br + \frac{3}{16}(br)^{2} + \frac{1}{48}(br)^{3})dr. (7.37)

Integrals of the form (I) and (III) have the same form as the integrals discussed above. To calculate the integral (II) first replace r′ = ar → r = r′/a to get

(II) = \frac{1}{2aπ^{2}q^{2}} ∫_{0}^{∞} sin(qr′/a)\frac{2α}{3π}| ln(r′)|αℏc(1−e^{−br′/a}(1+ \frac{11}{16} \frac{br′}{a} + \frac{3}{16}(\frac{br′}{a})^{2}+ \frac{1}{48}(\frac{br′}{a})^{3}))dr′ =

= \frac{1}{2aπ^{2}q^{2}} \frac{2α}{3π}αℏc ∫_{0}^{∞} sin(qr′/a)| ln(r′)|dr′

−\frac{1}{2aπ^{2}q^{2}} \frac{2α}{3π}αℏc(1 + \frac{11}{16}(−b\frac{d}{db}) + \frac{3}{16}(b^{2} \frac{d^{2}}{d^{2}b}) + \frac{1}{48}(−b^{3} \frac{d^{3}}{d^{3}b})) ∫_{0}^{∞} sin(qr′/a)e^{−br′/a}| ln(r′)|dr′.

(7.38)

Two integrals need to be performed to compute this term. They are

∫_{0}^{∞} sin(qr′/a)| ln(r′)|dr′  (7.39)

and

∫_{0}^{∞} sin(qr′/a)e^{−br′/a}| ln(r′)|dr′.  (7.40)

The integral

∫_{0}^{∞} sin(qr′/a)| ln(r′)|dr′ =
− ∫_{0}^{1} sin(qr′/a) ln(r′)dr′ +lim_{d→0}[∫_{0}^{∞} sin(qr′/a)e^{−dr′} ln(r′)dr′ −∫_{0}^{1} sin(qr′/a) ln(r′)dr′]. (7.41)

These integrals can be found in [13]:

∫_{0}^{1} sin(qx) ln(x)dx = −\frac{1}{q}[γ + ln(q) − ci(q)]  (7.42)

where

ci(x) = Ci(x) = − ∫_{x}^{∞} \frac{cos(t)}{t}  (7.43)

and

∫_{0}^{∞} e^{−br} sin(qr) ln(r)dr = \frac{1}{b^{2} + q^{2}} [b tan^{−1}(\frac{q}{b}) − qγ + \frac{q}{2} ln(q^{2} + b^{2})].  (7.44)

The quantity γ is the Euler constant.

Thus the first of the required integrals needed to compute the vacuum polarization contribution is

∫_{0}^{∞} | ln(r′)| sin(qr′/a)dr′ =

− ∫_{0}^{1} sin(qr′/a) ln(r′)dr′ + lim_{d→0}[∫_{0}^{∞} e^{−dr′} sin(qr′/a) ln(r′)dr′ − ∫_{0}^{1} ln(r′) sin(qr′/a)] =

−2 ∫_{0}^{1} ln(r′) sin(qr′/a)dr′ + lim_{d→0} ∫_{0}^{∞} e^{−dr′} sin(qr′/a) ln(r′)dr′ =

\frac{2a}{q}[γ + ln(q/a) − ci(q/a)] + \frac{a}{q}[−γ + \frac{1}{2} ln(q^{2}/a^{2})] =

\frac{a}{q}[γ + 3 ln(q/a) − 2ci(q/a)].  (7.45)

Returning to the original expression - the first term in (II) is

\frac{1}{2aπ^{2}q^{2}} \frac{2α}{3π}αℏc ∫_{0}^{∞} sin(qr′/a)| ln(r′)|dr′

= \frac{α^{2}ℏc}{3π^{3}q^{3}} [γ + 3 ln(q/a) − 2ci(q/a)]  (7.46)

where γ is the Euler constant. Or

\frac{1}{2π^{2}q^{2}} ∫_{0}^{∞} sin(qr)\frac{2α}{3π}| ln(ar)|αℏc = \frac{α^{2}ℏc}{3π^{3}q^{3}} [γ + 3 ln(q/a) − 2ci(q/a)].  (7.47)

We also have to compute the second term in (II). The integral that is needed is

−\frac{α^{2}ℏc}{3π^{3}aq^{2}} ∫_{0}^{∞} sin(qr′/a)| ln(r′)|e^{−br′/a}(1 + \frac{11}{16} \frac{br′}{a} + \frac{3}{16}(\frac{br′}{a})^{2} + \frac{1}{48}(\frac{br′}{a})^{3})dr′.  (7.48)

This can be evaluated using

− ∫_{0}^{∞} sin(qr′/a)| ln(r′)|e^{−br′/a}dr′ =
∫_{0}^{∞} sin(qr′/a) ln(r′)e^{−br′/a}dr′ − 2 ∫_{1}^{∞} sin(qr′/a) ln(r′)e^{−br′/a}dr′  (7.49)

by differentiation with respect to b. The integral

∫_{1}^{∞} sin(qr′/a) ln(r′)e^{−br′/a}dr′ =

\frac{1}{2i} ∫_{1}^{∞}(e^{(iq/a−b/a)r′} − e^{(−iq/a−b/a)r′}) ln(r′)dr′ =

−\frac{a}{2i} ∫_{1}^{∞} \frac{dr′}{r′} [\frac{e^{(iq/a−b/a)r′}}{iq − b} + \frac{e^{(−iq/a−b/a)r′}}{iq + b}]

−\frac{a}{2i}[\frac{E_{1}(\frac{b−iq}{a})}{iq − b} + \frac{E_{1}(\frac{b+iq}{a})}{iq + b}]

\frac{a}{2i(q^{2} + b^{2})}[(iq + b)E_{1}(\frac{b − iq}{a}) + (iq − b)E_{1}(\frac{b + iq}{a}] =

\frac{a}{q^{2} + b^{2}} [[q \frac{(E_{1}(\frac{b−iq}{a}) + E_{1}(\frac{b+iq}{a}))}{2} + b\frac{(E_{1}(\frac{b−iq}{a}) − E_{1}(\frac{b+iq}{a}))}{2i}] .  (7.50)

Thus

∫_{1}^{∞} sin(qr′/a) ln(r′)e^{−br′/a}dr′
= \frac{a}{q^{2} + b^{2}} [[q \frac{(E_{1}(\frac{b−iq}{a}) + E_{1}(\frac{b+iq}{a}))}{2} + b\frac{(E_{1}(\frac{b−iq}{a}) − E_{1}(\frac{b+iq}{a}))}{2i}] .  (7.51)

Combining the two integrals gives

−\frac{α^{2}ℏc}{3π^{3}aq^{2}} ∫_{0}^{∞} sin(qr′/a)| ln(r′)|e^{−br′/a}dr′ =

−\frac{α^{2}ℏc}{3π^{3}aq^{2}} × \{

\frac{a}{b^{2} + q^{2}} [b tan^{−1}(\frac{q}{b}) − qγ + \frac{q}{2} ln(q^{2}/a^{2} + b^{2}/a^{2})]+

\frac{a}{q^{2} + b^{2}} [q \frac{(E_{1}(\frac{b−iq}{a}) + E_{1}(\frac{b+iq}{a}))}{2} + b\frac{(E_{1}(\frac{b−iq}{a}) − E_{1}(\frac{b+iq}{a}))}{2i}]\} .  (7.52)

Putting all of the parts together gives the ln(r) contribution to the vacuum polarization contribution

(II)_{b} = \frac{α^{2}ℏc}{3π^{3}q^{3}} [γ + 3 ln(q/a) − 2ci(q/a)]+

−(1 + \frac{11}{16}(−b\frac{d}{db}) + \frac{3}{16}(b^{2} \frac{d^{2}}{db^{2}} ) + \frac{1}{48}(−b^{3} \frac{d^{3}}{db^{3}} ))

\frac{α^{2}ℏc}{3π^{3}aq^{2}} × \{

\frac{a}{b^{2} + q^{2}} [b tan^{−1}(\frac{q}{b}) − qγ + \frac{q}{2} ln(q^{2}/a^{2} + b^{2}/a^{2})]+

\frac{a}{q^{2} + b^{2}} [[q \frac{(E_{1}(\frac{b−iq}{a}) + E_{1}(\frac{b+iq}{a}))}{2} + b\frac{(E_{1}(\frac{b−iq}{a}) − E_{1}(\frac{b+iq}{a}))}{2i}]\} .  (7.53)

This needs to be added to (I) and (II) to get the full vacuum polarization integral. These integrals can be computed using the methods used for all of the other potentials.

The exponential integrals have simple derivatives

\frac{d}{dx}E_{1}(x) = −\frac{e^{−x}}{x} .  (7.54)

[1] R.B. Wiringa, V.G.J. Stoks, R. Schiavilla, Phys. Rev. C51(1995)38.

[2] V.G.J. Stoks, R.A.M. Klomp, C.P.F. Terheggen, J. J. de Swart, Phys. Rev. C49(1994)2950.

[3] R. Machleidt, Phys. Rev. C63(2001)024001.

[4] W. Glöckle, H. Witala, D. Hüber, H. Kamada, J. Golak, Phys. Rep. 274 , 107 (1986).

[5] S. Veerasamy, University of Iowa Thesis, 2011.

[6] R. Broucke, Communications of the Association for Computing Machinery, 16,(1973)254.

[7] B. D. Keister and W. N. Polyzou, J. Comp. Phys. 134(1997),231.

[8] J. Golak, W. Glöckle , R. Skibiński, H. Witała, D. Rozpkedzik, K. Topolnicki, I. Fachruddin, Ch. Elster, A. Nogga, Phys. Rev. C81(2010)034006.

[9] Private communication with Robert Wiringa.

[10] http://arxiv.org/abs/arXiv:0911.4173 .

[11] I. S. Gradshteyn and I. M. Ryzhik, Tables of Integrals, Series, and Products, Associated Press, (1965), page 711 eq. 6.621.

[12] N. Auerbach, J. Hüfner, A. K. Kerman, and C. M. Shakin, Reviews of Modern Physics, 44,(1972)48.

[13] I. S. Gradshteyn and I. M. Ryzhik, Tables of Integrals, Series, and Products, Associated Press, (1965), page 583 eq. 4.381,page 605 eq. 4.401.

FIG. 1: V_{1}(k).

36

FIG. 2: V_{5}(k).

FIG. 3: ΔV_{1 rational}(k), ΔV_{1 Chebyshev}(k).

FIG. 11: ΔV_{9b rational}(k), ΔV_{9b Chebyshev}(k).

FIG. 21: ΔV_{9a rational}(k), ΔV_{9a Chebyshev}(k).

37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61