Neutron-deuteron scattering cross-sections with chiral $NN$ interactions using wave-packet continuum discretization

In this work we present a framework that allows to solve the Faddeev equations for three-nucleon scattering using the wave-packet continuum-discretization method. We perform systematic benchmarks using results in the literature and study in detail the convergence of this method with respect to the number of wave packets. We compute several different elastic neutron-deuteron scattering cross-section observables for a variety of energies using chiral nucleon-nucleon interactions. For the interaction $\text{N2LO}_{\text{opt}}$ we find good agreement with data for nucleon scattering-energies $E_\text{Lab} \leq 70$ MeV and a slightly larger maximum of the neutron analyzing power $A_y(n)$ at $E_\text{Lab} = 10$ MeV and 21 MeV compared with other interactions. This work represents a first step towards a systematic inclusion of three-nucleon scattering observables in the construction of next-generation nuclear interactions.


I. INTRODUCTION
Nucleon-nucleon (N N ) and nucleon-deuteron (N d) scattering are prototypical processes for analyzing ab initio nuclear Hamiltonians [1,2]. While N N cross sections are straightforward to calculate and are nowadays routinely being used to calibrate modern N N interactions [3][4][5][6], the computation of three-nucleon (N N N ) scattering processes like N d scattering is much more demanding due to the presence of energy poles in the underlying equations, contributions from N N N interactions, and the existence of multiple reaction channels. In fact, it is a computationally challenging task to numerically solve the Faddeev equations [7] and its extensions [8][9][10][11][12] in a reliable and accurate way. That is why it was only in the late 1980's that realistic quantum scattering calculations [13] began to emerge. Due to this complexity, it has not yet been feasible to perform a simultaneous statistical analysis of N N and N N N interactions using N d and N N cross section data.
In this paper we present an implementation of the wave-packet continuum-discretization (WPCD) method [14] to solve the Faddeev equations with the chiral N2LO opt [15] interaction. The WPCD method is one of many bound-state techniques [16] for solving the multi-particle scattering problem. The main advantages of this method are: i) coarse-graining the continuum using a square-integrable basis smooths out all singularities and facilitates straightforward numerical solutions of the Faddeev equations for the scattering amplitude, ii) all on-shell energy dependence resides in a closed-form expression of the channel resolvent, and iii) once the wave-packet basis is antisymmetrized, which has to be done only once computationally, the computational cost of predicting scattering observables scales sublinearly with the number of scattering energies. This opens ways for efficient computation of coarse-grained N d predictions for several scattering energies. The accuracy of the solutions depends polynomially on the number of wave-packets used for discretizing the continuum.
To that end we also study the convergence of the WPCD results with respect to the number of employed wave-packet basis states.
We benchmark the WPCD results against published cross section results for the traditional Nijmegen-I N N interaction [17] and systematically compare and analyze different cross section observables at a variety of energies using the interactions Idaho-N3LO [18] and N2LO opt . Both interactions have a history of being routinely employed in ab initio studies of nuclear structure and nucleon-nucleus reactions. The latter one, N2LO opt , is a next-to-next-to-leading order chiral interaction optimized to reproduce N N scattering phase shifts, and yields an accurate description of low-energy N N scattering data up to 125 MeV scattering energy. More importantly, it also reproduces key nuclear properties such as the location of the oxygen neutron drip-line and calcium shell closures without having to invoke N N N interactions. N2LO opt also gives a rather good description of selected nuclear structure physics, transitions, and reaction data, see, e.g., Refs. [19][20][21]. Of course, a complete calculation requires N N N interactions and these correlations can play a pivotal role for obtaining realistic ab initio predictions of bound and continuum nuclear observables, see, e.g., Refs. [22][23][24][25][26][27][28][29]. From these observations it is therefore interesting to predict N d scattering observables with the N2LO opt interaction. In this work we pay particular attention to the low-energy neutron (n) analyzing power A y (n), where a long-standing puzzle [30] arXiv:2201.09600v1 [nucl-th] 24 Jan 2022 resides 1 , and the nd differential cross section dσ/dΩ at nucleon laboratory scattering energy E Lab = 64.5 MeV. The latter observable is known to depend sensitively on N N N interactions [31,32].
In Sec. II we present the formalism that we implemented to solve the Faddeev equations for elastic N d scattering and benchmark its convergence with respect to basis dimension. In Sec. III we present predictions for nd scattering cross sections using the N2LO opt interaction, and end with a summary and outlook in Sec. IV.

II. ELASTIC N d SCATTERING USING THE WPCD METHOD
In this section, we present i) the WPCD method for solving the N d Faddeev equations in momentum space (Sec. II A), ii) how to construct a WPCD-basis and its partial-wave expansion (Sec. II B), iii) our computational implementation for solving the resulting matrix equation (Sec. II C), and iv) a convergence analysis of the WPCD method (Sec. II D). All detailed expressions are relegated to the appendices A-E.

A. The Faddeev equations in momentum space
The Faddeev equations can be reduced to the Alt-Grassberger-Sandhas (AGS) equation [11], which for elastic N d scattering and without a N N N interaction can be written aŝ where E denotes the on-shell scattering energy and where we used the usual "odd-man-out" notation such that the index i here refers to the incoming nucleon relative to an antisymmetric state of a nucleon pair (jk), e.g., the deuteron, for unequal i, j, k ∈ {1, 2, 3}. Our goal is to calculate elastic cross-sections via the elastic transition operatorÛ i . The three operatorsÛ i are related via the permutation operatorsP ijk ≡P ijPjk : i is the full Hamiltonian,ĥ i =ĥ 0 +v i is the N N Hamiltonian,ĥ 0 is the kinetic energy of the pair, andĥ 0 i is the free Hamiltonian of the third nucleon relative to the pair. SinceÛ i for i = 1, 2, 3 in Eq. (1) are not independent it suffices to solve for only one of them, e.g., U 1 . For the most part we will also drop this subscript. This will hopefully avoid possible confusion with respect to subscripts denoting different basis states defined below.
In the WPCD method, the channel resolvent G(E) is diagonal and straightforward to evaluate analytically using a WPCD-basis of scattering states. This has the advantage of removing all complications from singularities that plague the Faddeev method formulated in a planewave basis. Such points are essentially averaged out when using wave packets to represent states in the continuum. Furthermore, the entire E-dependence of the scattering process resides in the channel resolventĜ(E) and multiple scattering energies can be accessed without inducing much computational overhead.
We end this section by linking the form of the AGS equation used in WPCD, Eq. (1), to its conventional formulation used as a starting point in standard Faddeev methods. Using thattĜ 0 ≡vĜ andv =Ĝ −1 0 enables us to replace the interactionv and the channel resolventĜ with a fully off-shell N N t−matrix and the free resolvent G 0 at the on-shell energy E. This latter replacement is necessary since the channel resolvent cannot be straightforwardly evaluated a priori using only a plane-wave basis. This also introduces an explicit energy-dependence in the N Nt-matrix and thereby in the integral kernel of the AGS equation. In addition to this, singularities arise in the representation of both these operators [2]. This can be dealt with using subtraction techniques. Such complications are avoided altogether in the WPCD method.

B. Setting up the WPCD basis
We define a free wave-packet (FWP) for a pair of particles with relative momenta p within some interval (bin) where |p is a plane-wave state with momentum p and normalization p |p = δ(p −p) p p . This normalization differs from the one used in Ref. [14]. Here, N i is a normalization constant. The function f (p) is a weighting function which allows us to define, for example, momentum wave-packets, f (p) = 1, or energy wave-packets, f (p) = p µ0 , where µ 0 is the reduced mass. In this work µ 0 (µ 1 ) denote the reduced mass for the two-body (three-body) system. The naming convention for the two kinds of wave packets indicate whether they correspond to eigenstates of the momentum operatorp, or the kinetic energy operatorĥ 0 , of the two-body system. In this work we use both kinds of wave packets since it is simpler to use and derive operator projections onto momentum FWPs, but the resolventĜ is evaluated in an energy wave packet basis, see App. B. In the three-body system we define momentum FWPs as where we use the bar-notation to denote wave packets with the momenta q ∈D j ≡ [q j , q j+1 ] of the third particle relative to the centre of mass (c.m.) of the pair.
In this work, we use the same number of wave packets, N WP , when discretizing the continuum of Jacobi momenta p and q. We discretize the continuum according to a Chebyshev grid, i.e., where we use α = 200 MeV and t = 1, which yields wave packets residing in momentum bins reaching momenta up to ∼ 10 GeV. The width of the momentum bins increases with i, such that the vast majority of the wave packets reside below momenta of ∼ 500 MeV, which is where we typically have the most relevant contributions from modern chiral N N -potentials. We work in a partial-wave representation of N N N states and introduce a spin-angular basis with total angular momentum J and isospin T , where L, S, J, and T denote the relative orbital angular momentum, spin, total angular momentum, and isospin, respectively, for the antisymmetric nucleon-pair system. The orbital angular momentum of the third (spin-1/2) nucleon relative to the c.m. of the pair systems is denoted with l and its total angular momentum is denoted with j. Each Jj-coupled channel has a total angular momentum J . We can therefore construct N N N partial-waves as All N N N partial-waves are equipped with a unique combination of good quantum numbers J and parity Π = (−1) L+l . In our calculations we explicitly break isospin T by including the charge dependence of the N N interaction in the 1 S 0 channel. The impact of this T = 3 2 − 1 2 isospin coupling on elastic N d scattering is very small [33]. On the other hand, the computational costs of including it is negligible.
The FWP states form a square-integrable basis, with appropriate long-range behavior to approximate scattering states [34]. It is also straightforward to represent matrix elements of the permutation operatorP and the N N potential operatorv 1 in a basis of such states. See App. B for details regarding this projection.
The elastic transition operator will be solved for in a basis of N N N scattering wave packets (SWP) defined as where |z i are scattering wave-packets (eigenstates) of the N N Hamiltonianĥ. The latter states can be approximated in a finite FWP-basis as where the (real) C ji coefficients are obtained via straightforward diagonalization of the N N Hamiltonian in a basis of FWPs |x i . The coefficients allow for straightforward transformation between FWP and SWP partialwave bases. From the diagonalization we obtain eigenvectors and eigenvalues, i.e., scattering wave packets |z i with eigenenergies i such thatĥ|z i = i |z i . The eigenenergies define the bin boundaries D i ≡ {E i , E i+1 } for the scattering wave-packets |z i [14]. We will refer to the (negative) energy bin corresponding to the deuteron bound state as |z i d and the corresponding N N N partialwaves with a deuteron channel as |α d . Although the FWP-basis is sub-optimal for describing bound states, it yields a very good approximation to scattering states which is more important here. For all computations in this work we use a spin-angular basis of positive and negative parity N N N partial waves with J ≤ 17/2 and J ≤ 3. This leads to 60 channels per N N N partial wave. As we will discuss in Sec. II D, we find that using N WP ≈ 125 wave-packets in both Jacobi momenta is more than sufficient for accurately computing low-energy elastic scattering observables with E Lab 100 MeV [2], which is the region we focus on in this work. There are exceptions and we will discuss those below.

C. Computational implementation
Naturally, we solve for the transition operatorÛ for each combination of N N N total angular momentum J and parity Π separately. We represent Eq. (1) in matrix form using a SWP-basis where A ≡ C T PVC. Here, we defined finite-dimensional matrices for the N N -potential matrix V and the permutation matrix P in a N N N momentum-FWP-basis. They are obtained using the expressions in App. B 2 and App. B 3, respectively. Note that V and C are block diagonal for momenta q in different binsD j and quantum numbers l and j. Once we have diagonalized the N N Hamiltonian, we construct an approximate SWP-basis and setup the (block-diagonal) matrix C of C ij coefficients in Eq. (9). The eigenvalues of G are easily obtained in the SWP-basis, see App. B 4. This is of key importance.
Formally, Eq. (10) is a matrix-equation that can be solved via inversion. However, straightforward inversion, or numerically stable equivalents, is unviable for realistic nuclear potentials since the matrix A is too large to be stored in memory for the basis sizes we require for convergence. Fortunately it is possible to store the matrices necessary to construct A in memory, i.e., C, V, and P. Indeed, P only has to be computed once and is very sparse (> 99%). We see that A is 100 times denser than P.
We solve Eq. (10) for the on-shell transition operator in the SWP-basis U i.e., the transition matrix elements corresponding to an incoming nucleon with on-shell momentum q ∈D j scattering elastically off a deuteron. We compute this amplitude by summing the first 20-30 terms of the Neumann (or Born) series where we have defined K(E) ≡ G(E)A. Note that G, and thereby also K, depend on the on-shell scattering energy E, and that since G is diagonal, A and K have identical densities. Thus, K must be re-computed in segments and on-the-fly for the repeated matrix-vector multiplications needed to generate the terms of the series above. We employ a Padé extrapolation [35] to handle a divergent Neumann series and we find that this rational approximant facilitates a convergent resummation in our case, see App. C. Next to the computational cost of initially constructing P, the cost of setting up the kernel K constitutes the numerical bottleneck in our current implementation of the WPCD method as it must be repeated several times. The product GA is trivial, which in turn makes it trivial to compute transition matrices at several different energies E with the WPCD method. The product C T PVC is a product of the sparse matrix P with the block-diagonal matrices C T and VC on either side. Note also that C and V have the same block-diagonal structure. For the product AK n we re-use the on-shell row(s) of the AK n−1 matrix product computed for the (n − 1)th term.
The (complex) on-shell transition amplitudes U (E) for spin− 1 2 -spin-1 scattering constitutes a 3×3 matrix. Once this matrix is computed in all relevant N N N partial waves, i.e., for J ≤ 17/2 and J ≤ 3 in our case, it is straightforward to obtain the 6×6 spin-scattering matrix for describing the elastic N d scattering cross sections at kinetic energy E Lab in the laboratory frame of reference, see App. D-E.

D. Convergence with respect to N WP
In the limit N WP → ∞, amplitudes computed using the WPCD method approaches results from the standard Faddeev method utilizing a plane-wave basis. This infinite limit cannot be reached in practice and all WPCD predictions that we present are based on solving Eq. (1) in a finite wave-packet basis. To analyze the convergence of predictions with respect to increasing N WP we computed nd scattering phase shifts, shown in Fig. 1 for the doublet and quartet spin-channels in the J Π = 1 2 + and 7 2 + N N N partial waves, respectively, using the Nijmegen-I N N interaction [17]. For this potential there exists published results [2] from a standard Faddeev calculation at E Lab = 13 MeV and this provides a valuable benchmark to ensure the correctness of our implementation. Detailed numerical inspection of the results reveal that we recover standard Faddeev results for all imaginary and real parts of the N N N phase shifts for J Π ≤ 7 2 ± within ∼ 1% using N WP 125 wave packets. We also observe that the magnitude of the imaginary part of the phase shifts is |Im(δ)| 10 −2 degrees for scattering energies below the deuteron breakup threshold. The convergence with increasing N WP is rather slow however, which is to be expected since a packetized basis corresponds to a coarse grained continuum representation across a wide range of energies simultaneously. In Fig. 1 it is nevertheless clear that the WPCD method yields highly accurate scattering phase shifts for E Lab 50 MeV already for N WP 75. See App. E for further information about how we computed phase shifts from the partial-wave scattering amplitudes U .
Predicting scattering observables is more interesting than scattering phase shifts since they are directly comparable to experimental data. To benchmark our WPCD computation of observables we compare with the results 2 from a standard Faddeev calculation [2] of the neutron analyzing power A y (n) at E Lab = 35 MeV using the Nijmegen-I N N potential, see Fig. 2. For the WPCDcalculations we varied the number of wave packets between 50 ≤ N WP ≤ 150. The convergence pattern is very similar to the one we observed for the phase shifts and for N WP 125 we hence claim convergence for this observable. Also in this calculation we included N N N partial-waves with J Π ≤ 17 2 ± and N N channels with J ≤ 3.
To further assess the convergence of the WPCD method with respect to N WP , we study a range of vector (A) and spherical tensor (T ) analyzing powers, spin transfer coefficients (K), and differential cross sections for 50 ≤ N WP ≤ 125 at E Lab = 3, 10, 65 MeV using the well-known chiral Idaho-N3LO interaction [18], see Fig. 3.
With this result we can establish that for most elastic nd scattering observables it is indeed enough to employ N WP 75 wave packets to obtain sufficiently accurate predictions for E Lab 70 MeV. If one can tolerate a WPCD method-error comparable to typical experimental errors of N d scattering data, then even N WP ≈ 50 will be enough for most low-energy N d predictions. Note that the number of wave packets dramatically impact the computational cost of the WPCD calculations since this scales as ∼ N 4 WP . Also, solving for all amplitudes with E Lab 100 MeV is merely ∼ 2 times slower than solving at a single scattering energy.
In Fig. 3 the wave-packet convergence of the tensor analyzing power T 21 stands out and exhibits a noticeable sensitivity to N WP . This observable is known to also depend more strongly on J max = 4 contributions of the N N interaction [2]. Observables that depend on finer details of the nuclear interaction will exhibit a slower convergence with respect to increasing N WP , due to the coarse-graining of WPCD. Fortunately, poorly converging predictions can be identified straightforwardly. We explored various modifications to the wave-packet distributions, e.g., increasing the density of wave-packets in the vicinity of the scattering energy, but this did not lead to any clear improvements.

III. PREDICTING nd-SCATTERING CROSS SECTIONS USING N2LOopt
In this section we present selected low-energy and elastic nd cross sections using the N2LO opt N N interaction and compare with neutron-deuteron (nd) as well as proton-deuteron (pd) data. We can neglect method uncertainties since we employ N WP = 125 wave-packets for all predictions, unless otherwise stated.
The world database of N d scattering cross sections contains mostly pd data from experiments with either polarized or non-polarized proton or deuteron beams. Indeed, nd scattering is difficult to perform. It is challenging to manipulate and focus electrically neutral particles. The neutron itself is unstable and does not make for a suitable target material on its own. Neutron detectors are also less efficient compared charged-particle detectors. Theoretically, we have the opposite situation. It is typically much easier to compute N d scattering cross sections without a Coulomb interaction [36,37]. Fortunately, Coulomb effects are only significant at low energies, e.g., below the deuteron breakup threshold, and for extremal scattering angles. As such, in most kinematic regions pd scattering data can be compared with theoretical N d scattering results without any Coulomb interaction. We will therefore use pd data in case nd data does not exist or is very scarce. To be clear, we do not include any Coulomb effects in our calculations. One can extend the WPCD method to incorporate such effects. Indeed, the challenge of treating a long-range interaction for small momenta is alleviated when using a square-integrable Coulomb wave-packet basis [14].
In Fig. 4 we show our predictions for the total nd scattering cross section with the N2LO opt interaction. The reproduction of experimental nd data is excellent up to E Lab ≈ 70 MeV. At this point we also begin to see a difference between the N WP = 100 and N WP = 125 calculations. At E Lab > 70 MeV, the inclusion of J > 3 N N -channels will have a percent-level effect on the predictions. We also note that E Lab ≈ 70 MeV corresponds to a relative momentum q ≈ 240 MeV of the incident neutron. This translates to a N N scattering energy of 125 MeV in the lab frame. The N2LO opt goodness-of-fit measure , i.e., the χ 2 /N datum with respect to N N scattering data, is ≈ 1 up to 125 MeV scattering energy. As such, it is reasonable to expect a gradual deterioration of the predictive power for E Lab > 70 MeV.
At energies below E Lab ≈ 50 − 100 MeV the effects of N N N interactions are typically smaller [22,27,32,43], and the bulk of low-energy N d scattering observables can be described quite well using only N N interactions. However, there exist a few scattering observables at these low energies that exhibit discrepancies due to missing N N N forces and (or) possibly fine-tuning effects, e.g., low-energy analyzing powers, the high-energy differential cross section minimum, and the nd doublet scattering length. The latter is known to correlate with the triton binding energy via the well-known Phillips line [44]. The nd scattering length can be computed using a boundstate formulation of the Faddeev equations [45] or via numerical extrapolation of the scattering amplitude to q → 0. Unfortunately, this limit is challenging to reach in the WPCD method with the Chebyshev distribution we employ. Resorting to a basis with an increased number wave-packets at small momenta will of course remedy this. However, to simultaneously maintain accurate scattering amplitudes for higher scattering energies will result in a needlessly large basis size. The prediction of the neutron analyzing power A y (n) at E Lab = 10 MeV with N2LO opt is shown in Fig. 5. For comparison we also include WPCD results using the Idaho-N3LO and Nijmegen-I potentials. All three potentials yield virtually the same result and the discrepancy with respect to 2 H(n, n) 2 H data [40] at the c.m. scattering angles θ cm ≈ 120 • , known as the A y -puzzle [30], persists also with N2LO opt . There is some tendency of a slight increase using this latter potential, but this is certainly not significant on an absolute scale. This result reflects that the low-energy interaction in the 3 Pchannels of N2LO opt , to which we know that A y is most sensitive [30], are similar to the ones in Idaho-N3LO and Nijmegen-I. A detailed calculation [32] to very high chiral orders suggests that the inclusion of leading N N N forces does not resolve the A y -puzzle. Instead, there are hints that the A y -puzzle could be resolved with subleading N N N forces [43]. Alternatively, the A y puzzle might vanish in a simultaneous N N +N N N analysis conditioned on N N and N d scattering data and informed by model discrepancies such as the truncation error in effective field theory.
Low-energy nd differential cross-section data is very well reproduced by N2LO opt , see the left panel of Fig. 6, and the results are identical to what is obtained using the Idaho-N3LO potential. At higher energies, however, there are some discrepancies with respect to data and the two employed potentials differ slightly in the vicinity of the cross section minimum. A previous study [31] concluded that the effects of N N N interactions are expected to be particularly noticeable in this angular region. Although the N2LO opt prediction lies marginally closer to the experimental data at θ cm ≈ 120 • the shape of the differential cross section is not correct.
In Fig. 7 we show a range of spin observables for E Lab = 13 − 70 MeV. Overall, N2LO opt and Idaho-N3LO describe the data rather well in this energy re-gion and the two different potentials produce virtually identical results. In the top row of Fig. 7 we present A y (n) for increasing values of E Lab . It is well known that at energies below ∼ 30 MeV nearly all N N interactions fail to describe the data for this observable [46]. As was discussed above, the N2LO opt interaction does not remedy the puzzle. Nevertheless, careful inspection reveals that the predictions for A y at E Lab = 21 MeV fits the data slightly better at large scattering angles when using N2LO opt . Unfortunately, discrepancies with respect to data persists for small scattering angles, i.e., at the minimum value of A y . In the second row of Fig. 7 we present low-energy neutron-to-neutron spin transfer (K) and neutron-to-deuteron correlation (C) observables. Previous studies [53] have found that the K y y spin transfer is most sensitive to the structure of the N N interaction in the 3 S 1 -3 D 1 and 1 P 1 channels. Since Idaho-N3LO and N2LO opt have very similar N N phase shifts below E Lab = 100 MeV in these channels it is not surprising to recover very similar results also for these observables. Of course, the former potential incorporates higher-order long-and short-range physics that modify the off-shell structure of the potential, but this does not appear alter the predictions much. Regarding the tensor analyzing powers presented in the third row, the discrepancy between theory and data for A xx (d) at θ cm ≈ 150 • persists for both potentials. Inclusion of modern N N N forces does not resolve this [43].

IV. SUMMARY AND OUTLOOK
In this work we presented a framework that allows to solve the Faddeev equations for elastic N d scattering using a newly developed code based on the WPCD method. We analyzed the convergence of the WPCD method, applied to chiral potentials, with respect to the number of basis wave packets N WP . We find negligible method errors when using N WP = 125 in the regime E Lab 70 MeV.
We studied different nd scattering observables up to E Lab = 70 MeV with the N2LO opt N N interaction and find a good overall reproduction of the scattering data. However, the A y -puzzle remains unsolved when applying the N2LO opt interaction. Compared to the Nijmegen-I and Idaho-N3LO N N interactions, we detect a minor increase in the maximum of the theoretical predictions of this observable at low energies. For other cross sections we find a good reproduction of experimental data and the results are virtually indistinguishable from the Idaho-N3LO interaction. The N2LO opt interaction prediction for dσ/dΩ at E Lab = 64.5 MeV is slightly closer to the data at the differential cross section minimum. This observable is typically associated with an increased sensitivity to N N N interactions.
N N N interactions in our calculations. Although some observables, that depend sensitively on finer details of the nuclear interaction, require more wave-packets to be accurately resolved one can obtain sufficiently accurate predictions for the vast majority of low-energy N d cross sections with N WP ≈ 50 wave packets, which will help to reduce the computational demands of the calculations to a level that allows a study of N d scattering observ-ables within a statistical analysis. Specifically, in the near term we plan to employ WPCD predictions to sample Bayesian posterior predictive distributions for N d scattering. Work in this direction, using frequentist methods, was also initiated in [54]. In the longer perspective, emulator methods based on perturbation theory [55] or eigenvector continuation [56] promise an efficient method for fast and accurate emulation of scattering observables [57][58][59][60] and will open new ways to systematically incorporate N d scattering observables in the construction and fitting process of next-generation N N and N N N interactions. In this section we present the permutation operator P 123 in a plane-wave partial-wave basis. The permutation operatorP ijk =P ijPjk performs two pairwise interchanges of particles: first j ↔ k followed by i ↔ j. Our derivation follows the steps presented in [29], as well as the notation and convention for the Jacobi momenta p and q. For this section we will use the indexing (ij) to denote the ij-pair system for the sake of clarity, rather than the odd-man-out notation. To complement [29], we use the (12)-subsystem as our initial states upon which the permutation operator acts. One can show that all representations ofP 123 are invariant under change of reference system. Furthermore, one can show [1,29] that for a basis that is antisymmetric under exchange of particles 2 and 3, we have P 123 = P 132 . This allows us to express projections ofP in Eq. (1) simply as P = 1 + 2P 123 .

(A7)
Here, we also defined (in the pair-system (12)) (A8) Note that the exact form of these relations depend on the choice of pair-system. See [29] for a summary of threebody kinematics.
Choosing to conserve (p,q) will restrict the bramomenta of an operator to the right of the permutation operator in the Faddeev equation due to the ensuing delta-functions (used in e.g. [2,29]). Likewise, (p,q) will restrict the ket-momenta of an operator to the left, while (κ, κ ) will restrict p of operators on both sides, and lastly (π, π ) will restrict q on both sides (used in e.g. [61,62]). In this work we have followed [29] and conserve (p,q). From this point on we drop the (12) subscript on momenta and get 12 p q ; L l L m L |pq; LlLm L 23 = ∞ 0 dp dq wherep = |p| andq = |q|. The integral is invariant under rotations, making it proportional to δ L L δ m L m L . This invariance allows us to simply average out m L , giving a factor 1 2L+1 . We now have freedom in the choice of axes. Choosingẑ p and the polar angle ofq to zero, we can simplify a spherical harmonic: By solving the remaining angular integrals we are left with 12 p q ; L l L |pq; where vectors are functions of (p , q , x) and x = cos(φ) (the angle from q to p ). Notice the change in notation of states on the left-hand side as we averaged with respect to m L . Inserting Eq. (A10) and Eq. (A5) back into Eq. (A2) gives 12 p q ; α |P 123 |pq; with the geometrical function G αα (p , q , x) in a form that is straightforward to implement algorithmically, and where we defined and where we used m L = 0. Given that we have S = S , m S = m S , L = L , and m L = m L , we used the orthogonality of Clebsch-Gordan coefficients, to set J = J and m J = m J . Previously, in [61] the angular dependence in G αα (p , q , x) was evaluated separately from the recoupling terms. This allows pre-calculation of the geometric recouplings before doing the angular integration. However, it turns out that keeping the angular dependence as above is both more numerically efficient and stabler with higher l and L [29]. As this function is the most computationally costly part of evaluating the integral of Eq. (A11), we mention some key optimizations one can use.
The simplest and most effective optimization is to calculate G αα (p , q , x) and store it in the computer memory in its entirety. From a computational viewpoint the function is 5-dimensional, which is still storable in the computer memory for the basis sizes and number of quadrature points we typically require.
Regardless of whether prestorage of G αα (p , q , x) is possible, we still wish to speed up the calculation of G αα (p , q , x). To this end, everything before the second sum in Eq. (A12) (i.e. all geometric recoupling) can easily be precalculated to improve computational performance. The second summation can be sped up by prestoring the three Legendre polynomials individually, which is usually still manageable and quite fast. In this section we present the expressions we employ for projecting operators to a wave-packet basis. Although Eq. (4) provides a definition of a three-body wave packet, it is written explicitly for plane-wave states. Identically, and in general, any three-body wave packet |Y ij (where we omit partial-wave indexing) can be defined from continuum states |ω p and |ω q of some operator (e.g.p,ĥ 0 , orĥ 1 ) with Jacobi momenta p ∈ D i and q ∈D j , where N ij is the wave-packet normalization constant, and f (p) andf (q) are the weighting functions defined in Sec. II B. The wave packet can be projected onto the continuum basis straightforwardly, ω p , ω q |Y ij = 1 N ij Di,Dj dp p dq q f (p )f (q ) × ω p , ω q |ω p , ω q where 1 Di (p) is the indicator function. From this it is easy to show that N ij = D iDj , where D i (D j ) is the width of D i (D j ). Note that for energy wave packets the width is expressed in energy as, for example, D i = E i+1 − E i for boundaries E i . The general form of Eq. (B 1) comes in handy for analytical derivation as it applies to both free and scattering wave packets alike. Going back to a FWP basis, a projection of a general N N N operator will look as follows, × p q ; α |Ô|pq; α . (B3)

Two-body free Hamiltonian
For our chosen normalization, the free N N Hamiltonian is given by (B4) Clearly, the free Hamiltonian is also diagonal in the FWP basis. Depending on the choice of wave packet, we get for the N N free Hamiltonian, The N N potentialv in a N N N partial-wave basis reduces to p q ; α |v|pq; α = δ γ γ δ Γ Γ δ(q − q) q 2 p |v n n |p , (B6) where γ denotes all the quantum numbers for the third nucleon relative to the pair system, i.e., γ = {l, j}, Γ = {J , T } denotes the coupled N N N quantum numbers, and the pair-system quantum numbers are jointly referred to as n = {L, S, J, T }. For our predictions we break total T isospin conservation, and the expressions below must be modified in an obvious way. We obtain the N N interaction in the N N N FWP-basis via Eq. (B3), and easily resolving the q-integral, This expression is straightforward to evaluate numerically using quadrature.

The permutation operator
The permutation operatorP 123 in a partial-wave basis, Eq. (A11), can be inserted into Eq. (B3). The deltafunctions are only non-zero when the Jacobi momentap andq fall within the bins D i andD j , which we express using the indicator function. Choosing momentum wave-packets gives The indicator function is discontinuous and an evaluation of Eq.(B8) using Gaussian quadrature in the p and q momenta yields poor convergence with an increasing number of quadrature points. Therefore, we transform the integral over p and q to polar coordinates using the procedure presented in Ref. [62]: system. The result iŝ Following [64], this can be expressed as the sum of two terms,Ĝ wherê are the bound-continuum (BC) and continuumcontinuum (CC) parts of the channel resolvent, respectively. Here we have defined eigenstates {|ψ α p,n , |ψ α p } ofĥ a with eigenenergies n < 0 for n ≤ N b and E p > 0, respectively, and eigenstates {|ψ α q } ofĥ 0 with eigenenergies E q . Equation (B15) can be projected onto a SWP basis {|Z α ij }, as the inner-products ψ α p,n , ψ α q |Z α ij and ψ α p , ψ α q |Z α ij are known analytically (Eq. (B2)), such that Z α i j |Ĝ(E)|Z α ij = δ i i δ j j δ α α R α ij (E) + Q α ij (E) , (B18) where R α ij (E) ≡ Z α ij |R(E)|Z α ij is given by and Q α ij (E) ≡ Z α ij |Q(E)|Z α ij is given by These integrals can be solved analytically and in the case of energy SWPs we get where and whereĥ 0 1 |x j =¯ j |x j . We also denoted the Heaviside step function with Θ. We do not distinguish between i ≶ 0 since this follows automatically from the operator being calculated, i.e. R α ij or Q α ij . In N d scattering there is only one N N bound state, the deuteron, such that there should only be one index i = i d where R α ij = 0, but here we have kept the expressions above general.
for any-dimensional variables x and y. The Neumann series of this equation is written, This series only converges if all so-called Weinberg eigenvalues η i of K satisfy |η i | < 1, and this is by no means guaranteed in nuclear physics. Indeed, analyzing the Weinberg eigenvalues for nuclear interactions reveals the non-perturbative character in many partial-waves, e.g., where we have bound states [65]. Using Padé approximants is a convenient method for resumming the terms of the Neumann series and extrapolating beyond its radius of convergence. See Refs. [35,66] for more details.
In brief, a Padé approximant of a meromorphic function f (z), which is analytic near z = 0, amounts to formulating the ratio of two polynomial functions P N (z) and Q M (z) of degrees N and M , respectively, such that f (z) = a 0 +a 1 z+a 2 z 2 +a 3 z 3 +. . . = P N (z) Q M (z) +O(z N +M +1 ), The advantage of this Padé approximant is that, contrary to a simple polynomial approximation which would only converge within some radius |z| < R, we can now approximate singularities in f (z). Finding the (unique) coefficients of the polynomials P N and Q N amounts to solving a system of polynomial equations. The solutions are effectively obtained by evaluating the following determinants built from the terms in the Neumann series, {a n } N +M For our studies we have only used "diagonal" Padé approximants where we use M = N ≈ 15 to ensure a convergent scattering amplitude.
convention [70] for the scattering plane: (E4) Uniquely identifying the phase shifts and mixing angles requires a convention for the ordering of eigenvectors. Below the deuteron breakup threshold we order the eigenvectors (which can be chosen to be real) to have a dominant and positive diagonal [71]. Above the threshold we will start getting imaginary components and it becomes necessary to use, for example, the continuity of eigenvectors to arrange U correctly to identify phase shifts [72].