Efficient emulator for solving three-nucleon continuum Faddeev equations with chiral three-nucleon force comprising any number of contact terms

We demonstrate a computational scheme which drastically decreases the required time to get theoretical predictions based on chiral two- and three-nucleon forces for observables in three-nucleon continuum. For a three-nucleon force containing N short-range terms all workload is reduced to solving N+1 Faddeev-type integral equations. That done, computation of observables for any combination of strengths of the contact terms is done in a flash. We demonstrate on example of the elastic nucleon-deuteron scattering observables the high precision of the proposed emulator and its capability to reproduce exact results.

Gaussian regulator has been applied to the minimal set of independent contact interactions.
This new approach can be straightforwardly utilised to regularize also three-nucleon (3N) forces. That new family of semilocal chiral potentials provides an outstanding description of the NN data.
Applications of the EFT approach in the form of chiral perturbation theory (ChPT) have resulted not only in the theoretically well grounded NN potentials but also for the first time have given a possibility to apply in practical calculations NN forces augmented by consistent 3N interactions, derived within the same formalism. Understanding of nuclear spectra and reactions based on these consistent chiral two-and many-body forces has become a hot topic of present day few-body studies [12].
The first nonvanishing contributions to the 3N force (3NF) appear at next-to-next-toleading order of chiral expansion (N 2 LO) [3,13] and comprise in addition to the 2π-exchange term two contact contributions with strength parameters c D and c E [14]. The difficult task to derive the chiral 3NF at N 3 LO was accomplished in [15,16]. At that order five different topologies contribute to 3NF. Three of them are of long-range character [15] and are given by two-pion (2π) exchange graphs, by two-pion-one-pion (2π − 1π) exchange graphs, and by the ring diagrams. They are supplemented by the short-range two-pion-exchange-contact (2π-contact) term and by the leading relativistic corrections to 3NF [16]. The 3NF at N 3 LO order does not involve any new unknown low-energy constants (LECs) and depends only on two parameters, c D and c E that parameterize the leading one-pion-contact term and the 3N contact term present already at N 2 LO. The c D and c E values need to be then fixed at this order, as at N 2 LO, from a fit to few-nucleon data. At the higher order, N 4 LO, in addition to long-and intermediate-range interactions generated by pion-exchange diagrams [17,18], the chiral N 4 LO 3NF involves thirteen purely short-range operators, which have been worked out in [19].
Since the advent of numerically exact three-nucleon continuum Faddeev calculations the elastic nucleon-deuteron (Nd) scattering and the deuteron breakup reaction have been a powerful tool to test modern models of the nuclear forces [20][21][22] and the question about the importance of 3NF has developed into the main topic of 3N system studies. That issue has been given a new impetus by the ChPT-based achievements and the possibility to apply consistent two-and many-body nuclear forces, derived within this framework, in 3N continuum calculations.
3 Using chiral 3NF in 3N continuum requires numerous time consuming computations with varying strengths of the contact terms in order to establish their values. They can be determined for example from the 3 H binding energy and the minimum of the elastic Nd scattering differential cross section at the energy (E lab ≈ 70 MeV), where the effects of 3NF start to emerge in elastic Nd scattering [23,24]. Specifically at N 2 LO, after establishing the so-called (c D , c E ) correlation line, which for a particular chiral NN potential combined with a N 2 LO 3NF gives pairs of (c D , c E ) values reproducing the 3 H binding energy, a fit to experimental data for the elastic Nd cross section is performed to determine the c D and c E strengths. Fine-tuning of the 3N Hamiltonian parameters requires an extensive analysis of available 3N elastic Nd scattering and breakup data. That ambitious goal calls for a significant reduction of computer time necessary to solve the 3N Faddeev equations and to calculate the observables. Thus finding an efficient emulator for exact solutions of the 3N Faddeev equations seems to be essential and of high priority.
In Ref. [25] we proposed such an emulator which enables us to reduce significantly the required time of calculations. We tested its efficiency as well as ability to accurately reproduce exact solutions of 3N Faddeev equations. In the present study we introduce a new computational scheme, based on the perturbative approach of [25], which even by far more reduces the computer time required to obtain the observables in the elastic nucleon-deuteron scattering and deuteron breakup reactions at any energy, and which is well-suited for calculations with varying strengths of the contact terms in a chiral 3NF. Before presenting this new emulator, for the reader's convenience we shortly outline the main points of the 3N Faddeev formalism and of the perturbative treatment of Ref. [25]. For details of the formalism and numerical performance we refer to [20,[26][27][28].
Neutron-deuteron (nd) scattering with nucleons interacting via NN interactions v N N and a 3NF V 123 = V (1) + V (2) + V (3) , is described in terms of a breakup operator T satisfying the Faddeev-type integral equation [20,26,27] T |φ = tP |φ + (1 + tG 0 )V (1) (1 + P )|φ + tP G 0 T |φ (1) The 2N t-matrix t is the solution of the Lippmann-Schwinger equation with the interaction v N N . V (1) is the part of a 3NF which is symmetric under the interchange of nucleons 2 and 3: V 123 = V (1) (1 + P ). 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 initial state |φ = | q 0 |φ d describes the free motion of the neutron and the deuteron with the relative momentum q 0 and contains the internal deuteron wave function |φ d . G 0 is the free threebody resolvent. The amplitude for elastic scattering leading to the final nd state |φ ′ is then given by [20,27] while the amplitude for the breakup reaction reads where the free breakup channel state | p q is defined in terms of the Jacobi (relative) momenta p and q.
We solve Eq. (1) in the momentum-space partial-wave basis |pqα , determined by the magnitudes of the Jacobi momenta p and q and a set of discrete quantum numbers α comprising the 2N subsystem spin, orbital and total angular momenta s, l and j, as well as the spectator nucleon orbital and total angular momenta with respect to the center of mass (c.m.) of the 2N subsystem, λ and I: The total 2N and spectator angular momenta j and I as well as isospins t and 1 2 , are finally coupled to the total angular momentum J and isospin T of the 3N system. In practice a converged solution of Eq. (1) using partial wave decomposition in momentum space at a given energy E requires taking all 3N partial wave states up to the 2N angular momentum j max = 5 and the 3N angular momentum J max = 25 2 , with the 3N force acting up to the 3N total angular momentum J = 7/2. The number of resulting partial waves (equal to the number of coupled integral equations in two continuous variables p and q) amounts to 142.
The required computer time to get one solution on a personal computer is about ≈ 2 h.
In the case when such calculations have to be performed for a big number of varying 3NF parameters, time restrictions become prohibitive. Fortunately, the perturbative approach of Ref. [25] leads to a significant reduction of the required computational time.
Let us consider a chiral 3NF at a given order of chiral expansion with variable strengths of its contact terms. The 3NF at N 2 LO has one parameter-free term (2π-exchange contribution) 5 and two short-range terms with strength parameters c D and c E . At N 3 LO there are more contributing parameter-free parts but again only two contact terms. At N 4 LO parameterfree contributions are supplemented by fifteen short-range terms with strengths: c D , c E , c E 1 , ..., c E 13 . All these contact terms are restricted to small 3N total angular momenta and to only few partial wave states for a given total 3N angular momentum J and parity π. For example for J π = 7/2 ± all matrix elements < pqα|V (1) |p ′ q ′ α ′ > proportional to c E 1 and c E 7 vanish, while the c D and c E terms are nonzero only for a restricted number of α, α ′ pairs (mostly these containing 1 S 0 and 3 S 1 − 3 D 1 quantum numbers) [13,14]. Bearing that in mind and taking into account the fact that contact terms yield a small contribution to the 3N potential energy compared to the leading, parameter-free part, it is possible to apply a perturbative approach in order to include the contact terms.
In the next step, the proper perturbative treatment is performed, solving first the second equation in set (7). Having determined α|∆T (θ)|φ from Eq.(8) the emulator solution of Eq. (9) is calculated (in the following this emulator will be denoted by E∆T ). That allows one to reduce the required computation time and to reproduce surprisingly well the exact predictions for neutron-deuteron (nd) elastic scattering as well as for nd breakup observables [25]. To be specific, taking set |β which includes all 2N states with the total 2N angular momenta j ≤ 2, leads to a reduction of the computing time by a factor of approximately 4 in comparison to the exact calculations. Note that it takes approximately 30 minutes on a personal computer to solve Eq. (1), provided that the V (θ 0 )(1 + P ) and V (θ i )(1 + P ) kernels, acting in (1 + tG 0 )V (θ)(1 + P )G 0 T (θ)|φ term of Eq. (1), are prepared in advance, with the strengths θ i = (c i = 1, c k =i = 0).
In spite of such a large reduction, the computational time can be even further decreased and calculation of 3N continuum observables made in a flash. This notion is based on the observation that among three kernel-terms in the second equation of set (7), it is possible (because of the smallness of ∆V (θ)) to neglect the term β|(1 + tG 0 )∆V (θ)(1 + P )G 0 ∆T (θ)|φ .
The resulting integral equation for β|∆T (θ)|φ : β|∆T (θ)|φ = β|(1 + tG 0 )∆V (θ)(1 + P )|φ permits one to transfer the linear dependence on the strengths c i from the ∆V (θ) on the ∆T (θ). Namely, let β|∆T i |φ be a solution of Eq.(10) for a set θ i = (c i = 1, c k =i = 0): Multiplying (11) by c i and summing over i one gets: and the solution of Eq. (10) is given by: In this way at a given energy the computation of observables in the elastic Nd scattering and deuteron breakup reaction for any combination of strengths c i of contact terms is reduced to solving once N + 1 Faddeev equations: one equation for T (θ 0 ) and N equations for ∆T i .
In the first step, solution for α(β)|T (θ 0 )|φ is found. Then Eq. (11) is solved for β|∆T i |φ , from which the α|T i |φ is calculated by: The above computations need to be done only once and then for any combination of the strengths c i α(β)|T (θ = (c i , i = 1, . . . , N))|φ is obtained by trivial summation: For a calculation of elastic scattering observables the required sum of the second and the third term in Eq. (2) is obtained by: The above computational scheme forms the new emulator, which will be denoted in the following by E∆T i . To check its quality and efficiency as well as to compare it with E∆T we have chosen, as in Ref. [25], the SMS N 4 LO + chiral potential of the Bochum group [11], with the regularization cutoff Λ = 450 MeV, and combined it with the chiral N 2 LO 3NF.
We solved the 3N Faddeev equation (1) exactly at two incoming neutron energies E = 70 and 190 MeV with that combination as well as with this NN potential supplemented with the parameter free 2π-exchange term of the N 2 LO 3NF (set θ 0 = (c D = 0, c E = 0)). The first energy was taken from a region where 3NF effects start to appear in 3N continuum observables and the second one from a range with well-developed 3NF effects [23,24]. The (according to the notation of Refs. [13,14]). The precision of both emulators at both energies is similar and amounts to ≈ 1 − 2 %, with E∆T being slightly more precise than To demonstrate the power of the new emulator we show in Fig. 2  well as the short-range D-term are important at both energies.  dash-double-dotted curves to changes by a 3NF parameter free + D and parameter free + E terms, respectively.