NN Scattering and Nuclear Uncertainties

Ab initio calculations in Nuclear physics for atomic nuclei require a specific knowledge of the interactions among their constituents, protons and neutrons. In particular, NN interactions can be constrained down to scale resolutions of Δr ~ 0.6 fm from the study of phase shifts below the pion production threshold. However, this allows for ambiguities and uncertainties which have an impact on finite nuclei, nuclear- and neutron-matter properties. On the other hand the nuclear many body problem is intrinsically difficult and the computational cost increases with numerical precision and number of nucleons. However, it is unclear what the physical precision should be for these calculations. In this contribution we review much of the work done in Granada to encompass both the uncertainties stemming from the NN scattering database in light nuclei such as triton and alpha particle and the numerical precision required by the solution method.


INTRODUCTION
One of the main goals in Theoretical Nuclear Physics for many years has been to achieve a sufficiently accurate ab initio solution of the Nuclear Many Body Problem from a reductionist perspective. Within the present context this means starting with the forces among the hadronic constituents, protons and neutrons, and solving the corresponding quantum mechanical problem. While this has been widely and openly recognized as an extremely difficult problem, it already represents a simplification as compared to the fundamental problem where the constituents are quarks and gluons building the nucleons and the interactions are deduced from the gauge principle in QCD. The nuclear problem schematically comprises two main steps (i) the determination of the basic interactions from spectroscopy and reactions at the few body level and (ii) a precise method of solution of the inferred interactions for the many body problem. The predictive power of the theory corresponds therefore to the relation between the input (nuclear two-, three-, four-body, and so on, forces) and the output nuclear binding energies, form factors and nuclear reactions, and the corresponding uncertainties.
The seminal paper of Yukawa [1] established the first theoretical evidence that the nuclear force has a finite range by the particle exchange mechanism. The first determination of the tensor force and its consequences for the deuteron were analyzed by Bethe [2,3]. The first χ 2 statistical analyzes of NN scattering data below pion production threshold started in the mid fifties [4] (an account up to 1966 can be traced from Arndt and Macgregor [5]). A modified χ 2 method was introduced [6] in order to include data without absolute normalization. The steady increase along the years in the number of scattering data with better precision generated incompatibilities and hence different criteria had to be introduced [7][8][9] to discard inconsistent data. For a comprehensive review up to 1977 see [10][11][12][13]. For a historical presentation before 1989 we recommend Machleidt [14].
Error analysis of NN phase-shifts for several partial waves became first possible when the Nijmegen group [15] carried out a partial wave analysis (PWA) fitting about 4,000 experimental np and pp data, after rejecting about 1,000 inconsistent data with a 3σ criterion. The analysis resulted in a value χ 2 /dof ∼ 1. In the fit the potential was an energy dependent square well of radius 1.4 fm, plus one-pion-exchange (OPE) and charge-dependent (CD) contributions starting at 1.4 fm, and a one-boson-exchange (OBE) piece operating below 2-2.5 fm. Unfortunately, the required energy dependence becomes messy for nuclear structure calculations. In the next decade a variety of NN (energy independent) potentials appeared in the literature fitting a large body of scattering data with χ 2 /dof ∼ 1 [15][16][17][18][19], but surprisingly error estimates on potential parameters were not made. While all these modern potentials share the local OPE and CD tail and include electromagnetic effects, the unknown short range components of these potentials display a variety of forms and shapes: local potentials [16], nonlocal ones with angular momentum dependence [17], energy dependence [15], or momentum dependence [16,18,19]. While in principle p−, L−, and E−non-localities are equivalent onshell (see e.g., Amghar and Desplanques [20] for a proof in a 1/M N expansion) they reflect truly different physical effects and generally one should consider them as independent quantities. Any specific choice results in a bias and hence becomes a source of systematic errors.
Error propagation from nucleon-nucleon data to threeand four-nucleon binding energies was pioneered in Adam et al. [21]. A rudimentary method based on coarse grained NN interactions was proposed [22,23] providing a first guess for error on bindings in nuclei and neutron and nuclear matter. The Granada analysis of the triton using hyper-spherical harmonics method was performed in Navarro Perez et al. [24]. The triton and the alpha particle were analyzed by solving the Faddeev equations for 3 H and the Yakubovsky equations for 4 He in [25], and in ab initio no-core full configuration calculations [26]. Theoretical uncertainties in the elastic nucleondeuteron scattering observables were calculated in Skibinski et al. [27].
While the history of the NN force and its applications to nuclear physics is rather long, uncertainty quantification has not been addressed seriously until recently (see e.g., [28] for a review prefacing a full volume of the ISNET community). There are several reasons why we think that stressing this aspect of the theory may be particularly useful and fruitful. One obvious one is to provide sensible error estimates in the theoretical calculations. The traditional way was to try out several schemes and compare the different results. Another, less obvious reason, is to address the many body nuclear problem within the realistic physical accuracy, rather than the computational accuracy as it has been the customary approach up to now. This applies in particular to the a priori accuracy of the solution of the nuclear many body problem, which may eventually be relaxed as to facilitate calculations not addressed before. However, this may occur at a high price; it is not unthinkable that any realistic attempt to quantify the theoretical uncertainties may end up with a lack of predictive power on the side of the theory.
We distinguish as usual in error analyses two sources of uncertainties: statistical errors stemming from the data uncertainties for a fixed form of the potential, and systematic errors arising from the different most-likely forms of the potentials. Assuming they are independent, the total uncertainty corresponds to adding both uncertainties in quadrature. In what follows it is advantageous to take the viewpoint of considering any of the different potentials as an independent but possibly biased way to determine the scattering amplitudes and/or phaseshifts. Because the biases introduced in all single potential are independent on each other, a randomization of systematic errors makes sense.
A prerequisite for such an analysis is to discern as much as possible between statistical and systematic uncertainties. The former correspond to the proper propagation of the experimental input while the latter is concerned with the model or scheme dependence of the calculation procedure. Systematic errors may include the genuine bias to describe the physics and truncation errors which are related to the approximate way the calculation is carried out. At the present stage, the model bias is the largest source of uncertainty.
After many years of tremendous efforts and steady progress, state of the art calculations suggest that considerable success can be expected if one includes the current knowledge of the two-, three-body forces and a variety of many body techniques are applied. Going beyond four-body forces has never been tried out, partly because of technical difficulties but also because of the appearance of α−clustering, based on the large stability and compactness of the 4 He nucleus, suggests that five body forces are marginal 1 .
As already said, a credible quantification of the accuracy of the theory requires a judicious determination of all sources of error in the final results, including both the experimental information needed to pin down the interactions as well as the convergence of the numerical procedure used to solve the many body problem. Given the formidable computational effort needed to implement accurately many body calculations-even for light nuclei-an a priori determination of the errors induced from input data would very helpful. This would set an useful accuracy goal and a limit beyond which all refinements in the numerics would not improve the theoretical accuracy of the output. The purpose of the present work is to review estimates on such limiting accuracy based on the imperfect knowledge of the basic two body interactions.
Unfortunately, the situation we face in strong interactions in general and in nuclear physics in particular is to compare and validate inaccurate theories on the basis of accurate data. No theoretical predictions outperforming experimental measurements in accuracy are easily found. To make our point and concern more clear let us take for instance the case of nuclear binding energies from a semi-empirical point of view, where a direct reference to nuclear forces is mostly avoided. Bindings are experimentally known to high accuracy, B = 0.01 − 10 KeV, whereas liquid-drop model inspired mass fit formulas yield a lower theoretical accuracy B = 0.6 MeV (see e.g., Toivanen et al. [29] and references therein). This suggests that already within such a simple picture the phenomenological theory is generally not expected to be more accurate in its predictions than experiment. Actually, according to the standard χ 2 /dof ∼ 1 criterion the previous results show that the theory is literally incompatible with data, and thus not even an error analysis based on uncertainty propagation may be undertaken. The situation is presumably less optimistic for the ab initio approach based entirely on the knowledge of (multiparticle) nuclear forces and a skillful solution of the nuclear many body problem. This provides a motivation to quantify the accuracy needed to solve the many body problem.

STATEMENT OF THE PROBLEM
Let us be more specific on the meaning of uncertainty quantification in nuclear physics. From a Hamiltonian describing A-nucleons, H A , with kinetic energy T = A i=1 p 2 i /2M and multi-nucleon forces V nN , where (2) one proceeds to solve the Schrödinger equation In the absence of useful and accurate QCD-ab initio determinations, phenomenological V 2N interactions are adjusted to NN scattering data and the deuteron, 2 H (A = 2), binding energy, while V 3N enter into the 3 H and 3 He (A = 3), bindings, V 4N in 4 He (A = 4), and so on. Thus, the theoretical predictive power flow is expected to be from light to heavy nuclei. For instance, in the case of the binding energy the problem of error propagation based on NN force variations corresponds to The meaning of the variation V NN is a bit subtle, since there are variations which are (scattering) equivalent and hence do not change the scattering observables. We are interested firstly in the NN scattering problem [30]. Quite generally we will consider non-relativistic scattering of two particles with masses m 1 and m 2 where H = H 0 + V and H 0 = p 2 /2µ is the kinetic energy and µ = m 1 m 2 /(m 1 + m 2 ) the reduced mass (we drop "NN" for simplicity). The S-matrix is defined as a boundary condition problem for E ≥ 0 where we have introduced the T-matrix which satisfies the scattering equation in operator form, where in the second equality we write the exact summation of the perturbative series. Other (complex) energy values are defined by analytical continuation. The T-matrix satisfies the reflection (6) and hence the unitarity condition, S(E + iǫ)S(E + iǫ) † = 1, follows also from V = V † in Equation (6). The phase-shift is defined in terms of the eigenvalues of the S-matrix, so that Sϕ α = e 2iδ α ϕ α and for rotational invariant interactions (we neglect spin to ease the notation) the scattering amplitude M(p ′ , p) is given by with Y lm (p) the spherical harmonics and in our convention dσ/d = |M(p ′ , p)| 2 the differential cross section. Any NN unitary transformation, U, transforms the Hamiltonian and hence the potential as V →Ṽ = UHU † − H 0 . For an infinitesimal transformation U = 1 + iη + . . . , where η is a small self-adjoint two-body operator, the scattering equivalent variation corresponds to the change V = i[η, H]. To see the effect on scattering, start with the LS equation in the form T −1 = V −1 − G 0 which upon a variation of the potential produces a variation of the T-matrix T = TV −1 VV −1 T and after some manipulation one gets so that sandwiching this expression between plane waves gives which vanishes in the on-shell limit E k = E k ′ = E and ǫ → 0. Thus, or equivalently for finite unitary transformations, using Equation (7), δ l,H (p) = δ l,UHU † (p). Given this general ambiguity the long lasting problem has been to decide which is the proper representation of the NN interaction based on NN scattering data. This is in essence the so-called inverse scattering problem which has been studied extensively in the past (see e.g., Chadan and Sabatier [31] and Newton [32] for reviews)] and requires additional strong assumptions to fix the particular form of the potential. For instance, assuming a local potential and complete knowledge of the phase-shifts in each partial wave it is possible to determine the solution uniquely provided the binding energies and long distance behavior of the corresponding bound states wave functions allocated by the potential are known. Clearly, these inverse scattering ambiguities have an impact on the solution of the many body problem, as was documented long time ago in nuclear matter [33] and in the triton and alpha particles [34], just to mention two prominent examples (see Srivastava and Sprung [35] for a review).
Much of the arbitrariness is reduced by invoking an underlying theoretical description in terms of hadronic degrees of freedom, which allows to compute V NN ( x) in terms of one-, two-,. . . , pion exchanges. which in turn may be related to the πN scattering process, involving coupling constants for vertex interactions. At present such a picture seems to hold down to NN separations of about the elementary radius, r c = 1.8 fm, below which composite and finite size effects start playing a role That means that, essentially, variations of the NN potential of are restricted at least to V NN ( x) = 0 for r ≥ r c ≈ 1.8 fm.

The Concept of a Potential
In order to properly formulate the uncertainties of the potentials it would be adequate to review first the meaning of a potential in nuclear physics. This is of utmost importance but also intriguing. On the one hand the potential is not an observable but on the other hand to our knowledge it is not practical to carry out ab initio calculations in Nuclear Physics at the hadronic level without potentials. Ultimately, one hopes to be able to provide a direct link between the uncertainties in the input data and propagate them to the output of the many body problem. As said, this is only possible by using non-observable nuclear potentials as intermediate steps.
From a classical (and macroscopic) point of view, potential and force can be measured directly by just determining the separation static energy between two infinitely heavy sources. Such a definition admits a direct extension to the quantum mechanical microscopic case and specifically to the NN interaction assuming interpolating composite local nucleon fields made out of three quarks. In essence, this is the approach followed in recent years in lattice QCD where many of the traditionally assumed features of the NN interaction seem to be confirmed [36][37][38]. A major drawback of this approach is that such a calculation determines the static NN energy which would become a physical observable if nucleons were infinitely heavy. The quantum mechanical problem needs adding kinetic energy contributions. Moreover, the fact that low energy NN scattering provides unnaturally large cross sections corresponds to an extreme fine tuning which is beyond the present lattice capabilities.

The Tensorial Structure
Assuming isospin invariance for the moment, the most general form of the NN interaction can be written as Okubo et al. [39] where p ′ and p denote the final and initial nucleon momenta in the CMS, respectively. Moreover, q = p ′ − p is the momentum transfer, P = (p ′ + p)/2 the average momentum, and S = ( σ 1 + σ 2 )/2 the total spin, with σ 1,2 and τ 1,2 the spin and isospin operators, of nucleon 1 and 2, respectively. The scalar functions appearing in the potential, Equation (11), depend on both initial and final momentum p and p ′ respectively. Because of rotational invariance we may thus form three independent invariants, such as p, p ′ and also q·P (which vanishes on-shell). Transforming to coordinate space in the variable r, conjugate to q, we have where we take P + 1 2 q|V|P − 1 2 q ≡ V(p ′ , p). The case where these functions depend only on the momentum transfer q = p ′ − p corresponds in coordinate space to a local potential, V(r, P) = V(r). Local potentials are appealing because they provide physical insight besides being directly manageable by means of a Schrödinger equation in configuration space. Moreover, attaching a field theoretical interpretation to the interaction, locality must be satisfied by heavy and point-like elementary nucleons which act as static sources, so that in this case the potential becomes the static energy between nucleons which is an unique observable defined by where we assume M N ≫ m π , E. Non-localities are expected to be weak because P/M N ≪ 1, and should have a larger influence at short distances (see e.g., Piarulli et al. [40] for an explicit implementation). The finite mass effects generate some ambiguity in the definition of the potential and, as we will see, are the largest source of uncertainties in nuclear physics. In any case, there is some freedom that can be used advantageously to choose-by means of suitable unitary transformations [41]-a convenient form of the potential to simplify the solution of the two-body problem, and to simplify a particular scheme of the many body problem. We remind, however, that this choice may be a source of bias and hence of systematic uncertainty.
L·ST 12 , (L·S) 2 T 12 (16) and are labeled as T, σ T,tT, τ z,σ τ z, l2T, l2σ T, lsT, and ls2T. The first five were introduced by Wiringa et al. [17]; the following two were included in Navarro Pérez et al. [42] to restrict CD to the 1 S 0 partial wave by following certain linear dependence relations between V T , V σ T , V l2T , and V l2σ T . The last two terms are required for the CD on the 3 P 0 , 3 P 1 , and 3 P 2 partial waves. To incorporate CD on P waves two more operators need to be added to the basis we used previously getting a total of 23 operators O n . As in our previous analysis we set V tT = V τ z = V σ τ z = 0 to exclude CD on the tensor terms and charge asymmetries. To restrict CD to the S and P waves parameters the remaining potential functions must follow The algebraic relation between the operator basis in momentum space and in configuration space is explicitly given in Navarro Perez and Ruiz Arriola [43] and several examples are displayed.

The Long Range Contributions
As mentioned above, the potential becomes an observable within a QFT setup for infinitely heavy hadronic sources. For the finite mass case one may use instead a perturbative matching procedure between a QFT with hadronic (and electro-magnetic fields) fields and the quantum mechanical problem, which should work at sufficiently long distances. The hadronic QFT calculable contribution is separated into two pieces, the strong (pion exchange) piece and the purely EM piece, The CD-OPE potential in the long range part of the interaction is the same as the one used by the Nijmegen group on their 1993 PWA [15] and reads (20) being f the pion coupling constant, σ 1 and σ 2 the single nucleon Pauli matrices, S 1,2 the tensor operator, Y m (r) and T m (r) the usual Yukawa and tensor functions, CD is introduced by the difference between the charged m π ± and neutral m π 0 pion mass by setting V OPE,pp (r) = V m π 0 ,OPE (r), The neutron-proton electromagnetic potential includes only a magnetic moment interaction where µ n and µ p are the neutron and proton magnetic moments, M n the neutron mass, M p the proton one and L · S is the spin orbit operator. The EM terms in the proton-proton channel include one and two photon exchange, vacuum polarization and magnetic moment, where Note that these potentials are only used above r c = 3 fm and thus form factors accounting for the finite size of the nucleon can be set to one. Energy dependence is present through the parameter where k is the center of mass momentum and α the fine structure constant. Table 1 lists the values used for the fundamental constants in our calculations and typically used since the benchmarking Nijmegen analysis.

Short Range Contributions
The short range contributions are fundamentally unknown and, despite some lattice QCD efforts [36][37][38]44], can only be determined indirectly and phenomenologically, mostly from NN scattering. Along the years some experience has been gathered about the size, shape, and range of the potentials in the bulk, at least in configuration space, so that refinements are made by a χ 2 minimization to pp and np scattering data (see below). Besides, the analysis of scattering data allows to obtain information on the lowest distance where the long range contributions can be trusted. We anticipate that they may be assumed to be valid for r c ≥ 1.8 fm when OPE and TPE contributions are included. This coincides a fortiori with the distance above which protons interact by Coulomb force as point-like particles, and also with the typical distance between nucleons in nuclear matter, d = ρ −1/3 = 1.8 fm for ρ = 0.17 fm −3 . Finally, there is the issue on which and how many parameters are needed to describe the short range force in a satisfactory manner. The primary 2013 Granada analysis has been carried out in terms of the so-called coarse grained potentials [45]. The coarse grain procedure samples the interaction with an optimal grain size, corresponding roughly to the reduced de Broglie wavelength r =h/p. For the maximum LAB energy, 350 MeV, this corresponds to r = 0.6 fm. Thus, we do not need to sample the potential functions V i (r) at all points, but rather in a grid of points, V i (r n ) given by r n = n r. We consider the V i (r n ) values as fitting parameters. The particular interpolations between these points are not physically relevant, because shorter scales than r cannot be probed by the scattering process below a maximal p = √ T LAB M N /2 ∼ 2 fm −1 . The number of grid points depends on the cut distance, r c , above which the functional form of the potential is known and corresponds to N = r c / r. Thus, the simplest case corresponds to r c = 1.8 fm and N = 3 grid points for any radial component, V i (r n ), in the operator basis. In the partial wave basis some refinements can be incorporated since the centrifugal barrier limits the sampling points below the barrier in the classically forbidden region, so that the estimate is Fernandez-Soler and Ruiz Arriola [46] and Ruiz Arriola and Ruiz de Elvira [47], where g S and g T are spin and isospin degeneracy factors. The counting of parameters for pp and np [48] yields about 40 "grained" points r n in the fit carried up to a maximum energy T LAB ≤ 350 MeV. This a priori estimate coincides in the bulk with the number of parameters which have traditionally been needed to fit data satisfactorily in the past. The previous argument suggests that including more parameters is not expected to improve significantly the fits to scattering data, but rather increase the correlations among the V i (r n ) parameters.
There are many possible ways to describe the interaction at the "grained" points. The simplest is to consider Dirac delta-shells located at the sampled points [49,50] We refer to Navarro Perez et al. [51] for a pedagogical presentation of coarse grained interactions which solve the Schrödinger equation by a discretized form [49,50] of the variable phase approach of Calogero [52]. This delta-shells decomposition implies a similar one at the partial waves level, so that one may use the partial wave strengths V JS LL ′ (r n ) as fitting parameters. This choice is rather convenient for least squares minimization as the low angular momentum partial wave components of the potential are largely uncorrelated, substantially speeding up the minimum search [53,54]. The transformation matrix from the V i (r n ) to the V JS LL ′ (r n ) basis can be found in Navarro Pérez et al. [42].

PARTIAL WAVE ANALYSIS
The NN scattering amplitude has five independent complex components which are a function of energy and scattering angle [55], We use the three unit vectors (k f and k i are relative final and initial momenta), For this amplitude the total spin S is conserved and in this case the partial wave expansion reads, where S is the unitary coupled channel S-matrix, and the C ′ s are Clebsch-Gordan coefficients, C l,S,J m,m s ,M = lmSM s |JM . The spins of the nucleon pair can be coupled to total spin S = 0, 1 and hence J = L ± 1 for unnatural parity, (−1) L+1 states and J = L for natural parity states. This amplitudes contains all measurable physical information and the relation to observable quantities such as differential cross sections and polarization asymmetries can be found in Hoshizaki [56] and Bystricky et al. [57].
In the Stapp-Ypsilantis-Metropolis (SYM) representation [4] the S-matrix is written in terms of the nuclear-bar phase shifts δ j±1 andǭ j . Dropping the bars for simplicity and denoting the phase shifts as δ J,s l,l ′ , for the singlet (s = 0, l = l ′ = J) and triplet uncoupled (s = 1, l = l ′ = J) channels the S matrix is simply e 2iδ J,s l,l , in the triplet coupled channel (s = 1, l = J ± 1, l ′ = J ± 1) it reads with ǫ J the mixing angle. The partial wave expansion provides an indirect way to find out the range of nuclear forces by truncating the expansion. According to the standard semi-classical argument (see e.g., [58]), for an impact parameter b = (J + 1/2)/p (p is the CM momentum) the no-scattering condition corresponds to b ≥ a, so that |δ J max | ≤ δ J max where maximal angular momentum is provided by J max ≈ pa with a the range of the force. For the Yukawa OPE interaction the exponential fall-off of the potential also means a similar behavior for the phase-shifts, so typically one takes S, P, D, and F waves as active if the condition is J+1/2 ≈ pr c with r c the separation distance.
We will review briefly the basics of scattering from a NN potential for completeness and to provide our notation. Details may be found in standard textbooks on scattering theory (see e.g., [59]). The generalization of the well-known Rayleigh expansion for spin S is where χ SMs is an eigenspinor with spin quantum numbers (S, M s ), and the functions Y lSJM (x) are the couplings of the spherical harmonics with the spinors χ SMs to total angular momentum J, The local (but angular momentum dependent) NN potential described in the previous section conserves spin S and total angular momentum J, but not the orbital angular momentum L. Therefore, the scattering wave function for spin S is expanded as (38) where the reduced radial wave functions u SJ l ′ l (r) satisfy the coupled channel differential equations and the reduced potential is defined as U(r) = 2µV(r). For regular potentials the boundary condition at the origin reads The integration of the equations can advantageously be done using the delta shell representation of the NN potential taking r = 0.6 fm for r ≤ r c (the coarse-grained and unknown part) and r = 0.1 fm for r ≥ r c (the known field theoretical part). The complete set of equations including Coulomb forces is provided in Navarro Pérez et al. [42]. The scattering boundary condition implies a similar asymptotic condition for the reduced radial wave functions. For the uncoupled case, l = J, one has for In the coupled triplet case, S = 1, the four wave functions u l ′ l (r), with l ′ , l = J − 1, J + 1, are coupled in pairs. The pair verifies the coupled equations On the other hand the pair verifies the same coupled equations by changing α → β. This is equivalent to say that the system (44,45) has two linearly independent solutions that we label as α and β solutions. Their asymptotic behavior can be expressed in terms of the eigen phase shifts as, This is known as the Blatt-Biedenharn (BB) parameterization in terms of the eigen phase shifts δ 1j j±1 and ǫ j . These are related to the nuclear-bar phase shifts by the following equations Unless otherwise stated, in this work the phase shifts will always be assumed to be the nuclear-bar ones. The Coulomb force is included exactly by replacing in the previous formulas the Bessel functions j l and y l by Coulomb functions F l and G l [59]. The inclusion of magnetic moments effect is complicated by their 1/r 3 behavior requiring about 1,000 partial waves [42].

STATISTICS
The statistical treatment we follow here is quite standard, and we list for the benefit of the newcomer to the field the main steps to be discussed in the following subsections. We first address the existing scattering data and then we formulate the nature of the problem and the standard χ 2 approach searching for the most likely potential. This requires discriminating between consistent and inconsistent data, something which can be formulated in terms of a self-consistent selection problem. After this, a direct statistically satisfactory result can be deduced and, more importantly, error propagation may legitimately be carried out in terms of the corresponding covariance matrix implementing statistical correlations. This allows in particular to determine the scattering phase-shifts with error bars reflecting directly the experimental uncertainties. More generally, it allows to transport these experimental errors to any observable based on the nucleon-nucleon potential. We will call these the statistical errors.

Scattering Data
Once we have defined the potential model and the scattering formalism we may proceed to determine the potential parameters V i (r n ) from the available np and pp scattering data and from the corresponding scattering observables which are obtained from the scattering amplitude [56,57] (see also Tables 2, 3 below for the notation). The compilation of the existing published data since 1950 till 2013 is described in detail in Navarro Pérez et al. [42] and comprises 8,124 fitting data including 7,709 experimental measurements and 415 normalizations provided by the experimentalists.

Statement of the Problem
The finite amount, precision and limited energy range of the data as well as the many different observables calls for a standard statistical χ 2 -fit analysis [62,63]. This approach is subjected to assumptions and applicability conditions that can only be checked a posteriori in order to guarantee the self-consistency of the analysis. Indeed, scattering experiments deal with counting Poisson statistics and for moderately large number of counts a normal distribution is expected. Thus, one hopes that a satisfactory theoretical description O th i can predict a set of N independent observed data O i given an experimental uncertainty O i as  We use the notation of Hoshizaki [56] and Bystricky et al. [57].
with i = 1, . . ., N and ξ i are independent random normal variables with vanishing mean value ξ i = 0 and unit variance Establishing the validity of Equation (54) is of utmost importance since it provides a basis for the statistical interpretation of the error analysis.

The Least Squares Minimization
If the ξ i are independent normal variables,then ν i=1 ξ 2 i represents a χ 2 distribution with ν degrees of freedom. Thus, under this hypothesis we may consider the standard χ 2 method, which in our case is defined as its estimated uncertainty and O th i (V k (r n )) are the theoretical results which depend on the fitting parameters V k (r n ), the values of the potentials at the sampled points r n . The least squares minimization has always a solution which may be a global or a local minimum, namely whereV k (r n ) the minimizing parameters. Basically, this minimization eliminates N Par parameters from the N Dat data and we are left with ν = N Dat − N Par degrees of freedom. The important aspect here is the statistical significance of the procedure. This can be checked a posteriori by looking at the residuals . According to the assumption underlying the χ 2 -method, the set of variables R 1 , . . ., R N par should be distributed as normal variables, i.e., they should look as N Par variables extracted from a normal distribution N(0, 1). For a finite sample the veracity of this hypothesis can only be established in probabilistic terms, so that we may estimate how likely or unlikely would it be to accept of reject the starting normality assumption. Technically, this can be done in a variety of ways (see e.g., [53,54,64]), but the most popular measure of goodness of a fit is the χ 2 -test which requires that the fit is accepted if with ν = N Dat − N Par . More elaborate tests may be applied and we refer to Navarro Perez et al. [53,54,64] for further details. In practice this means that for N Dat = 8000 and N Par = 50 we should get χ 2 min /ν = 1 ± 0.016 in order to validate Equation (54). Note that this is very different than the loose claims in the literature where χ 2 /ν ≈ 1 qualifies for a good fit, complemented with a visual inspection of the phase shifts. We emphasize that looking similar is not the same as statistical consistency. In fact, a direct fit to the full database with our model gives χ 2 min /ν = 1.41 which is 25σ away from the expected value. This clearly indicates either a bad model, inconsistent data, or both. A statistical measure of the probability that the theory is plausible is given by the p-value; assuming that the normality of residuals is correct it corresponds to the probability of obtaining results at least as extreme as the results actually observed [62,63]. Thus, the probability of having χ 2 min /ν = 1.41 for ν ∼ 7000 is p = 10 −20 , which clearly rules out that the theory describes the data within fluctuations.

Inconsistent vs. Consistent Data
The determination of theoretical uncertainties requires as a prerequisite the compatibility or consistency of all data. This is a strong condition which is not always fulfilled, particularly when the number of data becomes large. Most often, different experiments have different sources of errors and are mutually incompatible. Thus, while any statistical analysis benefits from a large amount of data, a side effect is the proliferation of inconsistent data. In that case it is obvious that no model will be able to simultaneously describe all the data in a satisfactory manner. To appreciate this point more clearly, assume two experiments which yield the measurements O exp1 ± O exp1 and O exp2 ± O exp2 . If the theoretical estimate is O th , we have Minimizing respect to O th we get , in which case we have two inconsistent measurements. The important question is whether both measurements are wrong or just only one. The term wrong here does not necessarily mean an incorrect measurement; it suffices if one or both errors O exp1 and O exp2 are unrealistically small. In case of a discrepancy one may re-analyze the experiment or simply ask the experts, an unfeasible strategy for the experiments performed in the time span 1950-2013 comprising the analysis. The advantage of the statistical method is that, for a large number of experiments, the systematic errors are also randomized and one may rule out some experiments in a kind of majority vote argument.
The case discussed previously corresponds to two different measurements of the same observable, say the differential cross section at the same energy and angle, and the generalization to any number of experiments is straightforward. However, in the case of experiments with close kinematics there is no simple way to decide between inconsistent data unless some continuity and smooth behavior is assumed in order to intertwine the two measurements. Here is where the model enters and statistical methods will never tell us if a given model is correct but rather if the model is inconsistent with the data. This is a kind of circular argument which can only be avoided by looking for models which congregate as many data as possible in a consistent way. Clearly, following this criterion, once one finds a good model, any improvement of the model should describe more data in a statistically significant fashion. The great advantage is that if there are reasons to intertwine theoretically the different measurements of all possible observables one may discuss the data consistency in a generalized way and be able to select between different observables.

Self-Consistent Data Selection
The self-consistent criterion for data selection was proposed by Gross and Stadler [19] and implemented in Navarro Pérez et al. [45]. The way data have been selected proceeds according to the following procedure: 1. Fit the model to all data. If χ 2 /ν < 1 you can stop. If not proceed further. 2. Remove data sets with improbably high or low χ 2 (3σ criterion). 3. Refit parameters for the remaining data. 4. Re-apply 3σ criterion to all data. 5. Repeat until no more data are excluded or recovered.
The effect of the selection criterion with our model is to go from χ 2 /ν| all = 1.41 to χ 2 /ν| selected = 1.05 with a reduction in the number of data from N Data = 8173 to N Data = 6713. While this seems a drastic rejection it is the largest self-consistent fit to date below 350 MeV. For this number of data this is not a minor improvement; in terms of a normality test, it makes the difference in p-value between having p = 10 −20 or p = 0.68.

Fitting Results
The set of 32 scattering observables which we use for the fits comprises a total of about 7000 selected measurements. It is interesting to decompose the contributions to the total χ 2 both in terms of the fitted observables as well as in different energy bins. The separation is carried out explicitly in Tables 2, 3 for pp and np scattering observables respectively and for the latest fit which includes also the pion-nucleon coupling constants [60,61] (see below). As we can see the size of the contributions χ 2 /N are at similar levels for most observables. Note that observables with a considerable larger or smaller χ 2 /N are also observables with a small number of data and therefore larger statistical fluctuations are expected (we remind that for N independent data we expect χ 2 /N ≈ 1 ± √ 2/N. Likewise, we can also break up the contributions in order to see the significance of different energy intervals, see Table 4. We find that, in agreement with the Nijmegen analysis (see [65,66] for comparisons with previous potentials), there is a relatively large degree of uniformity in describing data at different energy bins. We note also that the fit in the low energy region below 2 MeV gives the largest values for χ 2 /N.
From the optimal fitting parameters V α (r n ) with α = 1 S 0 , 3 P 0 , 3 S 1 , 3 D 1 , E 1 , . . . being the different partial waves in a given pp or np channel, we define (λ n ) α = 2µ ab V α (r n ) r which has units of fm −1 and ab = pp, np. In Table 5 we show the corresponding numerical values. It would be nice to see whether something can be said about the nn interaction. However, one remarkable feature of this and similar analyses is the fact that with the exception of S-waves the short distance parameters can be chosen to coincide in the pp and np systems with common partial waves. The fact that to this date it is not possible to do it for S-waves precludes to predict the nn interaction from the combined np and pp fit (see however a theoretical discussion in Calle Cordon et al. [67]). We compare the fit χ 2 /N| fit with the theoretical expectation χ 2 /N| th = 1 ± 2/N.

TABLE 5 |
Fitting delta-shell parameters (λ n ) JS l,l ′ (in fm −1 ) with their errors for all states in the JS channel for a fit with isospin symmetry breaking on the 1 S 0 partial wave parameters only and the pion-nucleon coupling constants f 2 0 , f 2 p , and f 2 c as fitting parameters We take N = 5 equidistant points with r = 0.6 fm. Wave

Covariance Matrix Error Analysis and Statistical Correlations
After the data selection and fitting, error propagation becomes applicable. Here we show the results for the conventional covariance error analysis which assumes small errors and where one first determines the uncertainty in the fitting parameters V i (r n ) which will be labeled generically as λ i for ease of notation 2 .
Expanding around the minimum values,λ i has where the N P × N P error matrix is defined as the inverse of the Hessian matrix evaluated at the minimum The correlation matrix between the fitting parameters λ i and λ j is given by We compute the error of the parameter λ i as Error propagation of an observable depending on the fitting parameters G = G(λ 1 , . . ., λ P ) is computed as The correlation matrix, Equation (63), has been evaluated in Navarro Perez et al. [53,54] where it has been found that for the potentials in the partial wave basis V JS l,l ′ (r n ) the different points r n are largely correlated within a given partial wave, whereas different partial waves are largely uncorrelated. This information allows to substantially speed up the minimum search as movements in the multidimensional space are thus independent and the approaching path to the minimum operates stepwise [53,54].

Phase-Shifts
The first useful application of error propagation regards scattering amplitudes and phase shifts. Extensive tables for the selected values T LAB = 1, 5, 10, 25, 50, 100, 150, 200, 250, 300, 350 MeV have traditionally been presented since the Nijmegen analysis as representative of the fits. These energy values corresponds to a grid of almost equidistant CM momenta p = √ T LAB M N /2 between 0 and 2 fm −1 . For illustration, Figure 1 compares, for low angular momentum, the phase shifts of the primary PWA in Navarro Pérez et al. [42] from a fit with fixed pion coupling constant, f 2 (blue bands), and the most recent ones [60] (red band) from a fit with charge symmetry breaking on the 3 P 0 , 3 P 1 , and 3 P 2 partial waves and in the pion coupling constants f 2 0 , f 2 p , and f 2 c .

DETERMINATION OF YUKAWA COUPLING CONSTANTS
The first determination of the coupling constant was carried out in 1940 by Bethe who obtained the value f 2 = 0.077 − 0.080 from the study of deuteron properties [3] and very close to the currently accepted value (see Table 1). Subsequent determinations based on a variety of processes can be traced from recent compilations [69,70]. A recent historical account has been given by Matsinos [71] where some newer determinations can be consulted according to his own eligibility criterium. For completeness we also quote recent studies based on piondeuteron scattering [72,73] or on the analysis of Roy equations for πN [74] where an upgrade of the corresponding scattering data is considered.
We note that what follows is a brief summary of the results presented in our previous papers where many more details may be found regarding the most influential observables, the dependence on the cut-off radius r c , the inclusion of two-pion exchange contributions or the energy range used in the fit or the evolution with the numerical values and precision along the years [60,61].
The πNN coupling constant is defined as the pion-nucleonnucleon vertex when the three particles are on the mass shell. The corresponding potentials would be There exist four pion nucleon coupling constants, f π 0 pp , −f π 0 nn , f π + pn / √ 2, and f π − np / √ 2 which coincide with f when up and down quark masses are identical and the electron charge is zero. In NN interactions we have access to the combinations, While there is no reason why the pion-nucleon-nucleon coupling constants should be identical in the real world, one expects that the small differences might be pinned down from a sufficiently large number of independent and mutually consistent data. Note that from np and pp analysis we would obtain f 2 p , f 2 0 , and f 2 c we may deduce the nn coupling using the previous equations f n = −f 2 0 /f p . We try to find out how many data would be needed by recalling that electroweak corrections scale with the fine structure constant α = 1/137 and the light quark mass differences. Thus, for the relative change around a mean value. These are naturally at the 1 − 2% level, a small effect. The question is on how many independent measurements N are needed to achieve this desired accuracy. According to the central limit theorem, for N FIGURE 1 | (Color online) Phase shifts obtained from a partial waves analysis to pp and np data and statistical uncertainties. Blue band from Navarro Pérez et al. [42] from a fit with fixed f 2 and orange band [60] from a fit with charge symmetry breaking on the 3 P 0 , 3 P 1 , and 3 P 2 partial waves and in the coupling constants f 2 0 , f 2 p , and f 2 c .
direct independent measurements the relative standard deviation scales as and δg ∼ g for N = 7000 − 10000. We cannot carry out these direct measurements of g but we can proceed indirectly by considering a set of mutually consistent NN scattering measurements The most recent analysis [60,61] based on the Granada-2013 database comprises 6713 published data. This allows: (i) to reduce the error bars, as expected and (ii) to discriminate between the three coupling constants (see Table 6). When charge dependence in 1 S 0 , P waves is allowed one has f 2 p = 0.0761(4), f 2 0 = 0.0790(9), f 2 c = 0.0772(5), The most remarkable consequence is that from the point of view of the strong interaction neutrons interact more strongly than protons.

SYSTEMATIC VS. STATISTICAL ERRORS: THE 6 GRANADA POTENTIALS
Within the phenomenological approach the estimation of systematic errors can be addressed by using different representations of the mid-range function below the separation distance r c while keeping the long range potential and the NN database. To this end we have analyzed 6 different potentials in Navarro Pérez et al. [75] which have been fitted to the same Granada 2013 database and have the same long distance components of the potential. First we have checked that the 6 Granada potentials are statistically acceptable. In fact, as it has been stressed in our previous works [53,54] one can globally slightly enlarge the experimental uncertainties by the so-called Birge factor [76] provided the residuals verify a normality test.
After this re-scaling the p-value becomes 0.68 for a 1σ confidence level and hence all potentials become statistically equivalent.
The results are summarized in Table 7. Thus, the overall spread between the various phenomenological models with χ 2 /dof ∼ 1 provides an estimate of the scale of the systematic uncertainty. A direct way of illustrating quantitatively the situation is by analyzing the corresponding phase shifts in the different analyses. Thus, for each energy and partial wave, one evaluates the phaseshifts δ (1) , . . ., δ (N) for a representative set of high-precision NN potentials V (1) , . . ., V (N) , and computes the average δ and standard deviation as a measure of the systematic uncertainty of the phaseshifts. In Figure 2 we show the results for four different situations. To provide some historical perspective, we show in the upper left panel the averaged phase shifts, i.e., the absolute (meansquare) errors for np partial wave phase shifts due to the different potentials fitting scattering data with χ 2 /dof ∼ 1 [15][16][17][18][19] as a function of the LAB energy, namely (CD Bonn) [78], Nijmegen (Nijm-I and Nijm-II) [15], Argonne AV18 [17], Reid (Reid93) [79], and the covariant spectator model [19]. As one naturally expects the average uncertainties grow with energy and decrease with the relative angular momentum which semi-classically corresponds to probing an impact parameter where p is the CM momentum, p = √ M N E LAB /2, making peripheral waves to be mostly determined from OPE. These analyses stop at the pion production threshold so that one probes distances larger than Note that the bumps or bulges at low energy in 1 S 0 and 3 S 1 channels in the top left panel are due to a unique potential which is an outlier at low energies. In particular, the authors believe that the outlier behavior is due to the use of an interpolating function used to approximate the potential between the values of laboratory energy at which phaseshifts are usually tabulated.
In the upper right panel of Figure 2 we show the errors obtained via the standard covariance-matrix method explained above and including correlations in the fitting parameters for the primary Granada 2013 analysis [45] which corresponds to the DS-OPE potential. Thirdly, in the lower left panel we show the case of the np phase shifts for the 6 Granada potentials [45,75,77]. Finally, in lower right panel we present the uncertainties for all the 7 pre-Granada potentials and the 6 Granada potentials simultaneously. We indicate the partial waves where charge dependence is allowed.  In Navarro Pérez et al. [75] we found similar statistical errors in all the Granada potentials, which are statistically validated with the same Granada-2013 database, i.e., if the phase-shift for potential V (i) in a given partial wave is δ (i) ± δ (i) stat , then However, we also found that the standard deviation of systematic errors obeys δ sys ≡ Std(δ (1) , . . . , δ (6) ) ≫ δ (i) stat . (77) In all the potentials, the tails above r = 3 fm (including CD-OPE and all electromagnetic effects) are the same, thus the discrepancies between the potentials at short distances dominate the uncertainties, rather than the np and pp experimental data themselves. This conclusion holds also when all high quality potentials are considered [75]. This counter-intuitive result relies not only on the specific forms of potentials which treat the midand short-range behavior of the interaction differently but also on the fact that the fits are mainly done to scattering amplitudes rather than to the phase-shifts themselves.

Low Energy Parameters
The effective range expansion was proposed by Bethe [80] in order to provide a model independent characterization of the scattering at low energies where the shape of the potential is largely irrelevant. The extension to higher partial waves reads (see e.g., [81]) 78) where α l is the scattering length, r l the effective range and v i,l the curvature parameters. In the case of coupled channels due to the tensor force one has that S JS = (M JS − i1)(M JS + i1) −1 with (M JS ) † = M JS a hermitian coupled channel matrix (also known as the K-matrix). At the level of partial waves the multipion exchange diagrams generate left hand cuts in the complex s-plane, which arise in addition to the NN elastic right cut and the πNN, 2πNN etc., pion production cuts. At low energies for |p| ≤ m π /2 we have [82] p l+l ′ +1 M JS l,l ′ (p) = −(α −1 ) JS l,l ′ + 1 2 (r) JS l,l ′ p 2 + (v) JS l,l ′ p 4 + . . . (79) which is the coupled channels effective range expansion. While at lowest orders explicit formulas where available in terms of wave functions, larger order and partial waves become rather cumbersome and no practical formula exists. Fortunately, the variable S-matrix approach of Calogero [52] offers a unique way to extract low-energy threshold parameters for a given NN potential which was extended to coupled channels [82] and applied to the Reid93 and NijmII potentials up to J ≤ 5. For the 6 Granada potentials these have also been extracted and we have found that the systematic uncertainties are generally at least an order of magnitude larger than statistical uncertainties [75]. In Table 8 where we provide the low energy parameters for (J ≤ 2).

Low Energy Constants
Alternatively, one may use effective interactions derived from a low momentum interaction where the coefficients can be The central value and statistical error bars are given on the first line of each partial wave and correspond to the mean and standard deviation of a population of 1020 parameters calculated with the Monte Carlo family of potential parameters described in Navarro Pérez et al. [83] using the DS-OPE potential [42,45]. The second line quotes the systematic uncertainties, the central value and error bars correspond to the mean and standard deviation of the 9 realistic potentials NijmII [16], Reid93 [16], AV18 [17], DS-OPE [42,45], DS-χTPE [48,77], Gauss-OPE [53], Gauss-χTPE, DS-BO, and Gauss-BO. For each partial wave we show the scattering length α and the effective range r 0 , both in fm l+l ′ +1 , as well as the curvature parameters v 2 in fm l+l ′ +3 , v 3 in fm l+l ′ +5 , and v 4 in fm l+l ′ +5 . For the coupled channels we use the nuclear bar representation of the S matrix. Uncertainties smaller than 10 −3 are not quoted.
identified with the phenomenological counter-terms of chiral effective field theory. To obtain such counter-terms we express the momentum space NN potential in the partial wave basis and use the Taylor expansion of the spherical Bessel function to get an expansion for the potential in each partial wave. Keeping terms up to fourth order O(p 4 , p ′4 , p 3 p ′ , pp ′3 , p 2 p ′2 ) corresponds to keeping only S-, P-, and D-waves along with S-D and P-F mixing parameters. Using the normalization and spectroscopic notation of Epelbaum et al. [84] one gets v JS 00 (p ′ , p) = C JS 00 + C JS 00 (p 2 + p ′2 )+ D 1 00 and each counter-term can be expressed as a radial momentum of the NN potential in a specific partial wave. Different methods have been proposed to quantify some of the uncertainties in these quantities [85,86]. Using the statistical uncertainties method and the corresponding systematic error estimates [87], the results are summarized in Table 9 for the 6 Granada potentials.

Scale Dependence and Correlations
While one normally uses a fixed value for the maximum energy in the fits (which in most NN studies has been 350 MeV), one may analyze the consequences of varying this fitting energy [88]. Denoting as the (running) maximal momentum it is clear that the fitting potential will change as is varied. Actually, these parameters may be mapped [54] into the so-called counter-terms which characterize the effective theories at small momenta [89]. We determined the two-body Skyrme force parameters arising from the NN interaction as a function of the maximal momentum in the fit. We found general agreement with the so-called V lowk interactions based on high quality potentials after high energy components have been integrated out [90,91].
In line with our remarks in section 5.7 let us note that, one major outcome of Navarro Pérez et al. [54] has been the fact that the counter-terms corresponding to volume integrals including OPE above 3 fm are weakly correlated, whereas those including OPE+TPE above 1.8 fm have larger but still moderate correlations. Thus, counter-terms in the partial waves basis would be efficient fitting parameters, unlike in the cartesian basis. As we have already discussed, using uncorrelated fitting parameters has the practical consequence of reducing the computational determination of the least squares minimization.

CHIRAL VS. NON-CHIRAL POTENTIALS
In common with the analysis presented in the previous sections, much of the early work on phase-shift analysis was undertaken long before the advent of QCD, so the NN potentials were at most considered to be derivable from Quantum Field Theory in purely hadronic terms. This implies in particular the One-Pion-Exchange potential, which has survived over the years, and the Two-Pion-Exchange which has been changing depending on the computational scheme since the first attempts in the early 50's (see e.g., Machleidt [14] for a historical review, in particular about the meson exchange picture).   After the appearance of QCD as a fundamental theory of strong interactions there emerged dedicated studies on the underlying quark dynamics in terms of quark cluster models, particularly concerning the origin of the nuclear core (see e.g., [92][93][94] and references therein). Despite the numerous attempts it is fair to say that these investigations did provide some microscopic and quantitative understanding of the short range components of the interaction but did not offer an alternative to the conventional partial wave analysis. Current QCD potentials determined on the lattice [36][37][38]44], are still less precise than phenomenological ones.
In the early 90's Weinberg [95] (see e.g., [96][97][98] for comprehensive reviews and references therein) proposed an Effective Field Theory (EFT) approach to NN scattering based on chiral symmetry directly inspired by QCD features, where the spontaneous breakdown of chiral symmetry underlies the would-be Goldstone boson nature of the pion. As compared to the phenomenological approaches, the attractive pattern of such an EFT was also the natural hierarchy of n-body forces and the possibility of making an a priori estimate of the systematic uncertainties in terms of a power counting to different orders. This happened at about the time when the phenomenological approach harvested its great success when the Nijmegen group obtained for the first time a statistically acceptable χ 2 /ν ∼ 1 by fitting and selecting np+pp scattering data. Comprehensive fits to data with chiral interactions have been made using the N2LO chiral potentials [99] to the Nijmegen database [15] for pp [100] and for pp+np [101] and the N3LO chiral potential [102] to the enlarged database [18] for np [103]. The newest generation of chiral potentials have already provided fits to the Granada-2013 database [40, 48, 77, 104-108].

Statistical Issues
Very recently chiral potentials to sixth order in the chiral expansion have been been claimed by the Bochum group to outperform the non-chiral potentials on the basis of the Granada-2013 database [107]. This was a major achievement of the chiral approach (see also [108] for a momentum space approach of the Idaho-Salamanca group). Another great advantage of the chiral approach is that the number of fitting parameters is substantially smaller than in the phenomenological approach. In no case, however, have the authors taken seriously the available statistical tests to verify a posteriori the normality of residuals.
Within the uncertainty quantification context, a critical analysis with an eye on the future developments has been put forward in Ruiz Arriola et al. [109] and Navarro Perez and Ruiz Arriola [43]. It has been suggested that a further order in the expansion, namely N5LO, might quite likely achieve the desired statistical consistency. At the present state, however, there are still some pending, hopefully manageable, issues which need to be resolved before the validation of the chiral approach to NN scattering can be declared without reservations.

The Chiral Tensorial Structure
For instance, the tensorial structure of the force requires phenomenologically that all allowed NN components should contribute to some extent to the total NN potential. Chiral perturbation theory proposes a hierarchy among the different components so that the chiral W Q component vanishes to N4LO, unlike all the phenomenological analyses so far [43]. In addition, the number of independent parameters in a scheme where W Q would be non-vanishing becomes comparable to the phenomenological potentials.

Peripheral Waves
One of the reasons why the coupling constants discussed in section 6 can be pinned down so accurately [60,61] is given by the fact that long distant physics is rather well-determined. From that point of view one expects that peripheral waves are rather sensitive to the shape of the potential and hence become independent of the short range components. This also provides a method to validate other analyses and in particular chiral potentials. A very vivid way of presenting the discrepancy is by comparing the phase-shifts in terms of the impact parameter variable [110] (see Equation 74) for every partial wave which provides a measure of the discrepancy with respect to a set of phase-shifts (see Figure 2 for a plot of different sets). The conclusion of Simo et al. [110] is quite unequivocal: In the range 2 fm ≤ b ≤ 5 fm the δ N4LO differ by more than 3σ when compared to the primary Granada 2013 analysis for F, G, and H waves, and become 1σ compatible with the spread of the 13 high quality potentials.

Perturbation Theory for Higher Partial Waves
The long distance character of chiral potentials suggests that one may determine the high peripheral partial waves in perturbation theory, as done explicitly in Entem et al. [111]. Actually, the low energy parameters discussed above in section 8.1 probe the longest distance features of a given partial wave. Going to N2LO one sees that, while there is some rough agreement between the perturbative and the full low energy parameters, the detailed comparison including both statistical and systematic errors do not agree. Using the perturbative version of the variable phase approach, a perturbative evaluation [43] in the context of chiral TPE (N2LO in the chiral expansion) was also undertaken and shown not to converge to the exact result within uncertainties, even at the largest angular momenta and hence for the most peripheral waves.

Coarse Graining Chiral Potentials
Chiral potentials can be combined with coarse graining in a statistically consistent way [48,48,77,104]. This allows for a reduction of parameters to about 30 since the separation distance can be made as small as r c = 1.8 fm without spoiling the statistical analysis. This approach assumes the chiral power counting for the potential above r c but not in the coarse grained region so that the all the potential components (including the chirally missing W Q ) are non-vanishing, and taking f 2 = 0.0075 has provided natural values for the chiral constants (c 1 , c 3 , c 4 ) = (−0.41 ± 1.08, −4.66 ± 0.60, 4.31 ± 0.17)GeV −1 for T LAB ≤ 350MeV [48,77]. In contrast, the canonical (Weinberg) power counting scheme applies to the full potential and only to at least N5LO provides all non-vanishing tensorial components (W Q = 0 at N4LO), in which case the number of parameters becomes comparable with the phenomenological approach. As emphasized in Navarro Perez and Ruiz Arriola [43], the end of the chiral roadmap in NN scattering based on the power counting will definitely occur when such a scheme becomes reliable enough to select and fit scattering data, without explicit reference to the phenomenological approach.

BINDING IN LIGHT NUCLEI: ERROR PROPAGATION
Much of the previous analysis may be used to analyze the impact of NN scattering uncertainties to binding energies. A precursor of this type of calculations was carried out in Adam et al. [21] where estimates on binding uncertainties were carried out using a statistical regularization of phases and a direct solution of the inverse scattering problem.

On-Shell vs. Off-Shell
NN Scattering data describe only the behavior of nucleons on-shell, i.e., with E p = p 2 + M 2 in the relativistic case. However, nuclear structure calculations usually need also the corresponding off-shell components so that when going from the NN scattering data to the binding energy calculation some extra information would be needed [35]. This ambiguity can be used in fact to our benefit, since ideally one would determine the offshellness from the determination of the finite nuclei properties. The successful attempts by Vary et al. are a good demonstration of that [112,113]

Computational vs. Physical Precision
Let us review the sources of numerical precision in the solution of the quantum-mechanical problem. In the simplest NN case, where we usually solve numerically the two-body Schrödinger equation, the precision is fixed by the precision in the wave function. In the positive energy situation corresponding to a scattering state we are rather interested in the determination of the scattering phase-shifts.
Within the few-body community there has been a trend to determine the quantum mechanical solution with an increasing pre-defined precision, say, a 1%. This is a pure conventional precision which has been a goal per se and, of course, good precision is not disturbing provided the computational cost does not scale up to an unbearable limit where the calculation becomes unfeasible. However, this does not correspond to the physical precision where all necessary effects are taken into account and which determines in fact the predictive power of the theory.

Monte Carlo Method
The normality property of the residuals has been exploited to extract the effective interaction parameters and corresponding counter-terms [54] and to replicate via Monte Carlo bootstrap simulation as a mean to gather more robust information on the uncertainty characteristics of fitting parameters [83]. We stress that the verification of normality, Equation (54), is essential for a meaningful propagation of the statistical error, since the uncertainty inherited from the fitted scattering data O exp i corresponds to a genuine statistical fluctuation. This allows to determine the 1σ error of the parameters p = p 0 ± p stat and hence the error in the potential (83) which generates in turn the error in the NN phase-shifs δ = δ(p 0 ) ± δ stat and mixing angles. Once the NN-potential is determined the few body problem can be solved for the binding energy, where Direct methods to determine p stat , V stat NN and E stat A proceed either by the standard error matrix or Monte Carlo methods (see e.g., [68]). In Navarro Pérez et al. [83] we have shown that the latter method is more convenient for large number of fitting parameters (typically N P = 40 − 60), and consists of generating a sufficiently large sample drawn from a multivariate normal probability distribution P(p 1 , p 2 , . . . , p P ) = 1 where E ij = (∂ 2 χ/∂p i ∂p j ) −1 is the error matrix. We generate M samples p α ∈ P with α = 1, . . ., M, and compute V NN (p α ) from which the corresponding scattering phase shifts δ(p α ) and binding energies E A (p α ) can be determined. Of course, one drawback of the MonteCarlo propagation method is that the object function, in this case the energy, needs to be evaluated a sufficiently large number of times which may be unduly time consuming. An analysis of statistical errors at the phase shift level shows that M = 25 may be sufficient to reproduce consistently the covariance matrix uncertainties from the MonteCarlo method.

The Deuteron
The deuteron is the simplest bound nuclear np system for which the theory has long been developed [114]. Its quantum numbers J P = 1 + correspond to the coupled 3 S 1 − 3 D 1 channel with reduced wave functions u(r) and w(r) respectively, so that we solve the bound state problem with E d = −B d = −γ 2 /2µ np , i.e., with p = iγ . At long distances For normalized states we list in Table 10 the asymptotic D/S ratio η, asymptotic S-wave amplitude A S , mean squared matter radius r m , quadrupole moment Q D , D-wave probability P D and inverse matter radius r −1 for some high quality potentials compared with two Granada potentials, DS-OPE [45], DS-TPE [77]. The PWA analysis indeed uses its binding energy as a fitting parameter, so that the quoted uncertainties are purely statistical. Unlike r m , Q D , or P D which require (small) meson exchange currents corrections before being compared to experimental data, A S and η are purely hadronic. As we see, both the DS-OPE [77] DS-TPE [77] provide smaller uncertainties than the experimental/recommended values for A S and η. To our knowledge, this is an unprecedented situation in Nuclear Physics. Similar trends are also observed for the corresponding deuteron charge, magnetic and quadrupole form factors (see e.g., [121] for a review) where DS-OPE [45] and DS-TPE [77,122] generate tiny uncertainties and offer an opportunity to discriminate meson exchange currents contributions.

Binding Energies for A = 3,4 Systems
The primary Granada DS-OPE potential which was used to fit and select np+pp scattering data uses Dirac delta-shells which are too singular in configuration space or have too long momentum tails, for instance in the deuteron [26], to be handled in few body calculations. Actually, this was the reason to design smooth SOG (Sum of Gaussian) potentials [53,75] referenced in section 7.
In Navarro Perez et al. [24] the triton binding energy was evaluated for the SOG-OPE Granada potential using the hyperspherical harmonics method with M ∼ 200 MonteCarlo replicas, and statistical distributions where also obtained yielding E t = 12 KeV. One motivation for such a calculation was to determine if the computational accuracy was unnecessarily better than the statistical accuracy inherited from the NN scattering data. Our points are illustrated in Table 11 from Navarro Perez et al. [24] where the numerical convergence regarding the number of partial waves is displayed. The error estimate clearly marks where the accuracy of the numerical calculation is larger than the physical accuracy.
The statistical uncertainty of experimental NN scattering data have also been propagated into the binding energy of 3 H and 4 He using the no-core full configuration method in a sufficiently large harmonic oscillator basis. The error analysis [26] Table 12 where we see a systematic underbinding with respect to the experimental values. A popular interpretation of this disagreement suggests that the influence of three-and four-body forces has been neglected. However, the contribution of three body forces depends on the definition of two body forces as we will discuss next.

The Tjon Line
Much of the error analysis which can and has been carried out in Nuclear Physics is probably best exemplified by the so called Tjon  [24] in the number of channels, N c , classified according to the orbital angular momentum of the pair L Pair and the spectator l spectator in the triton as the number of total accumulated channels, N Total , is increased. The potential used was Monte Carlo generated. A horizontal line is drawn when the change in Et is smaller than the statistical uncertainty Bt = 15(1) keV. We list binding energy E d , asymptotic D/S ratio η, asymptotic S-wave amplitude A S , mean squared matter radius rm, quadrupole moment Q D , D-wave probability P D , and inverse matter radius r -1 . line [34], a linear but empirical correlation between the triton and α-particle binding energies of the form where a, c depend on a family of NN potentials which have the same NN scattering phase shifts and deuteron properties. Thus, the slope may be schematically be written as a = (∂B α /∂B t )| B d . This empirical feature [124,125] comparing between phaseequivalent potentials has been corroborated by many calculations ever since [123,126,127]. It is remarkable that such a simple property has no obvious explanation. One clue would be the fact that the deuteron binding energy, B d = 2.2 MeV, is small compared to the triton and alpha bindings [128]. For small B d the alpha binding energy then would scale as B α = aB t + bB d + O(B 2 d ). The points along this line in the plane (B t , B α ) correspond to potentials with the same phase-shifts, verifying B α = a B t The points along a perpendicular line, B α = −1/a B t should correspond to potentials with very different phase-shifts. In particular the difference may be generated by a unitary transformation of the NN potential, V 2 → UV 2 U † , so that the bindings depend on U but the coefficients a and b do not depend on U [123]. On the other hand, a unitary transformation of the two-body potential implies a change in multi-nucleon forces, V 3 , V 4 , etc. and, one may actually fit E t with a suitable V 3 and E α with a suitable V 4 yielding for V 4 = 0 in the socalled on-shell limit the formula B α = 4B t − 3B d which works well [129,130].
Phase equivalent interactions produce a Tjon slope which is typically about B α / B t ∼ 5 − 6 both in the Faddeev-Yakubovsky [126] and in the no-core shell model [131]. For the Faddeev-Yakubovsky solutions of 3 H-4 He the results from five high quality potentials, i.e., with χ 2 /ν ∼ 1 at their time and the Granada SOG-OPE, in Table 12 give The extrapolation predicts the experimental binding of the alpha particle within uncertainties [25], since so that B α | stat ∼ 1MeV. Interestingly, this suggests a marginal effect of four body forces, for which independent estimates using approximate wave functions [132] give similar numbers, B α | 4N ∼ −100 KeV (see also Epelbaum [133] for a chiral scheme where this is argued to overestimate the result). Thus, we see that since B α | 4N ∼ B stat α the four-body force might be unobservable. While this is good news from the theoretical point of view, more detailed calculations might be needed to confirm this feature. Finally, let us also mention that along these lines, theoretical uncertainties of the elastic nucleon-deuteron scattering observables have been undertaken [27].

Moshinsky-Skyrme Parameters
Power expansions in momentum space of effective interactions were introduced by Moshinsky [134] and Skyrme [135] to provide significant simplifications to the nuclear many body problem in comparison with the ab initio approach, in which it is customary to employ phenomenological interactions fitted to NN scattering data to solve the nuclear many body problem. As a consequence of such simplifications effective interactions, also called Skyrme forces, have been extensively used in mean field calculations [136][137][138][139]. Within this framework the effective force is deduced from the elementary NN interaction and encodes the relevant physical properties in terms of a small set of parameters. However, there is not a unique determination of the Skyrme force and different fitting strategies result in different effective potentials (see e.g., [140] and [141]). This diversity of effective interactions within the various available schemes signals a source of statistical and systematic uncertainties that remain to be quantified. Fortunately the parameters determining a Skyrme force can be extracted from phenomenological interactions [88,142] and uncertainties can be propagated accordingly [54]. At the two body level the Moshinsky-Skyrme potential in momentum representation reads where P σ = (1 + σ 1 · σ 2 )/2 is the spin exchange operator with P σ = −1 for spin singlet S = 0 and P σ = 1 for spin triplet S = 1 states. These parameters correspond to radial moments of volume integrals of the potentials ∞ 0 d 3 xr n V i (r) which are increasingly insensitive to short distances.
As mentioned above different nuclear data can be used to constrain the Skyrme potential. The usual approach is to fit parameters of Equation (90) to doubly closed shell nuclei and nuclear matter saturation properties [136][137][138][139]. In Ruiz Arriola [142] the parameters were determined from just NN threshold properties such as scattering lengths, effective ranges and volumes without explicitly taking into account the finite range of the NN interaction; while in Navarro Perez et al. [88] the parameters were computed directly from a local interaction in coordinate space that reproduces NN elastic scattering data. In Navarro Pérez et al. [54] the latter approach was used to propagate statistical uncertainties into the Skyrme parameters. The quantification of the systematic uncertainties, which arise from the different representations of the NN interaction was discussed in Navarro Perez et al. [87]. The results, summarized in Table 13 clearly show, again, the dominance of systematic vs statistical errors.

Error Estimates for Heavy Nuclei and Nuclear Matter
Within the Skyrme effective interactions approach one can find a simple estimate of systematic errors due to the two body interaction uncertainty using (for a review see [139]) For nuclear matter at saturation, ρ 0 = 0.17fm −3 , our t 0 = 75MeV fm 3 implies We may implement finite size effects in light-heavy nuclei by using a Fermi-type shape for the matter density with R = r 0 A 1 3 , r 0 = 1.1 fm and a = 0.7 fm, Normalizing to the total number of particles A = d 3 xρ(x) we get values in the range depending on the value of A for 4 ≤ A ≤ 208.

COARSE GRAINED POTENTIAL RESULTS
Besides the aspect of uncertainty quantification which is the focus of the present work, we believe that the very idea of coarse graining proves useful in nuclear physics. This requires that special methods have to be developed for delta-shells interaction, which in our view are the most flexible ones which allow for selecting and fitting the largest NN database to date, but cannot be plugged directly in conventional computing codes dealing with nuclear structure and reactions, and hence smooth potentials (such as the SOG-Granada type potentials) need to be defined after the data selection process. This is similar to what happened with the energy dependence needed by the Nijmegen group which also led to subsequent high quality interactions. We discuss here some simple examples where delta-shells may be used directly.

Repulsive vs. Structural Core
Besides the well-accepted OPE mechanism for long distances and the mid-range attraction which is needed for nuclear binding, one of the traditional and well-accepted properties of the nuclear potential is the existence of a nuclear strongly repulsive core at about 0.5 fm. While this feature guarantees the stability of nuclei and nuclear matter against collapse it also complicates the solution of the many body problem, since the relative NN wave function must vanish below the core, therefore introducing a very strong short range correlation. At a practical level the existence of the core implies a vanishing of the wave function at about the core location, but something else is needed to determine the wave function below the core radius. The question is whether the repulsive core is indispensable from the analysis of collision experiments. However, in order to resolve the core in NN elastic scattering one needs a wavelength which corresponds to energies where there is a substantial in-elasticity and hence a complex optical potential is needed in order to deal with the absorption due to inelastic processes such as NN → NNπ. This point has been analyzed in Fernandez-Soler and Ruiz Arriola [46] and it has been  Errors quoted for each potential are statistical; errors in the last column are systematic and correspond to the sample standard deviation of the six previous columns. See main text for details on the calculation of systematic errors. Units are: t 0 in MeVfm 3 , t 1 , t 2 , W 0 , t U , t T in MeVfm 5 , and x 0 , x 1 , x 2 are dimensionless.
found that there exist two solutions, one corresponding to the usual repulsive core and the other one related to the socalled structural core, reminiscent of the composite character of the nucleon.

Coarse Graining Short Range Correlations
The Bethe-Goldstone equation [143,144] has been a way to describe short range correlations between nucleons inside the nucleus. In the nuclear medium the interaction produces no scattering due to the Pauli principle. Instead the relative wave function of a pair is modified in presence of the two-body interaction, generating high-momentum components above the Fermi momentum, p > p F . Using the delta-shell potential allows to simplify the problem of computing these high momentum components arising in an interacting nucleon pair in nuclear matter. This coarse graining of the Bethe-Goldstone equation has been explored in Ruiz Simo et al. [145,146] for back-to-back nucleons, with total center of mass momentum equal to zero. The formalism still has to be extended to other values of the center of mass.

Error Analysis of Nuclear Matrix Elements
The expected errors of harmonic oscillator nuclear matrix elements coming from the uncertainty on the NN interaction have been estimated in Amaro et al. [147] for the coarse grained (GR) interaction fitted to NN scattering data, with several prescriptions for the long-part of the interaction, including one pion exchange and chiral two-pion exchange interactions.

Shell Model Estimates
In a previous calculation [51], we showed how our approach is competitive not only as a way of determining the phase shifts but also compared to more sophisticated approaches to Nuclear Structure [148]. We computed the ground state energy of several closed-shell nuclei by using oscillator wave functions. In the case of 4 He, 16 O, and 40 Ca nuclei, our calculation reproduces the experiment at the 20 − 30%-level provided the phase-shifts are fitted up to energy E ≤ 100MeV [51]. This is a tolerable accuracy as we just intend to make a first estimate on the systematic uncertainties and then compute the change in the binding energy. For the A = 3, 4 nuclei we use the simple formulas, B( 4 He) = V 2 4 He = 6 1s| 1 2 V1 S 0 + V3 S 1 |1s , (96) where |1s is the Harmonic oscillator relative wave function with the corresponding oscillator parameter b fixed to reproduce the physical charge radius. The factors in front of the matrix elements are Talmi-Moshinsky coefficients corresponding in this particular case to the number of pairs interacting through a relative s-wave. Errors in the potential V are computed by adding individual contributions ( λ n ) JS l,l ′ in quadrature. By propagating the potential errors to Equation (95) we find B(3) 3 = 0.07 − 0.085MeV (97) depending on the fitting cut-off LAB energy, 100-350 MeV respectively, overestimating the Faddeev estimates given above. For the α−particle Equation (96) These systematic estimates using shell model are of the same order to the ones obtained above in the Skyrme interaction.

OUTLOOK
Despite the many years elapsed since the first NN partial wave analysis in 1957 and the huge theoretical and experimental efforts carried out, the nuclear force is poorly known still where it is most needed, namely in the mid-range regime which is relevant for ab initio calculation of nuclear binding energies. This is the explanation behind the relatively large uncertainties found in large scale calculations. During many years there has been a conformist attitude regarding these uncertainties, and in most papers a purely computational approach has prevailed, validating theoretical frameworks just on their numerical performance. Only in recent years the issue of uncertainties has been taken seriously, as it is actually the key to establish the predictive power of the theory. Clearly, the level of ambiguity we are dealing with in the evaluation of nuclear uncertainties of all sorts, statistical, systematic, and computational requires a rigorous treatment. In this work we have reviewed this topic from the perspective of the impact of the Granada NN database on the determination of the NN force and its consequences on nuclear binding. The main theoretical obstacle has to do with the great difficulty in providing a unique definition of the nuclear potential just from data. Quantum field theory at the hadronic level implies the existence of a long range interaction dominated by pion exchanges as the lightest particles and reduces the ambiguity. Lattice calculations of potentials may identify them with static energies assuming heavy quark-composite sources but their accuracy is at present not satisfactory. Chiral perturbation theory provides in addition several schemes based on a power counting which, while not fully satisfactory, may be and have been implemented in the NN sector and extended to multinucleon forces. The consistency among chiral multi-nucleon forces is theoretically very appealing and the use of potentials is possibly the only practical path toward a satisfactory solution of the nuclear many body problem. It should be stressed that the EFT point of view is the most suitable one since in principle one gets rid of the model dependence with a priori uncertainty estimates. However, a more detailed analysis reveals that there are issues regarding the necessary regularization of the theory, which effectively model the mid-range regime of the NN interaction. Moreover, the indispensability of the chiral scheme for NN scattering data remains to be proven, not to speak about its suitability for fitting and selecting a NN database itself. At a phenomenological level at the present stage the determination of the NN interaction below 1.8 fm (up to a phase equivalent unitary transformation) remains so far connected to a combination of an abundance of data in a variety of kinematics and observables with the corresponding experimental errors.
In our view, this unfortunate situation on the side of the hadronic theory will likely not necessarily improve neither with more and better experimental measurements nor with larger computational facilities, but with a better understanding on the essence of hadronic interactions and their range of applicability.