Perturbative calculations of nucleon-deuteron elastic scattering in chiral effective field theory

The nucleon-deuteron ($ Nd $) system provides an important testing ground for chiral nuclear forces. Not only does it test the prediction of chiral nucleon-nucleon ($ NN $) potentials for the three-nucleon ($ 3N $) system [1−7], but it also helps in understanding the role of $ 3N $ potentials by quantifying their importance via power counting [8−15]. While the triton bound state and its properties have been the natural choice, the scope of investigation is limited by the quantum numbers of the triton, e.g., the total angular momentum $ J = \frac{1}{2} $. Nucleon-deuteron scattering, especially neutron-deuteron ($ nd $) scattering, offers more probes of the $ 3N $ system without having to account for precision-level details of electromagnetic or weak interactions. In recent years, strict perturbative treatment of subleading-order interactions has been increasingly advocated in the development of effective field theories (EFTs) for nuclear physics. Perturbation theory on top of a nonperturbative leading order (LO) has the advantage of disentangling subleading interactions from the LO ones. For instance, two-pion exchange (TPE) potentials in $ NN $, although subleading, can become much stronger than the one-pion exchange (OPE) potential at intermediate momenta near the ultraviolet cutoff $ \sim \Lambda $. Perturbative calculations make it clear how TPEs are subtracted in the ultraviolet region by subleading contact interactions, thereby producing small corrections to on-shell amplitudes [16−19]. However, a perturbative treatment complicates the computation when the $ 3N $ continuum problem is already more involved than the bound-state problem. The main focus of this paper is to develop a technique to perform these perturbative calculations for chiral nuclear forces in the context of nucleon-deuteron scattering.

For the nuclear force, we use the power counting developed in Refs. [17, 18], and later modified by Ref. [20], to organize chiral $ NN $ forces at different orders. More specifically, we follow Refs. [17, 18] to treat the so-called nonperturbative-pion channels: $ {{}^{{{1}}} {{S}}_{{{0}}}} $, $ {{{}^{{{3}}} {{S}}_{{{1}}}}-{{}^{{{3}}} {{D}}_{{{1}}}}} $, and $ {{}^{{{3}}} {{P}}_{{{0}}}} $. For other partial waves where OPE is considered as a perturbation, the power counting laid out in Ref. [20] is adopted. This scheme was explained recently in Ref. [21], and the part relevant to this paper is reviewed in more detail in Sec. 3. The same power counting is also employed to study electroweak processes in Refs. [21−23]. We note that the power counting of Refs. [17, 18] has been examined with Bayesian analyses for $ NN $ scattering data [24, 25] and has been used to study the structure of the deuteron and triton [26]. In constructing the power counting for two-body potentials, these works adopt renormalization-group (RG) invariance as a guideline, which requires that the phase shifts be independent of the momentum cutoff, an arbitrary parameter of the ultraviolet regularization.

The OPE potential is the most important long-range nucleon force in chiral EFT, and it has a tensor component that behaves like $ 1/r^3 $ at short distances. For attractive singular potentials, such as the OPE tensor force in $ {{{}^{{{3}}} {{S}}_{{{1}}}}-{{}^{{{3}}} {{D}}_{{{1}}}}} $, $ {{}^{{{3}}} {{P}}_{{{0}}}} $, and $ {{{}^{{{3}}} {{P}}_{{{2}}}}-{{}^{{{3}}} {{F}}_{{{2}}}}} $, RG invariance requires a contact potential, often referred to as a counterterm, to appear at LO if OPE is considered nonperturbative in that partial wave [27−30], even though the naive dimensional analysis (NDA) adopted by Weinberg's power counting would not require one [31−33]. Because there are, in principle, an infinite number of attractive singular channels for OPE, one would have to invoke an infinite number of counterterms already at LO. This conundrum is avoided once we recognize that OPE does not need to be resummed nonperturbatively in the Lippmann-Schwinger or Schrödinger equation for sufficiently high orbital angular momentum [20, 34, 35] and that NDA is restored if OPE is treated in pure perturbation theory. We follow Ref. [20] in letting OPE enter at LO only in $ {{}^{{{1}}} {{S}}_{{{0}}}} $, $ {{{}^{{{3}}} {{S}}_{{{1}}}}-{{}^{{{3}}} {{D}}_{{{1}}}}} $, and $ {{}^{{{3}}} {{P}}_{{{0}}}} $ and at next-to-leading order (NLO) in all other waves. Not only does this development of two-body chiral forces serve as the foundation of our study of nucleon-deuteron scattering, but it also illustrates the intertwined logic among renormalization, power counting, and perturbation theory for subleading interactions.

A universal feature of EFTs is the increasing momentum power of higher-order interactions. Although this facilitates expansions for low-momentum initial and final states where momenta Q are well below the breakdown scale $ M_{\text{hi}} $, these higher-order interactions are not necessarily small for intermediate states with momenta up to the ultraviolet cutoff $ \Lambda \gtrsim M_{\text{hi}} $. Perturbative renormalization of subleading orders has been advocated as a reliable way to ensure that the resulting large contributions from intermediate states can be absorbed into low-energy constants (LECs). Numerous applications of strict perturbation theory in pionless and chiral EFTs can be found in the recent review in Ref. [36].

A key aspect of the power counting of chiral forces used in this paper is that the LO potentials are nonzero only in a limited number of $ NN $ partial waves. In the distorted-wave expansion, perturbation theory for subleading potentials is applied by directly evaluating matrix elements between the LO asymptotic states. By contrast, our technique solves a hierarchy of integral equations at subleading orders, all of which share the kernel from the LO equation but have a distinct driving term at each order. This approach to implementing perturbation theory for subleading-order interactions is in line with the methods developed in Refs. [37, 38] for pionless-EFT calculations of few-body systems. A similar framework has been developed in Ref. [21] to calculate the longitudinal response function of the deuteron perturbatively.

This paper is organized as follows: in Sec. 2, we describe the numerical framework for solving the Faddeev equation, focusing on the contour-deformation method. In Sec. 3, we give details of the perturbative treatment of the NLO potentials. We then present benchmark calculations to validate our methods in Sec. 4. The LO and NLO results for $ Nd $ elastic scattering are presented and discussed in Sec. 5, and we conclude with a summary in Sec. 6.

In our calculations, we expand the Faddeev equation in the Jacobi partial-wave basis [39]:

Here, p and q are the magnitudes of the Jacobi momenta: $ \vec{p} $ is the relative momentum of the subsystem (nucleons 1 and 2), and $ \vec{q} $ is the momentum of the "spectator" (nucleon 3) relative to the center of mass of the subsystem. The quantum numbers l, s, j, and t are, respectively, the orbital angular momentum, spin, total angular momentum, and isospin of the subsystem; λ and I are the orbital angular momentum and total angular momentum of the spectator; and J and T are the total angular momentum and total isospin of the $ 3N $ system. The Jacobi partial-wave basis is partially antisymmetrized for the $ NN $ subsystem, i.e., $ l+s+t = $ odd, and the parity is given by $ P=(-1)^{l+\lambda} $. This paper focuses on $ Nd $ elastic scattering, for which $ T=\dfrac{1}{2} $. For simplicity, we use the collective label α to denote these discrete quantum numbers. Our choice for the normalization of the Jacobi partial-wave basis is as follows:

We require the initial and final wave functions, $ \left| {\phi} \right\rangle $, describing configurations in which the nucleon and the deuteron are far apart. To project an initial or final wave function onto the partial-wave basis, we first enumerate all $ 3N $ channels that include a deuteron $ NN $ channel:

The orbital angular momentum is restricted to $ 0 $ or $ 2 $ (i.e., $ l_d = 0 $ or $ 2 $). For a given value of the total $ 3N $ angular momentum J, there may be multiple combinations of λ and I that yield a valid $ \alpha_d $ channel. The initial $ Nd $ state, with definite J, λ, and I, is then constructed as

where $ \varphi_{l_d}(p) $ denotes the deuteron wave function in the $ NN $ partial wave $ l_d $, and $ q_0 $ is the center-of-mass momentum of the incoming nucleon.

However, it is more customary to use the $ \lambda \Sigma $ basis for defining the phase shifts and mixing angles in $ Nd $ elastic scattering, where Σ is the channel spin, i.e., the total spin of the deuteron and the incoming or outgoing nucleon [40]:

The $ \lambda\Sigma $ and $ \lambda I $ bases are related as follows:

For nucleon-deuteron scattering, the Faddeev equation can be greatly simplified because the three nucleons are identical fermions. Particle-exchange operators are crucial for enforcing the fermionic nature of the nucleons. If $ P_{ij} $ denotes the exchange of the two nucleons labeled i and j, the $ 3N $ permutation operator P is the sum of the cyclic and anticyclic permutations of the three nucleons [39]:

In practice, the permutation operator P is projected onto the Jacobi partial-wave basis (1):

where

and

Here, $ P_k $ denotes the Legendre polynomial of order k; $ l_{1, 2} $ and $ l_{1, 2}' $ are intermediate angular momenta that are summed over; the final-state quantum numbers, such as $ l', \lambda', s', \cdots $, are collectively denoted by $ \alpha' $; and the coefficient $ g_{\alpha'\alpha}^{kl_1'l_2'l_1l_2} $ is defined as follows [2]:

where $ \hat{X}\equiv{2X+1} $; L and S denote, respectively, the total orbital angular momentum and the total spin allowed by the total angular momentum J; and f and $ f' $ are, again, the intermediate angular momenta used to facilitate recoupling.

The inhomogeneous Faddeev equation is diagrammatically illustrated in Fig. 1. We use chiral potentials in this study, and the $ 3N $ forces do not contribute up to NLO for renormalization purposes, as shown in Ref. [41] and verified in this work. Therefore, only two-body potentials are considered here. The yellow blob denotes the breakup amplitude T, which starts from an $ Nd $ initial state ϕ and ends with three free nucleons in the final state. The solid circle represents the full off-shell two-body t-matrix, which satisfies the Lippmann-Schwinger equation (LSE).

Figure 1.  (color online) Diagrammatic representation of the Faddeev equation. Solid lines represent nucleons; the solid circle denotes the two-body off-shell t-matrix; the blue half-circle denotes the deuteron; and the yellow blob denotes the Faddeev breakup amplitude T. Not all particle-exchange topologies are shown.

where $ V_2 $ is the two-body potential and $ G_0 $ is the free propagator. The propagation of the breakup process, denoted by $ T\left| {\phi} \right\rangle $, is given symbolically by the following equation [39, 2]:

The total energy of the $ 3N $ system, $ E_3 $, is related to the center-of-mass momentum $ q_0 $ of the incoming nucleon by:

where $ E_d = - B_d $ denotes the (negative) deuteron binding energy and $ m_N $ is the nucleon mass.

The Faddeev equation (14) is solved in the $ \lambda I $ basis to obtain the Faddeev breakup amplitude $ \left( pq\alpha \right) T\left|{\phi_{\lambda I}^J} \right\rangle $. In turn, the $ Nd $ elastic amplitude U is computed from the following relation [2]:

The basis transformation from $ \lambda I $ to $ \lambda \Sigma $ is carried out according to Eq. (6). The S-matrix for elastic scattering is directly related to $ U^J_{\lambda'\Sigma', \lambda \Sigma} $ by

We follow Ref. [40] in parameterizing the S-matrix, thereby defining the phase shifts and mixing angles.

Projecting the abstract operator equations (14), (13), and (16) onto the partial-wave basis yields integral equations for $ NN $ and $ 3N $ dynamics. The momentum-space integrals in these equations typically involve singularities. In numerical computations, we employ contour deformation to circumvent these singularities. Suppose $ p'' $ is the integration variable and that the original contour runs from $ 0 $ to $ \infty $ along the positive real axis; we deform the contour by rotating it clockwise by a small angle θ:

In Fig. 2, the rotated contour is illustrated by a ray at angle θ. For sufficiently large $ p'' $, the deformed contour returns counterclockwise to $ +\infty $ on the real axis. In practical calculations, we ensure that the integrand decays rapidly enough that the integral along the arc can be neglected. This technique is akin to the "complex scaling" method used in many coordinate-space and momentum-space calculations of scattering and reaction processes, e.g., in Refs. [42−47].

Figure 2.  (color online) Diagram illustrating the analytic structure of $ V^{1\pi}(p',p'') $ when the integration contours in both $ p' $ and $ p'' $ are rotated by an angle θ. The red lines trace the trajectories of the two branch points at $ p'' = p' \pm i m_{\pi} $, which arise from the endpoint singularities of the integral.

The implementation begins with the LSE (13):

where p ($ p' $) denotes the incoming (outgoing) relative momentum, and $ E_2 $ is the center-of-mass energy of the $ NN $ pair. We introduce the notation for coupled $ NN $ partial waves: $ \bar{l} = l $ for uncoupled channels, while for coupled channels $ \bar{l} $ takes the values $ j-1 $ and $ j+1 $. Not only is the $ p'' $ contour deformed by rotation, but the arguments of the full off-shell t-matrix $ t(p', p, E_2) $ — $ p' $ and p — also lie along the rotated axis. The most prominent singularity is the pole of the free $ NN $ propagator at $ p'' = \sqrt{m_NE_2} $. The deformed contour clearly avoids it. In our numerical integrations, these "rotated" complex variables are represented by real Gauss–Legendre mesh points multiplied by a complex phase.

Thus, the integral is approximated by a sum:

where $ {x_n} $ and $ {w_n} $ are the standard abscissae and weights on the positive real axis.

In the ChEFT construction of $ V_{l'\bar{l}'}(p', p'') $, the contact potentials are usually polynomials in $ p' $ and $ p'' $; therefore, they do not introduce any singularities in the contour integration. However, the pion-exchange components of $ V_{l'\bar{l}'}(p', p'') $ exhibit nontrivial singularities. The OPE potential $ V^{1\pi} $ admits the following integral representation in momentum space:

where $ m_{\pi} $ is the pion mass. For fixed $ p' $, $ V^{1\pi}_{l'\bar{l}'}(p', p'') $ has branch points in the complex $ p'' $-plane at $ \pm p' \pm \rm{i} m_\pi $, arising from the endpoint singularities of the integral, as illustrated in Fig. (2). As $ p' $ varies along the rotated axis, these branch points trace out boundaries that the $ p'' $ contour cannot cross, as indicated by the solid red lines in the figure. This configuration does not pose a problem because the $ p'' $ contour runs parallel to these boundaries.

To regularize the ultraviolet behavior of the potential, we introduce the following separable regulator:

where Λ denotes the ultraviolet momentum cutoff. For large values of p and $ p' $, the regulator takes the following asymptotic form:

The chosen value of θ is usually small enough that $ \cos{4\theta} \,(p^4/\Lambda^4) $ remains positive, which in turn ensures the proper ultraviolet regularization of $ t(p', p, E_2) $.

Although we do not need the on-shell t-matrix $ t(k, k, E_2) $ in this paper, where $ k = \sqrt{m_N E_2} $ is real, it can be calculated from the off-shell solution by using the known off-shell $ t(p'', p, E_2) $ as input on the right-hand side of Eq. (19). In this contour integration, the branch point $ (k, - \mathrm{i} m_\pi) $ does not cross the deformed contour provided that $ k< m_{\pi}/\tan{\theta} $, as illustrated in Fig. 3. This condition imposes an upper limit on the accessible values of k, which is nevertheless sufficiently high for ChEFT applications, where the on-shell momenta under consideration are typically $ \lesssim 3 m_\pi $.

Figure 3.  The analytic structure of $ V^{1\pi}(k,p'') $ with the $ p'' $ contour rotated by an angle θ. The solid dots indicate the end-point singularities of the integral.

The Faddeev equation in momentum space is given by:

Here, $ \bar{\alpha} $ denotes the channel coupled to α via the $ NN $ interaction, obtained by replacing l in α with $ \bar{l} $. ϕ is the collective symbol for the quantum numbers of the initial state, including $ q_0 $. The integral on the right-hand side of Eq. (25) exhibits three types of singularities that must be handled if all momenta are kept on the real axis, as detailed in Ref. [2]. First, the free propagator $ G_0 $ in the $ q'' $ integral has a pole when $ E_3 $ exceeds the $ 3N $ threshold. Second, $ t_{l'\bar{l}'}(p', \pi_1, z) $ possesses a branch point at $ z=0 $, corresponding to the $ q' $ satisfying $ E_3 - 3q'^2/(4m_N) = 0 $. This induces a singularity of $ T(p'q'\alpha'; \phi) $ as a function of $ q' $, which in turn produces a singularity of $ T(\pi_2\, q'' \alpha''; \phi) $ in the $ q'' $ integral. Third, in the $ {}^3S_1-{}^3D_1 $ channel, $ t_{l'\bar{l}'}(p', \pi_1, z) $ has an additional singularity: the deuteron pole at $ z = E_d $. If $ E_3 $ exceeds the nucleon-deuteron threshold, the deuteron pole corresponds to a singularity in $ q' $ where $ E_3 - 3q'^2/(4m_N) = E_d $, which again creates a singularity of $ T(\pi_2\, q'' \alpha''; \phi) $ in the $ q'' $ integral.

In our implementation, these singularities are avoided by rotating $ p' $, $ q' $, and $ q'' $ in Eq. (25) into the complex plane, as defined by Eq. (18), while the on-shell momentum $ q_0 $, the energies $ E_d $ and $ E_3 $, and the Legendre variable x remain real-valued. The argument $ \pi_2(q', q'') $ of $ T(\pi_2 q'' \alpha''; \phi) $ also lies on the rotated contour. Consequently, only one spline operation is required to map $ \pi_2(q', q'') $ onto the $ p'' $ mesh:

where $ p''_m $ denote the rotated discrete momenta defined by Eq. (20), and the spline functions $ S_m $ are given in Ref. [48]. Using the momentum mesh and spline functions transforms the integral equation (25) into a linear system, which we denote symbolically by

Here, the indices k, r, m, and n range over the momentum-mesh points associated with $ p' $, $ q' $, $ p'' $, and $ q'' $, respectively. The kernel K, the unknown vector Ψ, and the driving term D are defined as follows:

The Faddeev breakup amplitude T is used to compute the elastic scattering amplitude U. This is obtained through an additional integration:

$ U_{\lambda ' I', \lambda I}^{J}(q_0) $ is then transformed into $ U_{\lambda ' \Sigma', \lambda \Sigma}^{J}(q_0) $ via Eq. (6). The integration contour for $ q'' $ in Eq. (31) is also rotated; consequently, the argument $ \pi_2(q_0, q'') $ of T in the second line does not lie on the rotated contour for $ p' $. To evaluate $ T[\pi_2(q_0,q'')q''\alpha'';\phi_{\lambda I}^{J}] $, which is then fed into Eq. (31), we again use Eq. (25) by setting $ p' = \pi_2(q_0,q'') $. Therefore, two additional integration steps are required to compute $ U_{\lambda ' I', \lambda I}^{J}(q_0) $ after obtaining the initial solution $ T(p'_n q'_m \alpha'; \phi) $ via Eq. (25), where $ p'_n $ and $ q'_m $ are the rotated meshes defined by Eq. (20).

A key feature of our Faddeev-equation implementation is the perturbative treatment of subleading-order EFT interactions. We adopt the power counting of chiral nuclear forces as presented in Refs. [17, 18, 20]. Two types of soft scales arise in EFTs. The first comprises nucleon momenta, e.g., the initial and final momenta of the nucleon and deuteron, and the deuteron binding momentum $ \gamma_d \equiv \sqrt{m_N B_d} \simeq 46 $ MeV. The second consists of dimensionful parameters encoded in LECs, such as the $ NN $ scattering lengths and effective ranges. We use the pion mass $ m_\pi $ as a generic proxy for these scales. The breakdown scale is chosen to be the nucleon-delta mass splitting $ \delta = 293 $ MeV.

At LO, where the potential must be treated nonperturbatively — that is, where the Faddeev equation, Eq. (25), is solved exactly — we include OPE potentials in the $ {{}^{{{1}}} {{S}}_{{{0}}}} $, $ {{{}^{{{3}}} {{S}}_{{{1}}}}-{{}^{{{3}}} {{D}}_{{{1}}}}} $, and $ {{}^{{{3}}} {{P}}_{{{0}}}} $ partial waves, where contact potentials provide the necessary short-range interactions alongside the OPE. In all other partial waves, OPE is sufficiently weak to be treated as an NLO correction. For $ {{}^{{{1}}} {{S}}_{{{0}}}} $ and $ {{{}^{{{3}}} {{S}}_{{{1}}}}-{{}^{{{3}}} {{D}}_{{{1}}}}} $, the contact potential is expressed in momentum space as

where the superscript "(0)" indicates LO, "(1)" denotes NLO, and so on. We suppress the channel label on the LECs when there is no risk of confusion. At $ {{}^{{{3}}} {{P}}_{{{0}}}} $, the contact interaction takes the following form:

At NLO, OPE begins to contribute to additional channels, up to a maximum orbital angular momentum of $ l \leqslant 2 $. We also retain the higher partial waves coupled to these channels; accordingly, $ {{{}^{{{3}}} {{P}}_{{{2}}}}-{{}^{{{3}}} {{F}}_{{{2}}}}} $ and $ {{{}^{{{3}}} {{D}}_{{{3}}}}-{}^{{3}}{G}_{{3}}} $ are included. In addition, the $ {{}^{{{1}}} {{S}}_{{{0}}}} $ potential receives the following NLO correction [18]:

where $ C_0^{(1)} $ is the NLO correction to the LO LEC $ C_0^{(0)} $, and $ D_0^{(1)} $ is the LEC associated with the momentum-dependent $ {{}^{{{1}}} {{S}}_{{{0}}}} $ contact interaction.

In summary, our chiral potentials contain three undetermined LECs at LO and one at NLO. These LECs are fixed by fits to the phase shifts from the Nijmegen partial-wave analysis [49, 50] up to $ k=300 $ MeV, with k the center-of-mass momentum. In addition, the deuteron binding energy, $ B_d=2.225 $ MeV, is reproduced at LO and NLO.

A crucial feature of this power counting for implementing the Faddeev equation is the significantly smaller number of channels at LO compared to subleading orders, as shown in Table 1. We leverage this feature to implement a perturbative treatment of the subleading interactions. To this end, the full space of $ 3N $ channels, $ \mathscr{C} $, is decomposed into two sets. The first, $ \mathscr{A} $, consists of Jacobi partial waves whose two-body subsystem is subject to the LO interactions $ {{}^{{{1}}} {{S}}_{{{0}}}} $, $ {{{}^{{{3}}} {{S}}_{{{1}}}}-{{}^{{{3}}} {{D}}_{{{1}}}}} $, and $ {{}^{{{3}}} {{P}}_{{{0}}}} $. The remaining channels form the complementary set $ \mathscr{B} $:

For any channel β in $ \mathscr{B} $ and any other 3N channel γ, we have

As we shall see shortly, this seemingly trivial identity proves useful for perturbative calculations.

To treat the NLO potential perturbatively, we begin with the formal EFT expansions of the two-body potential V, the two-body t-matrix t, the three-body Faddeev breakup operator T, and the initial state ϕ:

Substituting these expansions into Eqs. (13) and (14) yields the following perturbative hierarchy:

where the $ t^{(1)} $ is given by

This result should be compared with the direct calculation of $ T^{(1)} $.

Equations (41), (42), and (43) share the same integral kernel.

but differ in their driving terms,

By Eq. (36), the operator $ t^{(0)}PG_0 $, acting on the right, annihilates any $ 3N $-channel state $ \beta \in \mathscr{B} $:

As a result, the kernel has a block-triangular structure:

where $ K_A $ acts only within $ \mathscr{A} $ and $ K_{AB} $ couples $ \mathscr{A} $ to $ \mathscr{B} $.

Here, and later when it does not cause confusion, we omit the labels for the momentum mesh. Upon discretization, $ K_A $ is typically smaller than the full K, thus saving computational resources. Accordingly, the unknown vector Ψ and the driving term D are also split into two parts:

where the unknown vectors at each order are defined as follows:

With this decomposition, we can write Eqs. (41) and (42) as:

and

The case corresponding to Eq. (43) is more involved.

where for $ \alpha' \in \mathscr{A} $,

Therefore, even at NLO we deal with a linear system that spans only $ \mathscr{A} $. While the integral kernel remains identical to that at LO, the main computational effort at NLO shifts to constructing the driving term. By mathematical induction, one finds that at all higher orders the same kernel is reused, although the driving terms become increasingly complex. Because of this feature, we refer to this approach to perturbative calculations as fixed-kernel perturbation theory (FKPT).

In all implementations in this work, $ \mathscr{A} $ comprises 6 channels for $ J=\frac{1}{2} $ and 8 for $ J\geqslant\frac{3}{2} $. The dimension of the full space of $ 3N $ channels $ \mathscr{C} $ varies with J as follows: 22 for $ J=\frac{1}{2} $, 38 for $ J=\frac{3}{2} $, 46 for $ J=\frac{5}{2} $, $ \cdots $. In general, the dimension of the LO kernel is reduced by a factor of $ \left(\frac{\text{Dim}_{\mathscr{A}}}{\text{Dim}_{\mathscr{C}}}\right)^2 $ compared to that of the full kernel; consequently, the former is at least an order of magnitude smaller than the latter.

Using $ T^{(0)}\left| {\phi^{(0)}} \right\rangle $, $ T^{(0)}\left| {\phi^{(1)}} \right\rangle $, and $ T^{(1)}\left| {\phi^{(0)}} \right\rangle $ as inputs, we can express the elastic scattering amplitude at LO and NLO as

In the chiral power counting adopted in this work, the $ {{{}^{{{3}}} {{S}}_{{{1}}}}-{{}^{{{3}}} {{D}}_{{{1}}}}} $ interaction vanishes at NLO, as noted in Table 1. As a result, the NLO correction to the deuteron wave function is also zero, which implies that the NLO corrections to the initial and final $ Nd $ states vanish as well:

Although exact unitarity is violated in perturbation theory, we can still extract the NLO phase shifts from $ U^{(1)} $ in a manner consistent with power counting. As detailed for $ NN $ coupled channels (see, e.g., Ref. [17]), the phase-shift parameters are expanded as

where n indexes the various phase-shift parameters, including mixing angles. The exact relation between U and $ \delta_n $ is provided in Ref. [40] and can be written symbolically as follows:

We expand both sides:

Using $ U^{(1)} $ from our perturbative calculations as input, together with the matrix $ \partial W/\partial \delta_n $, we can extract $ \delta_n^{(1)} $ by matching orders on both sides.

To determine the number of one-dimensional momentum mesh points, N, and the contour-rotation angle, θ, we compare several options. Table 2 and Table 3 show the effects of varying N and θ, respectively, on the $ ^2S_{\frac{1}{2}} $ phase shifts. The variations are generally at the subpercent level or smaller, which is more than sufficient for our purposes. For definiteness, we adopt $ N=48 $ mesh points and a contour-rotation angle of $ 10^\circ $.

To benchmark our contour-deformation implementation of the Faddeev equation, we compare our results with those obtained using the wave-packet continuum-discretization (WPCD) method [51]. We also note other studies in which WPCD has been successfully combined with various interactions to compute $ Nd $ scattering [52, 10, 53]. For the interaction, we use LO potentials with a cutoff $ \Lambda=400 $ MeV.

We calculate phase shifts and mixing angles at nucleon laboratory energies $ E_N = 3 $, $ 14 $, and $ 30 $ MeV for $ J^P = \dfrac{1}{2}^\pm $ and $ \dfrac{3}{2}^+ $. The results from the two methods are compared in Table 4 for $ J^P=\dfrac{1}{2}^\pm $ and in Table 5 for $ J^P=\dfrac{3}{2}^+ $. Overall, the discrepancies between the results of the two methods are less than one percent. This level of agreement is sufficient for the present study, as the subsequent NLO corrections and cutoff variations are expected to be substantially larger than the differences introduced by the numerical methods.

We also compare the differential cross sections and nucleon analyzing powers, $ A_y $, calculated using the two methods, as shown in Fig. 4. Following Ref. [2], we relate the partial-wave amplitude $ U^J $ to these observables. For these observables, the partial-wave sum of the scattering amplitude $ U^J $ includes contributions up to $ J^{P} \leqslant \frac{15}{2}^{\pm} $. The level of agreement for these observables is consistent with that found for the phase shifts in Tables 4 and 5.

Figure 4.  (color online) Comparison of the $ Nd $ elastic-scattering differential cross sections and nucleon analyzing powers. Results obtained with the contour-deformation and WPCD methods are shown as solid and dashed curves, respectively.

Implementing the FKPT method — our perturbative treatment of NLO potentials outlined in Sec. 3 — requires substantial extensions to the code beyond the standard Faddeev-equation solver. To validate this functionality, we note that perturbation theory can be implemented straightforwardly within the nonperturbative framework, albeit with additional computational cost. Starting from an auxiliary potential,

Here, x is an auxiliary multiplier sampled over an interval around 0. We can use the nonperturbative solver, i.e., Eqs. (25) and (31), to calculate the x-dependent scattering amplitude $ U(E_3, x) $ for arbitrary x. By numerically expanding $ U(E_3, x) $ about $ x = 0 $, we obtain the NLO scattering amplitude in perturbation theory:

This approach is computationally demanding, as it requires multiple evaluations of the amplitude at different values of x to carry out the Taylor expansion numerically. In practice, we use at least five x values to ensure the numerical stability of the NLO results. Compared with the FKPT method, the auxiliary-potential method requires at least an order of magnitude more computation time when both methods use the same number of mesh points N and the same set of truncated $ 3N $ channels $ \left| {\alpha} \right\rangle $. Furthermore, the auxiliary-potential method has a large memory footprint because it does not exploit the significantly smaller size of the LO kernel. While useful for benchmarking, these limitations motivate the adoption of the FKPT approach. For the sampled phase shifts, the two approaches yield consistent results, agreeing to at least six significant digits.

We use $ Nd $ elastic scattering to investigate two aspects of the power counting of chiral nuclear forces adopted in this paper. The LO and NLO potentials are explained in Sec. 3. One is the ultraviolet-cutoff dependence of the phase shifts, and the other is the effects of the NLO correction.

The ultraviolet cutoff Λ introduced in Eq. (23) is arbitrarily chosen, and observables must be independent of its value up to the EFT uncertainty allowed at the given order; this is a manifestation of RG invariance. A violation of RG invariance can be interpreted as a defect in the power counting, which can be remedied by promoting certain operators to lower orders. In Ref. [41], RG invariance of the triton binding energy and $ Nd $ scattering lengths was observed, leading to the conclusion that the $ 3N $ forces are not required up to NLO for renormalization purposes.

We verify this conclusion by inspecting how the phase shifts vary with Λ. Because ultraviolet divergences tend to be suppressed at higher orbital angular momentum, it suffices to study the S-wave phase shifts of $ Nd $ scattering. When Λ exceeds $ \sim 600 $ MeV in $ {{}^{{{3}}} {{P}}_{{{0}}}} $ and $ \sim 1000 $ MeV in $ {{{}^{{{3}}} {{S}}_{{{1}}}}-{{}^{{{3}}} {{D}}_{{{1}}}}} $, spurious $ NN $ bound states begin to develop [29]. In order to avoid unphysical breakup thresholds in the $ Nd $ system, we remove these spurious states using a method similar to that described in Ref. [29]. More concretely, we raise the energies of the spurious $ NN $ states to a large positive value, for example a few GeV, so that they lie far outside the EFT region and do not open undesired reaction channels. This is implemented by adding the following to the LO two-body potential:

where $ \left| {\Psi_s} \right\rangle $ denotes the normalized spurious state, and Ω denotes the artificial positive energy assigned to it. Figure 5 presents the phase shifts of $ ^2S_{\frac{1}{2}} $ and $ ^4S_{\frac{3}{2}} $ as Λ varies from 400 to 1600 MeV. Up to $ E_N = 108 $ MeV, the phase shifts exhibit convergence with respect to Λ.

Figure 5.  (color online) The doublet- and quartet-channel S-wave phase shifts as functions of $ E_N $ for various values of the cutoff Λ. The top and bottom panels correspond to the $ ^2S_{\frac{1}{2}} $ and $ ^4S_{\frac{3}{2}} $ channels, respectively. The "Bonn-B" results are taken from Ref. [54].

The angular distributions of the differential cross sections and the nucleon analyzing power $ A_y $ are computed to study how much the NLO corrections change relative to LO. These observables are computed with the partial-wave sum truncated beyond $ J^P=\dfrac{19}{2}^{\pm} $, at which point including higher waves alters the results by less than $ 1\% $. Using $ \Lambda = 400-1600 $ MeV, Figs. 6 and 7 compare the LO and NLO EFT predictions with available experimental data for these observables at their respective energies.

Figure 6.  (color online) The LO and NLO angular distributions at various $ E_N $ are shown. The bands are generated by varying Λ from 400 to 1600 MeV. The $ nd $ scattering data are taken from Ref. [55].

Figure 7.  LO and NLO nucleon analyzing powers at various $ E_N $ are shown. The bands are generated by varying Λ from 400 to 1600 MeV. The $ nd $ scattering data at $ E_N=3.0 $, 10.0 and 14.1, 21.0, and 67.0 MeV are taken from Refs. [56], [57] (for 10.0 and 14.1 MeV), [58], and [59]. The $ pd $ data at $ E_N=35.0 $ MeV are from Ref. [60].

For the angular distribution, we find that the LO results agree with the data better than the NLO ones at forward angles at most energies studied. We first note that the NLO correction is mostly driven by repulsion in the $ {{}^{{{3}}} {{P}}_{{{1}}}} $ channel, which consists of only the OPE potential at the NLO. As shown in Fig. 8, removing the $ {{}^{{{3}}} {{P}}_{{{1}}}} $ component from the NLO potentials results in the angular distributions nearly identical to the LO results. Because the tensor operator $ S_{12} $ has its largest matrix element in $ {{}^{{{3}}} {{P}}_{{{1}}}} $, compared with other perturbative channels, the OPE tensor force in $ {{}^{{{3}}} {{P}}_{{{1}}}} $ appears to be the strongest at NLO. However, it is not entirely clear to us whether other mechanisms exist for the $ {{}^{{{3}}} {{P}}_{{{1}}}} $ OPE to dominate the NLO correction to the differential cross section. Using the nucleon-delta mass splitting $ \delta = 293 $ MeV as the breakdown scale for our delta-less chiral forces, we estimate the EFT expansion error by powers of $ \mathrm{max}(q_0, \gamma_d, m_\pi)/\delta $. For instance, the EFT error for NLO is $ (m_{\pi}/\delta)^2\simeq 22\% $ at $ E_N = 9.0 $ MeV, which is able to explain the discrepancy between the EFT predictions and the data.

Figure 8.  (color online) The LO and NLO angular distributions, obtained using the chiral $ NN $ potential with $ \Lambda=800 $ MeV, are shown for various $ E_N $. The dotted curves show the NLO results with the $ {{}^{{{3}}} {{P}}_{{{1}}}} $ OPE turned off.

The LO, however, agrees worse than the NLO with data in terms of describing $ A_y $. In Fig. 7, the maximum of $ A_y $ at LO has the opposite sign relative the experimental data. Due to the lack of $ nd $ experimental data, we use proton-deuteron ($ pd $) data in Fig.7 at $ E_N=35 $ MeV. The NLO corrections reverse the wrong trend and move toward the data. This suggests that for the power counting we adopt, $ A_y $ can only be described well at quite high orders.

In Ref. [53], the EKM chiral potentials [61] are used to calculate $ Nd $ scattering, and the subleading potentials are treated nonperturbatively. Note that "NLO" in the EKM scheme refers to $ \mathcal{O}(Q^2) $ corrections to LO and therefore includes more physics than our NLO does, e.g., two-pion-exchange potentials. By contrast, the perturbative calculations in this work allow us to access a wider cutoff window to test RG invariance, as in the $ NN $ sector. The LO+NLO potentials used here are somewhat comparable to the EKM LO; both have OPE as the long-range force. In this work, however, OPE is nonperturbative only in $ {{}^{{{1}}} {{S}}_{{{0}}}} $, $ {{{}^{{{3}}} {{S}}_{{{1}}}}-{{}^{{{3}}} {{D}}_{{{1}}}}} $, and $ {{}^{{{3}}} {{P}}_{{{0}}}} $, and the setup of the contact interactions differs: our scheme includes a $ {{}^{{{3}}} {{P}}_{{{0}}}} $ term at LO and a momentum-dependent $ {{}^{{{1}}} {{S}}_{{{0}}}} $ term at NLO. While the NLO of this work and the EKM LO yield similar angular distributions, they differ in their description of the analyzing power: our NLO exhibits gradual improvements toward the experimental data. It will be worthwhile in future work to identify which factor drives this difference: the perturbative treatment of OPE in higher partial waves or the arrangement of contact interactions.

We have performed calculations of $ Nd $ elastic scattering using renormalization-group-invariant chiral nuclear forces up to next-to-leading order, as developed in Refs. [17, 18, 20]. In our implementation of the Faddeev equation, a contour-deformation technique is employed, avoiding subtractions of the singularities of the integral equation. Strict perturbative calculations for the NLO potential after treating the LO nonperturbatively are performed. We demonstrated that computational costs can be substantially reduced in perturbation theory by leveraging the fact that the LO potentials are restricted to just a few $ NN $ partial waves. This approach yields a hierarchy of linear equations characterized by a fixed kernel — determined solely by the LO $ NN $ channels and their matrix elements — with driving terms that incorporate the higher-order potentials. In general, the computation time and memory usage of FKPT are at least an order of magnitude smaller than those of the auxiliary-potential method. The WPCD method was first used to benchmark our techniques. Good agreement was found for the phase shifts, elastic scattering differential cross sections, and nucleon analyzing powers.

We used the FKPT method to calculate the $ Nd $ scattering phase shifts for the doublet and quartet S waves. The phase shifts were found to converge with respect to the cutoff Λ up to 1600 MeV. This finding is consistent with earlier work, where RG invariance of the triton binding energy and $ nd $ scattering lengths indicated that three-nucleon forces are not required for renormalization purposes up to NLO. In comparisons with experimental data, we found that although the NLO calculation underpredicts the differential cross section at forward angles relative to LO, it generally yields better agreement for the analyzing power $ A_y $.

In anticipation of studying $ Nd $ scattering at next-to-next-to leading order ({ N$ ^2 $LO}), we note that, with a similar power counting, the triton binding energy exhibits severe cutoff dependence at this order near the so-called exceptional cutoffs [26]. These cutoffs were first reported in Ref. [62] in the context of $ NN $ scattering and the issue was discussed further in Refs. [63−67]. Because $ Nd $ scattering will allow us to access more channels of the $ 3N $ system than the triton, we expect the techniques presented here to help us exam whether such cutoff variations persist in other $ 3N $ channels with different quantum numbers.

We would like to thank Zeyuan Ye, Kaifei Ning, and Jinlong Dang for their contributions in the early stages of developing the Faddeev code used in this paper.