Low-energy neutron-deuteron reactions with N 3 LO chiral forces

We solve three-nucleon Faddeev equations with nucleon-nucleon and three-nucleon forces derived consistently in the framework of chiral perturbation theory at next-to-next-to-next-to-leading order in the chiral expansion. In this first investigation we include only matrix elements of the three-nucleon force for partial waves with the total two-nucleon (three-nucleon) angular momenta up to 3 (5/2). Low-energy neutron-deuteron elastic scattering and deuteron breakup reaction are studied. Emphasis is put on Ay puzzle in elastic scattering and cross sections in symmetric-space-star and neutron-neutron quasi-free-scattering breakup configurations, for which large discrepancies between data and theory have been reported.


Introduction
A special place among few-body systems is reserved for the three-nucleon (3N) system, for which a mathematically sound theoretical formulation in the form of Faddeev equations exists, both for bound and scattering states. Over the past few decades algorithms have been developed to solve numerically three-nucleon Faddeev equations for any dynamical input which, in addition to nucleon-nucleon (NN) interactions, also involves threenucleon forces (3NFs) [1][2][3]. Using these algorithms and standard, (semi)phenomenological nucleon-nucleon interactions alone or supplemented by three-nucleon force model, numerous investigations of 3N bound states and reactions in the 3N continuum have been carried out. Highprecision nucleon-nucleon potentials such as the AV18 [4], a e-mail: witala@if.uj.edu.pl CD Bonn [5], Nijm I and II [6] NN forces, which provide a very good description of the nucleon-nucleon data up to about 350 MeV, have been used. They have also been combined with model 3N forces such as the 2πexchange Tucson-Melbourne (TM99) 3NF [7] or the Urbana IX model [8].
When realistic NN forces are used to predict binding energies of three-nucleon systems they typically underestimate the experimental bindings of 3 H and 3 He by about 0.5-1 MeV [10,9]. This missing binding energy can be corrected for by introducing a three-nucleon force into the nuclear Hamiltonian [9]. Also the study of elastic nucleondeuteron (Nd) scattering and nucleon induced deuteron breakup revealed a number of cases where the nonrelativistic description using only pairwise forces is insufficient to explain the data. The best studied case at low energies is the vector analyzing power in elastic Nd scattering for which a large discrepancy exists in the region of its max-imum around c.m. angles θ c.m. ∼ 125 • and for incoming nucleon energies below ∼ 20 MeV [2,11,12]. For the elastic scattering angular distribution at such energies, negligible effects of 3NF's have been found and theory based on realistic NN forces agrees well with the data [2,11].
That picture changes with increasing energy of the Nd system. Generally, the studied discrepancies between experiment and theory using only NN potentials become larger and adding a three-nucleon force to the pairwise interactions leads in some cases to a better description of the data. The elastic Nd angular distribution in the region of its minimum and at backward angles is the best known example [13,14]. The clear discrepancy in these angular regions at energies up to E lab, N ∼ 100 MeV between a theory using only NN potentials and the cross section data can be removed by adding standard models of threenucleon forces to the nuclear Hamiltonian. Such 3NFs are adjusted for a given NN potential to reproduce the experimentally observed binding energy of 3 H and 3 He [13,11,14]. At energies higher than ∼ 100 MeV current 3NFs only partially improve the description of cross section data and the remaining discrepancies, which increase with energy, indicate the possibility of relativistic effects. The need for a relativistic description of three-nucleon scattering was realized when precise measurements of the total cross section for neutron-deuteron (nd) scattering [15] were analyzed within the framework of nonrelativistic Faddeev calculations [16]. NN forces alone were insufficient to describe the data above ∼ 100 MeV. The effects due to the relativistic kinematics considered in [16] at higher energies were comparable in magnitude to the effects due to 3NFs. These results provided further motivation to study relativistic effects in the three nucleon continuum in a systematic way.
Subsequent studies of relativistic effects in the threenucleon continuum [17][18][19][20] revealed, that when the nonrelativistic form of the kinetic energy is replaced by the relativistic form and a proper treatment of the relativistic dynamics is introduced, the elastic scattering cross section is only slightly increased at backward angles and higher energies while spin observables are practically unchanged. These results led to the conclusion that discrepancies between data and theory at higher energies must reflect the action of 3NF's which have to be included in the nuclear Hamiltonian.
The main drawback of all those studies was inconsistency between applied NN interactions and 3NFs. With the advent of effective field theoretical methods in the form of chiral perturbation theory, it became possible to construct consistent two-and many-nucleon forces. In this way an exciting possibility to study few-nucleon systems and their reactions with consistent two-and many-nucleon interactions has emerged.
In [21], the above-mentioned inconsistency was removed and low-energy 3N continuum were investigated with chiral next-to-next-to-leading order (N 2 LO) NN and 3N forces. The NN interaction in that order, however, does not describe the NN experimental phase-shifts in a sufficiently wide energy range to allow application of those forces at higher energies.
In [22][23][24] and [25,26], precise two-nucleon potentials have been developed at next-to-next-to-next-to-leading order (N 3 LO) of the chiral expansion. They reproduce experimental NN phase-shifts [27,28] in a wide energy range and practically with the same high precision as realistic (semi)phenomenological NN potentials. The necessary work to derive the consistent chiral 3NF's at N 3 LO has been accomplished in [29,30] and [31]. At that order, six different topologies contribute to the 3NF. Three of them are of a long-and intermediate-range [29] and are given by two-pion (2π) exchange, two-pion-one-pion (2π-1π) exchange and the so-called ring diagrams. They are supplemented by the shorter-range 1π-contact and threenucleon-contact components, which appear first at N 2 LO, by the two-pion-exchange-contact (2π-contact) term as well as by the leading relativistic corrections to the threenucleon force [30].
The results of refs. [25,22,23,29,30] enable one to perform, for the first time, consistent calculations of threenucleon reactions at N 3 LO order of chiral expansion. The 3NF at this order does not involve any new unknown lowenergy constants (LECs) and depends only on two free parameters, c D and c E , which parametrize the strengths of the leading 1π-contact and the three-nucleon-contact terms. Their values need to be fixed (at given order) from a fit to a few-nucleon data. Among the few possible observables that have been used in this connection one can mention the triton binding energy, the nd doublet scattering length 2 a nd [21], the 4 He binding energy [32,33] along with the point proton rms radius [34], the properties of light nuclei, or the triton β decay rate [35]. Notice that the first three observables are known to be strongly correlated and therefore might not be the best choice for the determination of c D and c E . Application of N 3 LO 3NF in few-body calculations is challenging due to its very rich and complicated operator structure. The large number of terms in the 3NF at N 3 LO [29,30] requires an efficient method of performing partial-wave decomposition. Recently such a method, which runs under the name of automatized partial-wave decomposition (aPWD), was proposed in [36][37][38]. In that approach, the matrix elements in the 3N momentum-space partial wave basis for different terms contributing to N 3 LO 3NF are obtained in two consecutive steps. First, the spinmomentum and isospin parts of three-nucleon interactions are computed using symbolic algebra systems. The resulting momentum-dependent functions are then integrated numerically in five dimensions over angular variables. The major advantage of this method is its generality since it can be applied to any momentum-spin-isospin operator. Application of that method for higher angular momenta requires large computer resources. Therefore, in this first study of the 3N continuum with full N 3 LO chiral force, we restrict ourselves to low energies only. In that region of incoming neutron laboratory (lab.) energies below ∼ 30 MeV, the most challenging observables are the nd elastic scattering analyzing power and cross sections in symmetric space star and neutron-neutron quasifree-scattering configurations of the nd breakup reaction. The discrepancies between data and theory for these ob-servables could not be removed with standard NN and 3NFs [39].
Our paper is organized as follows. In sect. 2 we describe our method to determine the nuclear Hamiltonian by fixing the two parameters c D and c E in the chiral N 3 LO 3NF. This is achieved by first requiring reproduction of the 3 H binding energy which leads to pairs of allowed (c D , c E ) values. Using the experimental data for an additional 3N observable, which in our case is taken to be the doublet nd scattering length 2 a nd , fixes completely the nuclear Hamiltonian at N 3 LO. Based on the resulting Hamiltonian, we discuss in sect. 3 some results for lowenergy elastic nd scattering observables while in sect. 4 the results for selected low-energy nd breakup configurations are presented. We summarize and conclude in sect. 5.

Determination of nuclear Hamiltonian at N 3 LO
Neutron-deuteron scattering with neutrons and proton interacting through a NN interaction v NN and a 3NF V 123 = V (1) (1 + P ), is described in terms of a breakup operator T satisfying the Faddeev-type integral equation [1][2][3] T |φ = tP |φ The two-nucleon t-matrix t is the solution of the Lippmann-Schwinger equation with the interaction v NN . V (1) is the part of a 3NF which is symmetric under the interchange of nucleons 2 and 3. The permutation operator P = P 12 P 23 + P 13 P 23 is given in terms of the transposition operators, P ij , which interchange nucleons i and j. The incoming state |φ = |q 0 |φ d describes the free nd motion with relative momentum q 0 and the deuteron state |φ d . Finally, G 0 is the resolvent of the three-body center-of-mass kinetic energy. The amplitude for elastic scattering leading to the corresponding two-body final state |φ is then given by [2,3] while for the breakup reaction one has where |φ 0 is the free three-body breakup channel state. The nuclear Hamiltonian at N 3 LO of the chiral expansion is fixed by specifying the values of LECs c D and c E which parametrize the strengths of the leading 1π-contact and the three-nucleon-contact terms. To determine them we follow the approach of ref. [21] and use the experimental triton binding energy E( 3 H) and the nd doublet scattering length 2 a nd as two observables from which c D and c E can be obtained. The procedure can be divided into two steps. First, the dependence of E( 3 H) on c E for a given value of c D is determined. The requirement to reproduce the experimental value of the triton binding energy yields a set of pairs c D and c E . This set is then used in the calculations of 2 a nd , which allows us to find the values of c D and c E describing both observables simultaneously. As already emphasized above, using the triton binding energy and the nd doublet scattering length is probably not the optimal way to fix the parameters in the 3NF due the strong correlation between these two observables. We will discuss this issue in the next two sections and present results obtained by relaxing the condition to reproduce 2 a nd .
We compute the 3 H wave function using the method described in [9], where the full triton wave function |Ψ = (1 + P )|ψ is given in terms of its Faddeev component ψ, which fulfills the Faddeev equation The doublet scattering length 2 a nd is calculated using (c D ,c E ) pairs, which reproduce the correct value of E( 3 H).
To this end, we solve the Faddeev equation (1) for the auxiliary state T |φ at zero incoming energy [40]. We refer to [2,3,41] for a general overview of 3N scattering and for more details on the practical implementation of the Faddeev equations.
In this first study, where the full N 3 LO 3NF is applied, we restrict ourselves to nd reactions at low energies, E lab, n < 20 MeV. At such low energies it is sufficient to include NN force components with a total two-nucleon angular momenta j ≤ 3 in 3N partial-wave states with the total 3N system angular momentum below J ≤ 25/2. For the 3NF it is sufficient to incorporate its matrix elements with j ≤ 3 and J ≤ 5/2.
Here and in what follows, we employ the N 3 LO chiral NN potential of refs. [22,23]. From among five versions corresponding to different sets of cut-off parameters used to regularize the Lippmann-Schwinger equation and in spectral function regularization, namely (450, 500) MeV, (450, 700) MeV, (550, 600) MeV, (600, 500) MeV, and (600, 700) MeV, we applied for the present study two N 3 LO chiral NN potentials with cut-off sets (450, 500) MeV and (450, 700) MeV, denoted in the following by 201 and 204, respectively. Only for these two sets of cut-offs were we able to determine the LECs c D and c E using our procedure. to the two cut-off sets used in the present study. Namely in [43] an application of N 3 LO 3NF, however with relativistic 1/m corrections and short-range 2π-contact term omitted, also led to unnaturally large c D values for all five cut-off combinations. We hope that new generations of chiral forces with other regularization schemes will cure this problem [44]. We also plan to use other 3N observables, for example triton β-decay rate instead of 2 a nd , to fix values of LECs c D and c E .

Low-energy elastic nd scattering
At low energies of the incoming neutron, the most interesting observable is the analyzing power A y for nd elastic scattering with polarized neutrons. Theoretical predictions of standard high-precision NN potentials fail to ex-plain the experimental data for A y . The data are underestimated by ∼ 30% in the region of the A y maximum which occurs at c.m. angles Θ c.m. ∼ 125 • . Combining standard NN potentials with commonly used models of a 3NF, such as, e.g. the TM99 or Urbana IX models, removes approximately only half of the discrepancy with respect to the data (see fig. 2). When instead of standard forces chiral NN interactions are used, the predictions for A y vary with the order of chiral expansion [22,23]. In particular, as reported in ref. [21], the NLO results overestimate the A y data while N 2 LO NN forces seem to provide quite a good description of them (see fig. 2). Only when N 3 LO NN chiral forces are used, a clear discrepancy between theory and data emerge in the region of A y maximum, which is similar to the one for standard forces. This is visualized in fig. 2, where bands of predictions for five versions of the Bochum NLO, N 2 LO and N 3 LO potentials with different cut-off parameters used for the Lippmann-Schwinger equation and the spectral function regularizations are shown [23]. Such behaviour of A y predictions at different orders in the chiral expansion can be traced back to a high sensitivity of A y to 3 P j NN force components and to the fact, that  only at N 3 LO of chiral expansion the experimental 3 P j phases [27,28], especially the 3 P 2 -3 F 2 ones, are properly reproduced [43]. It is interesting to study whether the consistent chiral N 3 LO 3NF's can explain the low-energy A y -puzzle. In the present investigation, we, for the first time include all contributions to N 3 LO 3NF: long-range contributions comprising 2π exchange, 2π-1π exchange, ring components and relativistic 1/m corrections together with short range 1π-contact, three-nucleon-contact and 2π-contact terms. In fig. 3 we show by dash-dotted (blue) line the results for A y based on the values of the c E and c D parameters which reproduce the triton binding energy and 2 a nd scattering length. It turns out that adding the full N 3 LO 3NF does not improve the description of A y . On the contrary, adding the chiral N 3 LO 3NF lowers the maximum of A y with respect to the chiral N 3 LO NN prediction, shown by the solid (red) line, thus, increasing the discrepancy between the theory and the data.
In order to check the restrictiveness of the requirement to reproduce, in addition to the 3 H binding energy, also the experimental value of 2 a nd , we show in fig. 3 also a band of predictions for (c E , c D ) pairs from fig. 1(a) and (b). Even after relaxing the requirement to reproduce 2 a nd , the A ypuzzle cannot be explained by the N 3 LO NN and 3NF.
It is interesting to see how different components of the N 3 LO 3NF contribute to A y . Taking in addition to the NN N 3 LO chiral force only the 2π-exchange term with leading 1π-contact and three-nucleon-contact terms (these three topologies appear for the first time at N 2 LO) lowers the maximum of A y (see fig. 4, solid (cyan) line). When, in addition, the short-range 2π-contact component is included, the value of A y practically remains unchanged (dash-dotted (magenta) line in fig. 4). This shows that contributions of the 2π-contact term are negligible at those energies. The long-range 2π-1π exchange and ring terms lower significantly the maximum of A y (in fig. 4 dotted (maroon) and dashed (green) lines, respectively).  Finally, inclusion of the relativistic 1/m contribution leaves the maximum of A y practically unchanged (dashdouble-dotted (blue) line in fig. 4). It should be pointed out that when taking into account the 1/m corrections to the N 3 LO 3NF, one should also include the corresponding relativistic corrections in the NN force and, in addition, also relativistic corrections to the kinetic energy, which are formally of the same importance. This would considerably complicate the calculation. In our present work, we do not take into account such corrections and employ the standard nonrelativistic framework. This seems to be justified in view of the low energies considered in this paper and the very small effects caused by relativistic 1/m corrections to the 3NFs found in this study. Last but not least, we emphasize that the contributions of the individual 3NF topologies to the A y puzzle are not observable and depend, in particular, on the regularization scheme and employed NN forces.
It is important to address the question of uniqueness of our approach to determine the constants c D and c E . To this aim, we checked how taking instead of 2 a nd a different nd observable would influence determination of c D and c E . The low-energy elastic nd scattering cross section is an observable which seems to be reasonably well described by standard theory [47]. In fig. 5 we show (orange) bands of predictions for the nd elastic scattering cross section at E lab, n = 6.  tions for the full N 3 LO chiral force with constants c D and c E fixed by requirement that the doublet nd 2 a nd scattering length is also reproduced. For comparison to standard potential cross sections in fig. 5 also the CD Bonn potential results are shown by solid (blue) lines. The backward angle nd elastic scattering cross section data are properly described by standard, high-precision NN potentials [47].
To fix values of c D and c E it would be desirable to have forward angle cross section data. Assuming that in this angular region the data will be properly described by our theory indicates that replacing 2 a nd by cross section would lead to consistent c D and c E values in both approaches.

Low-energy nd breakup
Among numerous kinematically complete configurations of the nd breakup reaction the SST and QFS configurations have attracted special attention. The cross sections for these geometries are very stable with respect to the underlying dynamics. Different potentials, alone or combined with standard 3NFs, lead to very similar results for the cross sections [39] which deviate significantly from available SST and neutron-neutron (nn) QFS data. At low energies, the cross sections in the SST and QFS configurations are dominated by the S-waves. For the SST configuration, the largest contribution to the cross section comes from the 3 S 1 partial wave, while for the nn QFS the 1 S 0 partial wave dominates. Neglecting rescattering, the QFS configuration resembles free NN scattering. For free, low-energy neutron-proton (np) scattering one expects contributions from 1 S 0 np and 3 S 1 force components. For free nn scattering, only the 1 S 0 nn channel is allowed. This suggests that the nn QFS is a powerful tool to study the nn interaction. The measurement of np QFS cross sections have revealed good agreement between the data and theory [48], thus confirming the knowledge of the np force. For the nn QFS it was found that the  theory underestimates the data by ∼ 20% [48]. The large stability of the QFS cross sections with respect to the underlying dynamics means that, assuming correctness of the nn QFS data, the present day 1 S 0 nn interaction is probably incorrect [39,49,50].
Also the chiral N 3 LO forces with all components of the 3NF included are not an exception and cannot explain the discrepancy between the theory and data found for the SST configuration [51] (fig. 6). The solid (black) line shows the cross section when only NN chiral N 3 LO force is active. Adding the full N 3 LO 3NF with c D and c E pairs reproducing the experimental binding energy of 3 H and nd doublet scattering length leads to dash-double-dotted (blue) line. At 13 MeV, it lies only slightly below the NN potential prediction indicating only small 3NF effects at this energy.
It is interesting to see how the SST cross section depends on the choice of parameters (c D , c E ) which enter the N 3 LO nuclear Hamiltonian. In fig. 6, the SST cross sec-tions at E lab, n = 13 MeV are shown for a number of c D and c E pairs which reproduce only the experimental binding energy of 3 H (taken from fig. 1(a) and (b)). For the 201 N 3 LO nuclear Hamiltonian (see fig. 6(a)) decreasing the value of c D leads to big changes of the SST cross section. Starting from c D = 13.78, which reproduce also 2 a nd , and decreasing it to c D = 9 leads to only small changes of the SST cross sections. Further lowering of c D down to c D = −3 reduces the cross section and the discrepancy to nd data at 13 MeV is drastically increased. If we continue to reduce the c D value the SST cross section rises, however, it remains always below the pure NN prediction. For the 204 N 3 LO nuclear Hamiltonian the changes of the SST cross section are not so drastic and decrease of the c D reduces the cross section (see fig. 6(b)). Thus, in spite of the strong sensitivity of the SST cross sections to values of c D and c E , it is not possible to describe the available experimental data for the SST nd cross sections at 13 MeV even allowing for pairs of (c D ,c E ) which do not reproduce 2 a nd . and [52,53], respectively. The (black) squares are proton-deuteron (pd) data of ref. [54].  fig. 7(a)). For the 204 N 3 LO nuclear Hamiltonian decreasing c D leads to the increase of the QFS cross section (see fig. 7(b)). The values of c D and c E which reproduce the 3 H binding energy and 2 a nd lead only to a slight increase of the QFS cross section with respect to the N 3 LO NN prediction and thus to small 3NF effects.
Recent efforts towards the derivation and implementation of the N 3 LO 3NF allowed us, for the first time, to apply the full chiral N 3 LO Hamiltonian to the lowenergy nd elastic scattering and breakup reactions. The nuclear Hamiltonian at that order of the chiral expansion is unambiguously given after fixing the two constants c D and c E which determine the strengths of the 1π-contact and three-nucleon-contact components of the N 3 LO chiral 3NF. We determined these low-energy constants by requiring reproduction of the binding energy of 3 H and the doublet nd scattering length 2 a nd . We found indications that using low-energy nd elastic scattering cross section instead of 2 a nd would probably lead to similar values of these parameters.
It turns out that applying the full N 3 LO 3NF with specific cut-off parameters used in this study cannot explain the low-energy A y -puzzle. Contrary to the 3NF effects found for A y with standard NN potentials combined with 3NF models such as TM99 or Urbana IX, where the inclusion of the 3NF decreased the discrepancy to data by about ∼ 50%, the chiral N 3 LO 3NF combined with the NN potential of ref. [22] lowers the maximum of A y increasing the discrepancy. It should, however, be emphasized that the low-energy 3N A y is a fine-tuned observable which is very sensitive to changes in 3 P j NN force components as well as to P -waves in the Nd system [55,56]. Thus, the disagreement with the data must be interpreted with considerable caution. Our result suggests the lack of some spin-isospin-momenta structures in the N 3 LO 3NF. However, possible inaccuracies in lowenergy 3 P j NN phase-shifts cannot be excluded. The 3NF derived in the standard formulation of chiral perturbation theory based on pions and nucleons as the only explicit degrees of freedom is known to miss certain significant intermediate-range contributions of the Δ(1232) resonance at N 3 LO, which, to some extent, are accounted for only at N 4 LO and higher orders [57,58]. It would therefore be interesting, to apply the recently derived N 4 LO 3NF [57,58] in calculations of nd reactions together with subleading contributions to the three-nucleon contact interactions [59]. The short-range 3N forces at N 4 LO which contribute to Nd P-waves may solve the A y -puzzle in a trivial way.
We found that cross sections in kinematically complete SST and QFS nd breakup configurations at low energies are quite sensitive to the values of c D and c E . For their values fixed by the experimental binding energy of 3 H and 2 a nd only small 3NF effects were found in these configurations. Large discrepancies with the data remain in these configurations.
For the SST geometry at 13 MeV, there is a serious discrepancy between theory and two independent nd data sets of refs. [51,53] as well as between theory and protondeuteron (pd) data of ref. [54]. While the nd data lie ∼ 20% above the theory, the pd data lie ∼ 10% below theory and ∼ 30% below nd data. Recent pd calculations with Coulomb force included show practically negligible effects of the proton-proton Coulomb force for this config-uration [60]. The observed large splitting between the nd and pd data indicates either that there are large isospinbreaking effects or that the data are not consistent.
Higher-energy nd reactions, in which clear evidence of large 3NF effects was found, call for applications of the full N 3 LO force. Studies of the cut-off dependence of N 3 LO NN chiral interaction in higher-energy nd elastic scattering revealed preference for larger cut-off values [43]. The use of lower cut-offs would preclude applications of N 3 LO chiral dynamics in that interesting region of energies. It is important to address the issue of reducing finite-cut-off artifacts and increasing the accuracy of chiral nuclear forces prior to applying the chiral N 3 LO Hamiltonian at higher energies. In addition, one needs to explore different possibilities to determine the LECs entering the 3NF in view of the known strong correlations between, e.g. the 3 H and 4 He binding energies and the nd doublet scattering lengths, see [61] for a related discussion. Last but not least, more effort should be invested into providing a reliable estimation of the theoretical uncertainty at a given order in the chiral expansion. Work along these lines is in progress. Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.