Introduction

Since Einstein’s theory of special relativity in the early twentieth century, there are still issues when considering systems that contain more than two nucleons, as in such systems, pair nucleons are influenced by the presence and motion of the other nucleons. As the three- and four-nucleon bound and scattering problems can be numerically solved with controlled errors, they provide an ideal theoretical laboratory for investigating the relativistic effects in the few-nucleon systems. Several techniques are developed to study the relativistic effects in few-body systems. Among them, the Faddeev–Yakubovsky method provides an exact numerical treatment of three and four-nucleon systems.

The inputs for the relativistic Faddeev–Yakubovsky equations are the fully-off-shell (FOS) relativistic 2N t-matrices, which can be obtained with two different approaches. In the first approach, the FOS relativistic 2N t-matrices are obtained directly from the nonrelativistic 2N t-matrices applying a two-step process. In the first step, using an analytical relation proposed by Coester et al., the relativistic right-half-shell (RHS) t-matrices are obtained from the nonrelativistic RHS t-matrices1. In the second step, the FOS relativistic t-matrices are obtained from the RHS t-matrices by solving a first resolvent equation proposed by Keister et al.2. This approach has been successfully implemented in few-body bound and scattering calculations in a three-dimensional (3D) scheme3,4,5,6,7,8, without using a partial wave (PW) decomposition.

In a second approach, the relativistic FOS t-matrices can be calculated from the solution of the relativistic Lippmann–Schwinger (LS) equation using the relativistic 2N potentials. The input relativistic NN potentials can be obtained from the nonrelativistic potentials by solving a quadratic equation, using an iterative scheme proposed by Kamada and Glöckle9. We have recently implemented this iterative technique in a 3D scheme to calculate the matrix elements of relativistic two-body (2B) potentials for spin-independent Malfliet-Tjon (MT) potential as a function of the magnitude of 2B relative momenta and the angle between them. To do so, we formulated the quadratic operator relation between the nonrelativistic and relativistic NN potentials in momentum space leading to a 3D integral equation10. We successfully implemented this iterative approach to calculate the matrix elements of boosted 2B potential from the MT potential to study the relativistic effects in a 3B bound state11. Our numerical results showed that the relativistic effects lead to a 2% reduction in 3B binding energy using MT potential. Our exact and detailed numerical studies of relativistic effects in 3B bound states using spin-independent MT potential demonstrates that direct integrations in the 3D scheme can be utilized to achieve the same results obtained using a PW method and paves the path for an extension to realistic interactions that have a more complicated spin–isospin dependence. Considering modern nucleon–nucleon (NN) potentials by including spin and isospin degrees of freedom and calculating the relativistic NN potentials from realistic NN potentials is the task we address in this paper. This is the first step toward our goal for a fully relativistic treatment of triton and Helium-3 bound state properties and the long-term interest in studying the scattering problems at the few-GeV energy scale in a 3D scheme. We show that the representation of the quadratic equation in momentum helicity basis states leads to a single and four coupled 3D integral equations for NN singlet and triplet spin state, correspondingly. The single and coupled integral equations are solved using the mentioned iterative scheme, and the matrix elements of relativistic NN potentials are obtained from CD-Bonn potential12. Our numerical analysis indicates that the calculated relativistic potential reproduces the deuteron binding energy and differential and total cross-sections of neutron-proton (np) elastic scattering with very high accuracy.

The motivation for using the 3D scheme and implementing a direct integration method is to replace the discrete angular momentum quantum numbers of a PW representation with continuous angle variables and consider all partial wave components to infinite order, independent of the energy scale of the problem. Consequently, 3D representation avoids the very involved angular momentum algebra for permutations, transformations, and few-nucleon forces, and in contrast to the PW approach, the number of equations in the 3D representation is energy independent. The 3D scheme is successfully implemented in a series of few-body bound and scattering states calculations by different few-body groups, from Ohio-Bochum collaboration4,5,7,13,14,15,16,17,18,19,20,21,22,23,24,25 to Tehran26,27,28,29,30,31,32 and Kraków33,34,35,36,37,38,39,40,41 groups.

In the “Relativistic NN potentials in a momentum helicity representation” section, we present the 3D formalism for the relationship between relativistic and nonrelativistic NN potentials. By projecting the quadratic relation between nonrelativistic and relativistic NN potentials in momentum helicity basis states, we obtain the matrix elements of relativistic NN potentials as a function of the magnitude of 2N relative momenta, the angle between them, and the helicity eigenvalues. We derive a single integral equation for NN total spin state \(s=0\) and four coupled integral equations for \(s=1\). In the “Calculation of relativistic NN interactions” section, we present our numerical results for the matrix elements of relativistic NN potential obtained from CD-Bonn potential in different spin and isospin channels. In the “Numerical tests for the relativistic NN potentials” section, we test the obtained relativistic potential by calculating and comparing deuteron binding energy and differential and total cross-sections of np elastic scattering with corresponding nonrelativistic results. Finally, a conclusion and outlook are provided in the “Conclusion and outlook” section.

Relativistic NN potentials in a momentum helicity representation

In this section, we show how to obtain the matrix elements of relativistic NN interactions in a 3D scheme from the nonrelativistic interactions by solving a nonlinear equation derived by Kamada and Glöckle. The relativistic interactions are designed to accurately reproduce the NN bound and scattering observables. To this aim and to check the accuracy of obtained relativistic interactions, as we show in the “Numerical tests for the relativistic NN potentials” section, one needs to solve the homogeneous and inhomogeneous LS integral equations (19), (21), and (22) in a momentum helicity representation to calculate relativistic deuteron binding energy and the scattering amplitudes to obtain the differential and total cross-sections.

The relativistic and nonrelativistic NN potentials, i.e., \({{\mathscr {V}}}_r\) and \(V_{nr}\), are related together by a quadratic operator equation as9

$$\begin{aligned} V_{nr}=\frac{1}{4m} \biggl (\omega ({\hat{p}}){{\mathscr {V}}}_r+ {{\mathscr {V}}}_r \omega ({\hat{p}})+{{\mathscr {V}}}_r^2 \biggr ), \end{aligned}$$
(1)

where m is the mass of nucleons, p is the relative momentum of two nucleons (\({\hat{p}}\) is the operator), and \(\omega (p)=2E(p)=2\sqrt{m^2+p^2}\). To calculate the relativistic NN potential \({{\mathscr {V}}}_r\) from a nonrelativistic potential \(V_{nr}\) in a 3D representation, we present Eq. (1) in momentum helicity basis states. The antisymmetrized momentum helicity basis states for a 2N system with total spin and isospin s and t, and the relative momentum \({\mathbf{p}}\) are introduced as16

$$\begin{aligned} | {\mathbf{p}}; \hat{{\mathbf{p}}}s \lambda ; t \rangle ^{\pi a}\equiv \frac{1}{\sqrt{2}} \biggl (1-\eta _{\pi }(-)^{s+t} \biggr ) \ |{\mathbf{p}};\hat{{\mathbf{p}}}\ s\lambda \rangle _{\pi }\ |t\rangle , \end{aligned}$$
(2)

where \(\hat{{\mathbf{p}}}\) is the unit momentum operator, \(\lambda \) is the eigenvalue of the helicity operator \({{\mathbf{s}}\cdot \hat{{\mathbf{p}}}} \), with the parity eigenvalues \(\eta _{\pi }=\pm 1\) and eigenstates \(|{\mathbf{p}};\hat{{\mathbf{p}}}\ s\lambda \rangle _{\pi }=\frac{1}{\sqrt{2}} (1+\eta _{\pi } P_{\pi }) |{\mathbf{p}};\hat{{\mathbf{p}}}\ s\lambda \rangle \). The 2N helicity basis states are orthogonal and normalized as

$$\begin{aligned}{}& {}^{\pi' a}\langle {\mathbf{p}}';\hat{{\mathbf{p}}}' s' \lambda ' ; t' | {\mathbf{p}}; \hat{{\mathbf{p}}}s \lambda ; t \rangle ^{\pi a}= \biggl (1-\eta _{\pi }(-)^{s+t}\biggr )\ \delta _{\eta _{\pi '} \eta _{\pi }} \delta _{s's} \delta _{t't} \ \biggl \{\delta ({\mathbf{p}}'-{\mathbf{p}})\delta _{\lambda \lambda '}+\eta _{\pi }(-)^{s} \delta ({\mathbf{p}}'+{\mathbf{p}})\delta _{\lambda '-\lambda } \biggr \}, \\&\sum _{s\lambda \pi t}\int d{\mathbf{p}}\ | {\mathbf{p}}; \hat{{\mathbf{p}}}s \lambda ; t \rangle ^{\pi a}\ \frac{1}{4}\ \ ^{\pi a}\langle {\mathbf{p}}; \hat{{\mathbf{p}}}s \lambda ; t |=1. \end{aligned}$$
(3)

The matrix elements of NN nonrelativistic and relativistic potentials in 2N helicity basis states, introduced in Eq. (2), are given as

$$\begin{aligned} V^{\pi s t}_{nr,\lambda \lambda '}({\mathbf{p}},{\mathbf{p}}')\equiv & {} \ ^{\pi a}\langle {\mathbf{p}}; \hat{{\mathbf{p}}}s \lambda ; t |V_{nr} | {\mathbf{p}}' ; \hat{{\mathbf{p}}}' s \lambda ' ; t \rangle ^{\pi a}, \end{aligned}$$
(4)
$$\begin{aligned} {{\mathscr {V}}}^{\pi s t}_{r,\lambda \lambda '}({\mathbf{p}},{\mathbf{p}}')\equiv & {} \ ^{\pi a}\langle {\mathbf{p}}; \hat{{\mathbf{p}}}s \lambda ; t |{{\mathscr {V}}}_r| {\mathbf{p}}' ; \hat{{\mathbf{p}}}' s \lambda ' ; t \rangle ^{\pi a}. \end{aligned}$$
(5)

Representation of the quadratic relation between nonrelativistic and relativistic NN potentials, given in Eq. (1), in 2N helicity basis states is as

$$\begin{aligned} \ ^{\pi a}\langle {\mathbf{p}}; \hat{{\mathbf{p}}}s \lambda ; t |4m V_{nr} | {\mathbf{p}}' ; \hat{{\mathbf{p}}}' s \lambda ' ; t \rangle ^{\pi a}= \,& {} \ ^{\pi a}\langle {\mathbf{p}}; \hat{{\mathbf{p}}}s \lambda ; t |\omega ({\hat{p}}){{\mathscr {V}}}_r| {\mathbf{p}}' ; \hat{{\mathbf{p}}}' s \lambda ' ; t \rangle ^{\pi a} \\&+\, \ ^{\pi a}\langle {\mathbf{p}}; \hat{{\mathbf{p}}}s \lambda ; t |{{\mathscr {V}}}_r\omega ({\hat{p}})| {\mathbf{p}}' ; \hat{{\mathbf{p}}}' s \lambda ' ; t \rangle ^{\pi a} \\&+\, \ ^{\pi a}\langle {\mathbf{p}}; \hat{{\mathbf{p}}}s \lambda ; t |{{\mathscr {V}}}_r \cdot {{\mathscr {V}}}_r| {\mathbf{p}}' ; \hat{{\mathbf{p}}}' s \lambda ' ; t \rangle ^{\pi a}. \end{aligned}$$
(6)

The first and the second terms on the right-hand side of Eq. (6) can be evaluated straightforwardly. For evaluation of the third term, the completeness relation of Eq. (3) should be inserted. By these considerations, Eq. (6) reads as

$$\begin{aligned} \frac{4m V^{\pi s t}_{nr, \lambda \lambda '}({\mathbf{p}},{\mathbf{p}}')}{\omega (p) + \omega (p')}= & {} {{\mathscr {V}}}^{\pi s t}_{r,\lambda \lambda '}({\mathbf{p}},{\mathbf{p}}') + \frac{1}{ \omega (p) + \omega (p')} \frac{1}{4} \sum _{\lambda ''=-1}^{+1}\int d{\mathbf{p}}'' \ {{\mathscr {V}}}^{\pi s t}_{r,\lambda \lambda ''}({\mathbf{p}},{\mathbf{p}}'') {{\mathscr {V}}}^{\pi s t}_{r,\lambda ''\lambda '}({\mathbf{p}}'',{\mathbf{p}}'). \end{aligned}$$
(7)

By considering the following properties of NN potentials, one can obtain the negative helicity eigenvalue components of the potential from the positive ones as

$$\begin{aligned} V^{\pi s t}_{nr, -\lambda \lambda '}({\mathbf{p}},{\mathbf{p}}')= \,& {} \eta _\pi (-)^s \ V^{\pi s t}_{nr, \lambda \lambda '}(-{\mathbf{p}},{\mathbf{p}}'), \\ V^{\pi s t}_{nr, \lambda -\lambda '}({\mathbf{p}},{\mathbf{p}}')=\, & {} \eta _\pi (-)^s \ V^{\pi s t}_{nr, \lambda \lambda '}({\mathbf{p}},-{\mathbf{p}}'). \end{aligned}$$
(8)

The above relations are also valid for the relativistic potential \({{\mathscr {V}}}^{\pi s t}_{r,\lambda \lambda ''}({\mathbf{p}},{\mathbf{p}}'')\). By using the properties of Eq. (8), one can show that the second term of Eq. (7), for \(\lambda ''=-1\) and \(\lambda ''=+1\) are equal together

$$\begin{aligned} \int d{\mathbf{p}}'' \ {{\mathscr {V}}}^{\pi s t}_{r,\lambda -1}({\mathbf{p}},{\mathbf{p}}'') {{\mathscr {V}}}^{\pi s t}_{r,-1 \lambda '}({\mathbf{p}}'',{\mathbf{p}}')=\, & {} \int d{\mathbf{p}}'' \ \eta _\pi (-)^s \ {{\mathscr {V}}}^{\pi s t}_{r,\lambda 1}({\mathbf{p}},-{\mathbf{p}}'') \eta _\pi (-)^s \ {{\mathscr {V}}}^{\pi s t}_{r,1 \lambda '}(-{\mathbf{p}}'',{\mathbf{p}}') \\= \,& {} \underbrace{\eta _\pi ^2 (-)^{2s}}_1 \int d{\mathbf{p}}'' \ {{\mathscr {V}}}^{\pi s t}_{r,\lambda 1}({\mathbf{p}},-{\mathbf{p}}'') {{\mathscr {V}}}^{\pi s t}_{r,1 \lambda '}(-{\mathbf{p}}'',{\mathbf{p}}') \\= \,& {} \int d{\mathbf{p}}'' \ {{\mathscr {V}}}^{\pi s t}_{r,\lambda 1}({\mathbf{p}},{\mathbf{p}}'') {{\mathscr {V}}}^{\pi s t}_{r,1 \lambda '}({\mathbf{p}}'',{\mathbf{p}}'). \end{aligned}$$
(9)

For 2N singlet spin state, \(s=0\), Eq. (7) leads to a single integral equation to obtain the matrix elements of relativistic potential \({{\mathscr {V}}}^{\pi 0 t}_{r,00}({\mathbf{p}},{\mathbf{p}}')\)

$$\begin{aligned} \frac{4m V^{\pi 0 t}_{nr, 00}({\mathbf{p}},{\mathbf{p}}')}{\omega (p) + \omega (p')}= & {} {{\mathscr {V}}}^{\pi 0 t}_{r,00}({\mathbf{p}},{\mathbf{p}}') + \frac{1}{ \omega (p) + \omega (p')} \frac{1}{4} \int d{\mathbf{p}}'' \ {{\mathscr {V}}}^{\pi 0 t}_{r,00}({\mathbf{p}},{\mathbf{p}}'') {{\mathscr {V}}}^{\pi 0 t}_{r,00}({\mathbf{p}}'',{\mathbf{p}}'). \end{aligned}$$
(10)

For triplet spin states, \(s=1\), Eq. (7) leads to four coupled integral equations, corresponding to helicity eigenvalues \(\lambda ,\lambda '=0,+1\), to calculate the matrix elements of relativistic potentials \({{\mathscr {V}}}^{\pi 1 t}_{r,00}({\mathbf{p}},{\mathbf{p}}')\), \({{\mathscr {V}}}^{\pi 1 t}_{r,01}({\mathbf{p}},{\mathbf{p}}')\), \({{\mathscr {V}}}^{\pi 1 t}_{r,10}({\mathbf{p}},{\mathbf{p}}')\), and \({{\mathscr {V}}}^{\pi 1 t}_{r,11}({\mathbf{p}},{\mathbf{p}}')\)

$$\begin{aligned} \frac{4m V^{\pi 1 t}_{nr, \lambda \lambda '}({\mathbf{p}},{\mathbf{p}}')}{\omega (p) + \omega (p')}=\, & {} {{\mathscr {V}}}^{\pi 1 t}_{r,\lambda \lambda '}({\mathbf{p}},{\mathbf{p}}') \\&+\, \frac{1}{ \omega (p) + \omega (p')} \biggl ( \frac{1}{2} \int d{\mathbf{p}}'' \ {{\mathscr {V}}}^{\pi 1 t}_{r,\lambda 1}({\mathbf{p}},{\mathbf{p}}'') {{\mathscr {V}}}^{\pi 1 t}_{r,1\lambda '}({\mathbf{p}}'',{\mathbf{p}}') +\frac{1}{4} \int d{\mathbf{p}}'' \ {{\mathscr {V}}}^{\pi 1 t}_{r,\lambda 0}({\mathbf{p}},{\mathbf{p}}'') {{\mathscr {V}}}^{\pi 1 t}_{r,0\lambda '}({\mathbf{p}}'',{\mathbf{p}}') \biggr ) . \end{aligned}$$
(11)

Equations (10) and (11) should be solved for each value of 2N total isospin \(t=0, 1\). For the numerical solution of single and coupled integral equations, i.e., Eqs. (10) and (11), by choosing momentum vector \({\mathbf{p}}'\) parallel to the \(z-\)axis, the azimuthal angular dependence of the matrix elements of nonrelativistic and relativistic potentials can be factored out as an exponential phase as

$$\begin{aligned} V^{\pi s t}_{nr, \lambda \lambda '}({\mathbf{p}},p'{\hat{z}})=\, & {} e^{{i\lambda ' \phi } } V^{\pi s t}_{nr, \lambda \lambda '}(p,p',x) , \\ {{\mathscr {V}}}^{\pi s t}_{r,\lambda \lambda '}({\mathbf{p}},p'{\hat{z}})= \,& {} e^{{i\lambda ' \phi } } {{\mathscr {V}}}^{\pi s t}_{r,\lambda \lambda '}(p,p',x) , \\ {{\mathscr {V}}}^{\pi s t}_{r,\lambda '' \lambda '}({\mathbf{p}}'',p'{\hat{z}})=\, & {} e^{{i\lambda ' \phi ''} } {{\mathscr {V}}}^{\pi s t}_{r,\lambda '' \lambda '}(p'',p',x''), \end{aligned}$$
(12)

where \(x \equiv \hat{{\mathbf{p}}}\cdot \hat{{\mathbf{p}}}'\) and \(x'' \equiv \hat{{\mathbf{p}}}'' \cdot \hat{{\mathbf{p}}}'\). One can show that the matrix elements of \({{\mathscr {V}}}^{\pi s t}_{r,\lambda \lambda ''} ({\mathbf{p}},{\mathbf{p}}'')\) can be obtained as

$$\begin{aligned} {{\mathscr {V}}}^{\pi s t}_{r,\lambda \lambda ''} ({\mathbf{p}},{\mathbf{p}}'' ) = \frac{{{\mathscr {V}}}^{\pi s t}_{r,\lambda \lambda ''}(p,p'',\gamma )}{d_{\lambda '' \lambda }^s (\gamma )} \sum _{\lambda ^*=-s}^{+s} e^{i\lambda ^*(\phi -\phi '')} d_{\lambda ^* \lambda } ^s (x) \ d_{\lambda ^* \lambda ''} ^s (x'') , \end{aligned}$$
(13)

where \(\gamma = x x'' + \sqrt{1-x^2} \sqrt{1-x''^2} \cos {(\phi -\phi '')}\). By these considerations, Eqs. (10) and (11) can be written as

$$\begin{aligned} \frac{4m V^{\pi 0 t}_{nr, 00}(p,p',x)}{\omega (p) + \omega (p')}= & {} {{\mathscr {V}}}^{\pi 0 t}_{r,00}(p,p',x) + \frac{1}{ \omega (p) + \omega (p')} \frac{1}{4} \int _{0}^{\infty } dp'' p''^2 \int _{-1}^{+1} dx'' \ {{\mathscr {V}}}^{{\pi 0 t}, 0}_{r,0 0} (p,p'',x,x'') \ {{\mathscr {V}}}^{\pi 0 t}_{r, 0 0}(p'',p',x'') , \end{aligned}$$
(14)
$$\begin{aligned} \frac{4m V^{\pi 1 t}_{nr, \lambda \lambda '}(p,p',x)}{\omega (p) + \omega (p')}=\, & {} {{\mathscr {V}}}^{\pi 1 t}_{r,\lambda \lambda '}(p,p',x) \\&\,+ \frac{1}{ \omega (p) + \omega (p')} \int _{0}^{\infty } dp'' p''^2 \int _{-1}^{+1} dx'' \\&\times \,\biggl \{ \frac{1}{4} {{\mathscr {V}}}^{{\pi 1 t, \lambda '}}_{r,\lambda 0} (p,p'',x,x'') \ {{\mathscr {V}}}^{\pi 1 t}_{r, 0 \lambda '}(p'',p',x'') + \frac{1}{2} {{\mathscr {V}}}^{{\pi 1 t, \lambda '}}_{r,\lambda 1} (p,p'',x,x'') \ {{\mathscr {V}}}^{\pi 1 t}_{r, 1 \lambda '}(p'',p',x'') \biggr \} , \end{aligned}$$
(15)

where

$$\begin{aligned} {{\mathscr {V}}}^{{\pi s t, \lambda '}}_{r,\lambda \lambda ''} (p,p'',x,x'')= & {} \sum _{\lambda ^*=-s}^{+s} d_{\lambda ^* \lambda } ^s (x) \ d_{\lambda ^* \lambda ''} ^s (x'') \int _0^{2\pi } d\phi '' \ e^{i (\lambda ^*- \lambda ' ) (\phi -\phi '')} \frac{ {{\mathscr {V}}}^{\pi s t}_{r,\lambda \lambda ''}(p,p'',\gamma ) }{d_{\lambda '' \lambda }^s (\gamma )}. \end{aligned}$$
(16)

The \(\phi ''\) integration in the interval \([0,2\pi ]\) can be simplified to \([0,\pi /2]\) interval by

$$\begin{aligned} \int _0^{2\pi } d\phi '' \ e^{im (\phi -\phi '')} {\mathscr{F}}\bigl (\cos (\phi -\phi '')\bigr )=\, & {} \int _0^{2\pi } d\phi '' \ e^{im \phi ''} {{\mathscr {F}}}\bigl (\cos \phi ''\bigr ) \\= \,& {} 2 \int _0^{\pi /2} d\phi '' \cos m\phi '' \biggl [ {{\mathscr {F}}}\bigl (\cos \phi '' \bigr ) + (-)^m {{\mathscr {F}}}\bigl (-\cos \phi '' \bigr ) \biggr ] . \end{aligned}$$
(17)

Calculation of relativistic NN interactions

To calculate the matrix elements of relativistic potentials for different spin–isospin (s, t) channels, we solve the integral equations (14) and (15) using an iterative method proposed by Kamada and Glöckle9. The iteration starts with \({{\mathscr {V}}}^{(0){\pi s t}}_{r,\lambda \lambda '}(p,p',x) = \dfrac{4m V^{\pi s t}_{nr, \lambda \lambda '}(p,p',x)}{\omega (p) + \omega (p')}\) and stops when the maximal difference between the matrix elements of the relativistic potential \({{\mathscr {V}}}^{{\pi s t}}_{r,\lambda \lambda '}(p,p',x)\) obtained from two successive iterations drops below \(10^{-6}\) MeV fm\(^3\). After each iteration, to obtain the matrix elements \({{\mathscr {V}}}^{{\pi s t, \lambda '}}_{r,\lambda \lambda ''} (p,p'',x,x'')\) that appears in the kernel of Eqs. (14) and (15), we need to perform the azimuthal angle integration of Eq. (16). To speed up the convergence of the iteration in solving Eqs. (14) and (15), in some (s, t) channels even to be able to reach the convergence, we use the weighted average of relativistic potential obtained from two successive iterations as

$$\begin{aligned} {{\mathscr {V}}}_r^{(n)} \longrightarrow \frac{\alpha {{\mathscr {V}}}_r^{(n)} + \beta {{\mathscr {V}}}_r^{(n-1)} }{\alpha +\beta }; \ n=1,2,... \end{aligned}$$
(18)

The input to our calculations is a 3D form of CD-Bonn potential in momentum helicity representation obtained from the summation of partial wave matrix elements of the potential up to total angular momentum \(j_{max} = 20\). For the discretization of the continuous momentum and angle variables, we use the Gauss-Legendre quadrature. A combination of hyperbolic and linear mapping with 120 mesh points is used for the momentum variables, and for azimuthal and polar angle variables, a linear mapping with 40 mesh points is used. In our calculations, we use the nucleon mass \(m = \frac{2 m_p m_n}{ m_p + m_n} = 938.91852\) MeV and the conversion factor \(\hbar c = 197.327053 \ {\text {MeV fm}} = 1\), where proton and neutron masses are \(m_p = 938.27231\) MeV and \(m_n = 939.56563\) MeV12. The number of iterations needed to reach convergence in the matrix elements of relativistic potential calculations in different spin–isospin channels, for different values of the weight averaging parameters \(\alpha \) and \(\beta \), is shown in Table 1. Our numerical analysis indicates that for \((s=0, \ t=0)\) channel, the convergence can be reached only for \(\alpha =\beta =1\), whereas for \((s=0, \ t=1)\) and \((s=1, \ t=0)\) the fastest convergence can be reached by \(\alpha =2, \ \beta =1\). For \((s=1, \ t=1)\) the fastest convergence can be reached by \(\alpha =3, \ \beta =1\). It should be mentioned that Kamada and Glöckle have used \(\alpha = \beta = 1\) in their calculations to obtain relativistic potentials from AV18, CD-Bonn, and Nijm I, II potentials9. In Figs. 1, 2, 3, and 4, the matrix elements of relativistic potential obtained from CD-Bonn potential in different spin–isospin channels are compared with corresponding nonrelativistic potentials. The differences between nonrelativistic and relativistic potentials are also shown. As we can see, while the nonrelativistic and relativistic potentials show similar structures, the difference between them is significant and is in the same order of magnitude as the potentials.

Table 1 The number of iterations \(N_{iter}\) to reach convergence in the solution of Eqs. (14) and (15) for the calculation of relativistic potential in different spin and isospin channels from CD-Bonn potential, as a function of the weight averaging parameters \(\alpha \) and \(\beta \) defined in Eq. (18).
Full size table
Figure 1
figure 1

The matrix elements of relativistic NN potential (left panel), nonrelativistic NN potential (middle panel), and their differences (right panel), calculated for CD-Bonn potential in \((s=0, \ t=0)\) channel, as a function of 2N relative momenta \(p=p'\) and the cosine of the angle between them x.

Full size image
Figure 2
figure 2

The same as Fig. 1, but for \((s=0, \ t=1)\) channel.

Full size image
Figure 3
figure 3

The same as Fig. 1, but for \((s=1, \ t=0)\) channel.

Full size image
Figure 4
figure 4

The same as Fig. 1, but for \((s=1, \ t=1)\) channel.

Full size image

Numerical tests for the relativistic NN potentials

In this section, we present two numerical tests for NN bound and scattering states, which show the validity of our formalism and the accuracy of calculated relativistic potentials in the 3D scheme.

Deuteron binding energy and wave function

To test the accuracy of calculated relativistic NN potential in the \((s=1,\ t=0)\) channel, we calculate the deuteron binding energy and wave function for both nonrelativistic CD-Bonn and the obtained relativistic potentials. The relativistic form of the homogenous LS equation for describing deuteron binding energy \(E_d = m_d - 2m\) and wave function \(\psi _{M_d}(p)\) is given by the following coupled integral equations for wave function components \(\psi _0\) and \(\psi _1\)

$$\begin{aligned} \psi _{M_{d}}(p) = \frac{2\pi }{ m_d - \omega (p) } \int _{0}^{\infty }dp'\ p'^{2} \left\{ \frac{1}{2} {{\mathscr {V}}}_{r,M_{d} 1 }^{+110}(p,p') \psi _{1}(p') + \frac{1}{4} {{\mathscr {V}}}_{r,M_{d}0}^{+110}(p,p') \psi _{0}(p')\right\} , \end{aligned}$$
(19)

where

$$\begin{aligned} {{\mathscr {V}}}_{r,M_{d}\Lambda '}^{+110}(p,p') = \int _{-1}^{1}dx\ {{\mathscr {V}}}_{r,M_{d}\Lambda '}^{+110}(p,p',x) \ d^{1}_{M_{d}\Lambda '}(x) . \end{aligned}$$
(20)

In the nonrelativistic form, the free propagator is replaced by \((E_d - \frac{p^{2}}{m})^{-1}\), and the relativistic NN potential is replaced by the nonrelativistic potential \(V_{nr,M_{d}\Lambda '}^{+110}(p,p')\)17. In Table 2, we present our numerical results for deuteron binding energies obtained from both CD-Bonn and relativistic potentials. The relative percentage difference of \(0.06 \ \%\) indicates an excellent agreement between relativistic and nonrelativistic deuteron binding energies. In Fig. 5, we show the deuteron wave function components calculated for both relativistic and nonrelativistic CD-Bonn potentials. As we can see, the constructed relativistic potential reproduces the deuteron binding energy and wave function obtained by CD-Bonn potential with high accuracy.

Table 2 Calculated deuteron binding energy for relativistic and nonrelativistic CD-Bonn potentials and their relative percentage difference.
Full size table
Figure 5
figure 5

The absolute value of deuteron wave function components \(\psi _0\) and \(\psi _1\), calculated for relativistic (r) and nonrelativistic (nr) CD-Bonn potentials, as a function of the relative momentum p.

Full size image

np elastic scattering

For the second numerical test, we calculate the differential and total cross-section of np elastic scattering for the relativistic potential constructed from CD-Bonn potential. To study describe the relativistic np elastic scattering in momentum helicity space, the relativistic form of inhomogeneous LS equations for 2N t-matrices in singlet and triplet spin states can be obtained as

$$\begin{aligned} T^{\pi 0 t}_{00}(p,p',x)=\, & {} {{\mathscr {V}}}^{\pi 0 t}_{r, 00}(p,p',x) \\&+\, \frac{1}{4}\int _{0}^{\infty }dp''p''^2\int _{-1}^{+1}dx''\ \frac{{{\mathscr {V}}}^{\pi 0 t, 0}_{r,0 0} (p,p'',x,x'')\ T^{\pi 0 t}_{00}(p'',p',x'')}{\omega (p_{0})-\omega (p'')+i\varepsilon }, \end{aligned}$$
(21)
$$\begin{aligned} T^{\pi 1 t}_{\lambda \lambda '}(p,p',x)=\, & {} {{\mathscr {V}}}^{\pi 1 t}_{r, \lambda \lambda '}(p,p',x) \\&+\, \frac{1}{2} \int _{0}^{\infty }dp''p''^2 \int _{-1}^{+1} dx''\ \frac{{{\mathscr {V}}}^{\pi 1 t, \lambda '}_{r,\lambda 1} (p,p'',x,x'')\ T^{\pi 1 t}_{1\lambda '}(p'',p',x'')}{\omega (p_{0})-\omega (p'')+i\varepsilon } \\&+\, \frac{1}{4}\int _{0}^{\infty } dp'' p''^2 \int _{-1}^{+1} dx'' \ \frac{{{\mathscr {V}}}^{\pi 1 t, \lambda '}_{r,\lambda 0} (p,p'',x,x'')\ T^{\pi 1 t}_{0\lambda '}(p'',p',x'')}{\omega (p_{0})-\omega (p'')+i\varepsilon }. \end{aligned}$$
(22)

The on-shell momentum \(p_0\) for an incident projectile energy \(E_{lab}\) is defined by \(p_0 = \sqrt{ \frac{ m_p^2 E_{lab} (E_{lab} + 2m_n ) }{ (m_p + m_n)^2 + 2 m_p E_{lab} } }\)12. By replacing the relativistic free propagator \((\omega (p_{0})-\omega (p'')+i\varepsilon )^{-1}\) with the nonrelativistic one as \((\frac{p_{0}^{2}}{m} - \frac{p''^{2}}{m} +i\varepsilon )^{-1}\) and relativistic potential \({{\mathscr {V}}}_{r}\) with the nonrelativistic potential \(V_{nr}\), Eqs. (21) and (22) yield the nonrelativistic form of LS equations of Ref.16. The relativistic differential cross section of np elastic scattering as a function of the incident projectile energy is given by

$$\begin{aligned} \frac{d\sigma }{d\Omega } = (2\pi )^{4} \left( \frac{\omega (p_{0})}{4} \right) ^{2}\frac{1}{4} \sum _{m'_{s_{1}}m'_{s_{2}}m_{s_{1}}m_{s_{2}}} \left| T^{phys}_{m'_{s_{1}}m'_{s_{2}}m_{s_{1}}m_{s_{2}}}(p_{0},p_{0},x') \right| ^{2} , \end{aligned}$$
(23)

where the matrix elements of the on-shell physical t-matrices can be obtained from the on-shell relativistic NN t-matrices of Eqs. (21) and (22) as

$$\begin{aligned} T^{phys}_{m'_{s_{1}}m'_{s_{2}}m_{s_{1}}m_{s_{2}}}(p_{0},p_{0},x)=\, & {} \frac{1}{4}\sum _{{\pi s t}}(1-\eta _{\pi }(-)^{s+t})\ C^2 \left( \frac{1}{2}\frac{1}{2}t,m_{t_{1}}m_{t_{2}} \right) \\&\times\, C \left( \frac{1}{2}\frac{1}{2}s,m'_{s_{1}}m'_{s_{2}}\lambda '_{0} \right) \ C \left( \frac{1}{2}\frac{1}{2}s,m_{s_{1}}m_{s_{2}}\lambda _{0} \right) \ \sum _{\lambda '}d^{s}_{\lambda '_{0}\lambda '}(x)\ T^{{\pi s t}}_{\lambda '\lambda _{0}}(p_{0},p_{0},x), \end{aligned}$$
(24)

where \(m_{s_{i}}\) and \(m_{t_{i}}\) are the spin and isospin projection of single nucleons along the quantization \(z-\)axis, and the coefficients C are the Clebsch-Gordan coefficients. In Fig. 6, we show the differential cross sections of np elastic scattering, calculated for relativistic and nonrelativistic CD-Bonn potentials, for the projectile energies \(E_{lab}=50, 96, 143\) and 200 MeV. Total cross-sections can provide a more detailed comparison between relativistic and nonrelativistic potentials. In Table 3, we present our numerical results for total cross-sections of np elastic scattering, obtained by relativistic and nonrelativistic CD-Bonn potentials, as a function of the incident projectile energy \(E_{lab}\). As we can see, the maximum relative percentage difference is less than \(0.01 \ \%\), indicating that the relativistic total cross-sections are in excellent agreement with the corresponding nonrelativistic cross-sections.

Figure 6
figure 6

The np differential cross sections calculated for relativistic (solid magenta lines) and nonrelativistic CD-Bonn (dashed red lines) potentials as a function of scattering angle \(\theta _{c.m.}\) in the CM frame, for the incident projectile energies \(E_{lab}=50, \ 96, \ 143\), and 200 MeV. The experimental data are from Refs.42,43 for \(E_{lab} = 50\) MeV (Exp. 1 and Exp. 2, respectively), form Refs.44,45,46 for \(E_{lab} = 96\) MeV (Exp. 1, Exp. 2, and Exp. 3, respectively), from Ref.47 for \(E_{lab} = 142.8\) MeV (Exp. 1) and from Refs.48,49,50 for \(E_{lab} = 200\) MeV (Exp. 1, Exp. 2, and Exp. 3, respectively).

Full size image
Table 3 The total cross section of np elastic scattering as a function of the incident projectile energy \(E_{lab}\), calculated for relativistic (\(\sigma _{r}\)) and nonrelativistic (\(\sigma _{nr}\)) CD-Bonn potentials.
Full size table

Kamada and Glöckle have shown in Ref.9 that the obtained relativistic potential from AV18 potential, in a PW decomposition, reproduces the nonrelativistic phase shifts with five significant figures with projectile energy in the domain \((1-350)\) MeV. As one can see in Table 3, our nonrelativistic and relativistic cross-section results are also in perfect agreement with five significant figures with an incident projectile energy in the broader domain (0.001–750) MeV. So we are convinced that the 3D formulation and calculations for relativistic potential provide the same accuracy as a PW calculation.

Moreover, in a prior study for calculating relativistic potentials from spin-independent MT potential10, which has no spin and isospin complexity of CD-Bonn, we obtained a relative percentage difference of \(0.06\%\) between nonrelativistic and relativistic deuteron binding energies and a maximum relative percentage difference of \(0.007\%\) in two-body total elastic scattering cross-sections, which can be compared with \(0.06\%\) and \(0.01\%\) relative percentage differences obtained in this study for CD-Bonn potential. This comparison indicates calculating relativistic NN interactions from realistic interactions in a 3D scheme provides almost the same accuracy as a spin-independent calculation.

Conclusion and outlook

In this paper, the quadratic equation, which connects the relativistic and nonrelativistic NN potentials, is formulated in momentum helicity space as a single and four coupled three-dimensional integral equations for 2N singlet and triplet spin states. In our numerical calculations, we implement the CD-Bonn potential to obtain the matrix elements of the relativistic potential as a function of the magnitude of 2N relative momenta, the angle between them, and spin and isospin quantum numbers. The quadratic integral equations are solved using an iterative scheme. Our numerical results indicate that calculated relativistic NN potential from the CD-Bonn potential reproduces 2N observables for deuteron binding energy and the differential and total cross sections of np elastic scattering with high accuracy. The implementation of relativistic NN potentials in the relativistic description of triton binding energy and wave function is currently underway.