Minimally non-local nucleon-nucleon potentials with chiral two-pion exchange including ∆ ’s

We construct a coordinate-space chiral potential, including ∆-isobar intermediate states in its two-pion-exchange component up to order Q 3 ( Q denotes generically the low momentum scale). The contact interactions entering at next-to-leading and next-to-next-to-next-to-leading orders ( Q 2 and Q 4 , respectively) are rearranged by Fierz transformations to yield terms at most quadratic in the relative momentum operator of the two nucleons. The low-energy constants multiplying these contact interactions are ﬁtted to the 2013 Granada database, consisting of 2309 pp and 2982 np data (including, respectively, 148 and 218 normalizations) in the laboratory-energy range 0–300 MeV. For the total 5291 pp and np data in this range, we obtain a χ 2 /datum of roughly 1.3 for a set of three models characterized by long- and short-range cutoﬀs, R L and R S respectively, ranging from ( R L , R S ) = (1 . 2 , 0 . 8) fm down to (0 . 8 , 0 . 6) fm. The long-range (short-range) cutoﬀ regularizes the one- and two-pion exchange (contact) part of the potential.


I. INTRODUCTION
The nucleon-nucleon (N N ) interaction is a basic building block in nuclear physics as it makes it possible to describe nuclear structure and nuclear reactions. If the forces were known accurately and precisely, the nuclear many-body problem would become a large-scale computation where precision and accuracy are defined in terms of the preferred numerical method. However, the lack of direct knowledge of the forces among constituents at separation distances relevant for nuclear structure and reactions drastically changes the rules of the game. Indeed, the use of a large but finite body of scattering data below a given maximal energy to provide constraints on the interaction transforms the whole setup into a statistical inference problem, based on the conventional least χ 2 -method. This fact was recognized already in 1957 [1] (see Ref. [2] for an early review) and, after many years, culminated in the admirable Nijmegen partial wave analysis (PWA) of 1993 [3], based on the crucial observations that charge-dependent one-pion-exchange (CD-OPE), tiny but essential electromagnetic and relativistic effects, and a judicious selection of the scattering database could actually provide a satisfactory fit with χ 2 /datum ∼ 1 for a total number of data consisting, as of 1993, of 1787 pp and 2514 np (normalizations included) at the 3 σ level. These criteria have set the standard for PWA's and the design of high quality phenomenological potentials [4][5][6][7][8][9][10][11][12]. The inference point of view is mainly phenomenological and requires a balanced interplay between which data qualify as constraints and which models provide the most likely description of the data. None of these choices is free of prejudices and they are actually intertwined; a circumstance that should be kept in mind when assessing the reliability and predictive power of the theory aiming at a faithful representation of the input data and their uncertainties.
The quantum mechanical nature of the PWA with a given cutoff in energy leads to inverse scattering ambiguities which increase at short distances (see, for example, Refs. [13,14] and references therein). Remarkably, a universal and model-independent low-energy interaction arises when unobserved high energy components above the cutoff are explicitly integrated out of the Hilbert space preserving the scattering amplitude [15,16]. While this V low−k framework is an extremely appealing setup based on Wilsonian renormalization, to date this universal interaction has not been determined from data directly and one has to proceed via a fitted and bare N N interaction since off-shellness is required [17]. However, inferring a N N interaction from data, is not the full story, and three-nucleon, and possibly higher multi-nucleon, interactions are needed to describe residual contributions to nuclear binding energies [18]. As is well known, their strength and form are also affected by the chosen off-shell behavior of the N N interaction and a universal V low−k three-nucleon interaction remains to be found.
In an ideal situation all steps in the inference process, including the scattering data selection itself, should be carried out with the "true" theory, which for nuclear physics is quantum chromodynamics (QCD), the fundamental theory of interacting quarks and gluons. Assuming, as we do, that the theory is correct, QCD would just tell us which experiments are right and which are wrong, or whether the reported uncertainties are realistic with a given confidence level on the side of the experiment. At the same time one would set constraints on the QCD parameters such as the light quark masses and Λ QCD , or equivalently the pion mass m π and the pion weak decay constant F π . While there has been impressive progress in bringing lattice QCD simulations for light quarks closer to nuclear physics working conditions (see Refs. [19,20] and references therein), we do not yet envisage, at least not in the near future, the realization of conditions that would allow one to establish, on QCD grounds, the correctness of the about 8000 currently available published pp and np scattering data below pion production threshold. Instead, already in the early 90's the phenomenological analysis carried out by the Nijmegen group made it possible to pin down the pion masses with a precision of 1 MeV from their PWA of pp and np data [21].
In practice, we must content ourselves with an approximation scheme to the true theory in conjunction with a phenomenological approach. This specifically means assuming a sufficiently flexible parametrization of the interaction in terms of the relevant degrees of freedom which does not overlook some relevant physical feature. In what follows it is instructive to briefly review both the process and criteria taken into account to select a consistent database as well as the QCD-based theory used to describe it. Our aim is to make the reader aware of all the fine details which are needed in order to credibly falsify the theoretical model, QCD grounded or not, against the data and keep an open mind about the out-coming result.
On the theoretical side, we will assume along with Weinberg [22] that there is a chiral effective field theory (χEFT) capable of systematically describing the strong interactions among nucleons, ∆-isobars, and pions, as well as the electroweak interactions of these hadrons with external (electroweak) fields. In the specific case of two nucleons, the requirements imposed by χEFT can be incorporated into a non-relativistic quantum mechanical potential, constructed by a perturbative matching, order by order in the chiral expansion, between the on-shell scattering amplitude and the solution of the Schrödinger equation (see, for example, the review paper by Machleidt and Entem [23]). Such a theory provides the most general scheme accommodating all possible interactions compatible with the relevant symmetries of QCD at low energies, in particular chiral symmetry. By its own nature, χEFT needs to be organized within a given power counting scheme and the resulting chiral potentials can conveniently be separated into long-and short-distance contributions, the latter (short-distance ones) featuring the needed counter-terms for renormalization. At leading order in the chiral expansion one has the venerable one-pion-exchange (OPE) potential which, as already mentioned, emerges as a universal and indispensable long-distance feature for an accurate description of proton-proton and neutron-proton scattering data [3]. Higher orders in the chiral expansion incorporate the two-pion-exchange (TPE) potential [24], due to leading and sub-leading πN couplings (the sub-leading couplings c 1 , c 3 , and c 4 can consistently be obtained from low energy πN scattering data). The inclusion of TPE allows one to reduce the short-range cutoff separating long-and short-distance contributions, which helps in reducing the impact of details in the unknown shortdistance behavior of the potentials. Nonetheless, we will note in Sec. IV that uncertainties are dominated by this diffuse separation between short and long distances.
There are many practical advantages deriving from a χEFT that explicitly includes ∆-isobar degrees of freedom, the most immediate one being a numerical consistency between the values of the low-energy constants c 1 , c 3 and c 4 inferred from either πN or N N scattering. Such a theory also naturally leads to three-nucleon forces induced by TPE with excitation of an intermediate ∆ (the Fujita-Miyazawa three-nucleon force) as well as to two-nucleon electroweak currents (see for example Ref. [25]). In addition, there are rather strong indications from phenomenology that ∆ isobars play an important role in nuclear structure and reactions. An illustration of this are the three-nucleon forces involving excitation of intermediate ∆'s, needed to reproduce the observed energy spectra and level ordering of lowlying states in s-and p-shell nuclei or the correct spin-orbit splitting of P-wave resonances in low-energy n-α scattering (for a review, see Ref. [18]). Another illustration is the relevance of electroweak N -to-∆ transition currents in radiative and weak capture processes involving few-nucleon systems [26], specifically the radiative captures of thermal neutrons on deuteron and 3 He [27,28] or the weak capture of protons on 3 He (the so-called hep process) [29]. It is for these reasons that in the present work we construct a minimally non-local coordinate space chiral potential, that includes ∆ intermediate states in its TPE component-it is described in detail in Sec. II. Such a coordinate-space representation offers many computational advantages for ab initio calculations of nuclear structure and reactions, in particular for the type of quantum Monte Carlo calculations of s-and p-shell nuclei very recently reviewed in Ref. [18].
On the experimental side, there are currently ∼ 8000 published pp and np scattering data below pion production threshold corresponding to 24 different scattering observables, including differential cross sections, spin asymmetries, and total cross sections [30,31], see Ref. [12] for updated pp and np abundance plots in the (E lab , θ cm ) plane. However, not all of these data are mutually compatible and a decision has to be made as to which are more likely to be correct. In principle, the N N scattering amplitude can be determined uniquely, provided a complete set of experiments is given-a rare situation for the case under consideration. Therefore, a theoretical model is needed to provide a smooth energy dependence which allows one to interpolate between different energy values, and helps in deciding on the mutual consistency of nearby data in (E lab , θ cm ) plane. The PWA carried out in Granada parametrizes [10] 1 the interaction, for inter-nucleon distances r less than 3 fm, in terms of a set equidistant delta-shells separated by ∆r = 0.6 fm (in other words, a coarse-grained parametrization), while retaining only the OPE component for r > 3 fm. The choice of ∆r corresponds to the shortest de Broglie wavelength at about pion production threshold, and consequently all the data are weighted with their quoted experimental uncertainty. The result of the analysis has been a 3 σ self-consistent database comprising a total of 6713 pp and np scattering data. More details on the data analysis specific to our potential are presented in Sec. III. One important aspect of the Granada PWA is the correlation pattern among the fitting parameters, namely different partial waves are mostly uncorrelated which, together with the large number of selected data, speaks in favor of a lack of bias in the selection process. Actually the correlation length which decides on the specific form of the potential should be smaller than the distance ∆r = 0.6 fm in the coarse-grained parametrization.
Chiral potentials have been subjected to PWA and confronted to pp and np scattering data up to lab energy of 350 MeV. Within the χEFT framework the Nijmegen group used the TPE potential [24] to carry out pp [6] and np + pp [8] analyses determining the chiral constants c 3 and c 4 from these data while constraining c 1 from πN data. Taking the chiral constants from πN analyses, Entem and Machleidt [32] used a next-to-next-to-next-to-leading order (N3LO or Q 4 , Q generically specifying the low momentum scale) chiral potential to fit pp and np scattering data up to lab energy of 290 MeV. The resulting χ 2 /datum were 1.1 for 2402 np data and 1.50 for 2057 pp data, and consequently a global χ 2 /datum of 1.28. The chiral TPE potential [24] was also used within the coarse grained framework to determine the chiral constants in Ref. [11] with a global χ 2 /datum of 1.07, based on 6713 pp and np scattering data.
Other available chiral potentials [33,34] have not been confronted to scattering data directly but rather to phase shifts obtained in the Nijmegen analysis (the recent upgrade [35] of Ref. [33] relies on the same procedure, while in Ref. [34] a study of peripheral phase shifts is carried out with two-and three-pion exchange potentials up to order Q 5 ). As we will show in Sec. IV, there is a substantial difference between fitting scattering data and fitting phase shifts mainly because of the existing correlations among the many partial waves and mixing angles. Actually, a good χ 2 -fit to phase shifts may yield quite a bad χ 2 in a fit to data. Moreover, the spread in phase-shift values among different high-quality potentials fitting the same data reflects the differences in the potential representation and turns out to be larger than the estimated statistical errors (compare Fig. 1 of Ref. [36] with Fig. 3 of Ref. [37]). The consequences of these larger errors have been discussed in Ref. [38].
The previous comments address the use of chiral potentials to fit selected N N scattering databases which have been obtained from phenomenological representations of the interactions. An obvious question which comes to mind is whether chiral potentials, being credible and general low energy representations of QCD in the N N sector, should be used themselves to select the database. Within the coarse grained framework the impact of chiral interactions on the selection of the database has also been studied in Ref. [11]. The result was that a larger number of data were rejected but at the same time the number of parameters was reduced. This poses the interesting question on what is the meaning of improvement-a particularly critical issue when the potential itself (chiral or not) must be tested against the selected data. Obviously an incorrect model will appear to be correct if a sufficiently large number of data is discarded. However, the theory with just delta-shells+OPE is more general than that with delta-shells+(OPE+TPE), and hence data selection based on the former is more reliable. In any case, the results of Ref. [11] show also that the long range part of the next-to-next-to-leading order (N2LO or Q 3 ) chiral potential can indeed fit the delta-shells+OPE selected data satisfactorily with a χ 2 /datum of 1.07, when the potential is taken to be valid for inter-nucleon distances ranging from 1.8 fm outwards.
The present paper is organized as follows. In the next section we describe the potential, while in Sec. III we provide a brief discussion of the data fitting. In Sec. IV we report the χ 2 values obtained in the fits as well as the values for the low-energy constants that characterize the potential, and show the calculated phase shifts for the lower partial waves (S, P, and D waves) and compare them to those from recent PWA's. There, we also provide tables of the pp, np and nn effective range parameters and of deuteron properties, including a figure of the deuteron S and D waves. Finally, in Sec. V we summarize our conclusions. A number of details are relegated to Appendices A-E.

II. POTENTIALS
The two-nucleon potential includes a strong interaction component derived from χEFT up to next-to-next-to-nextto-leading order (N3LO or Q 4 ) and denoted as v 12 , and an electromagnetic interaction component, including up to terms quadratic in the fine structure constant α (first and second order Coulomb, Darwin-Foldy, vacuum polarization, and magnetic moment interactions), and denoted as v EM 12 . The v EM 12 component is the same as that adopted in the Argonne v 18 (AV18) potential [5]. The component induced by the strong interaction is separated into long-and short-range parts, labeled, respectively, v L 12 and v S 12 . The v L 12 part includes the one pion-exchange (OPE) and two pion-exchange (TPE) contributions, illustrated in Fig. 1: panel (a) represents the OPE contribution at leading order (LO or Q 0 ); panels (b)-(g) represent the TPE contributions at next-to leading order (NLO or Q 2 ) without and with ∆-isobars in the intermediate states; lastly, panels (h)-(p) represent sub-leading TPE contributions at next-to-next-to leading order (N2LO or Q 3 ). The NLO and N2LO loop corrections contain ultraviolet divergencies, which are isolated in dimensional regularization and then reabsorbed into contact interactions by renormalization of the associated low energy constants (LEC's) [39,40]. Additional loop corrections at NLO and N2LO only lead to renormalization of OPE and contact interactions [39,41], and will not be discussed any further here.
πN [46] and L (2) πN ∆ [43]. Note that relativistic 1/MN -corrections (MN is the nucleon mass) included in L πN Lagrangian are not considered here. In particular the contributions of diagrams (i), (k) and (n) are neglected. The LO, NLO, and N2LO terms are well known, and explicit expressions for them can be found in Refs. [39,40,[42][43][44]. The LO and NLO terms depend on the the pion decay amplitude F π , and the nucleon and N -to-∆ axial coupling constants, respectively g A and h A = 3 g A / √ 2 (this value for h A is from the large N c expansion or strong-coupling model [45], and is in good agreement with the value inferred from the empirical ∆-width). The sub-leading N2LO terms also depend on the LEC's c 1 , c 2 , c 3 , and c 4 and the combination of LEC's (b 3 + b 8 ), respectively from the second order πN and πN ∆ chiral Lagrangians L (2) πN [46] and L (2) πN ∆ [43]. The values of these LEC's, as determined by fits to πN scattering data [43], and of the masses and other physical constants adopted in the present study are listed in Tables I and II. In the static limit, the momentum-space LO, NLO, N2LO terms are functions of the momentum transfer k; hereafter, we define k = p −p and K = (p +p)/2, where p and p are the initial and final relative momenta of the two nucleons.
Coordinate-space expressions for the TPE terms are obtained by using the spectral function representation [44], however with no spectral cutoff 2 , v l,TPE with O l=1,...,6 denoted as c, τ, σ, στ, t, tτ . Those corresponding to diagrams (b)-(d) and (h)-(k) in Fig. 1 are known in closed form (see, for example, Ref. [44]) and are listed in Appendix A for completeness; the remaining ones corresponding to diagrams (e)-(g) and (l)-(p) have been derived in terms of a parametric integral, and they too are given in Appendix A. The radial functions v l L (r) are singular at the origin (they behave as 1/r n with n taking on values up to n = 6, see Refs. [47,48] for analytical expressions), and each is regularized by a cutoff of the form where in the present work three values for the radius R L are considered R L = (0.8, 1.0, 1.2) fm with the diffuseness a L fixed at a L = R L /2 in each case. The potential v L 12 , including the well known OPE components at LO regularized by the cutoff in Eq. (2.3), then reads in coordinate space

5)
O σT 12 = σ 1 · σ 2 T 12 , and O tT 12 = S 12 T 12 , and T 12 = 3 τ 1z τ 2z − τ 1 · τ 2 is the isotensor operator. The terms proportional to T 12 account for the charge-independence breaking induced by the difference between the neutral and charged pion masses in the OPE. However, this difference is ignored in the NLO and N2LO loop corrections which have been evaluated with m π = (2 m π + + m π 0 ) /3. Additional (and small) isospin symmetry breaking terms arising from OPE [49] and TPE [50] and from OPE and one-photon exchange [51,52] have also been neglected.
The potential v S 12 includes charge-independent (CI) contact interactions at LO, NLO and N3LO, and chargedependent (CD) ones at LO and NLO, in momentum- where S 12 (k) = 3 σ 1 · k σ 2 · k − k 2 σ 1 · σ 2 , C S and C T are the LO LEC's in standard notation, while C i=1,...,7 and D i=1,...,15 are generally linear combinations of those in the "standard" set, as defined, for example, in Ref. [23]. In the NLO and N3LO contact interactions terms proportional to K 2 and K 4 , which would lead to p 2 and p 4 operators in coordinate space (p −→ −i∇ is the relative momentum operator), have been removed by a Fierz rearrangement, for example with m = 2 or 4. Of course, mixed terms of the type k 2 K 2 or K × k cannot be Fierz-transformed away. In the potential v S,CD 12 (k, K) only terms up to NLO, involving charge-independence breaking (proportional to T 12 ) and charge-symmetry breaking (proportional to τ 1z + τ 2z ), are accounted for. The associated LEC's, while providing some additional flexibility in the data fitting discussed below (especially C IV 0 in reproducing the singlet nn scattering length), are not well constrained.
A couple of comments are now in order. The first is that strict adherence to power counting would require inclusion of additional one-loop as well as two-loop TPE and three-pion exchange contributions at order Q 4 . These contributions have been neglected, since they are known to be small (see, for example, Ref. [23]). Furthermore it is the D i LEC's at Q 4 that are critical for a good reproduction of phase shifts in lower partial waves, particularly D-waves, and a good fit to the N N database [23] in the 0-300 MeV range of energies considered in the present study.
The second comment is in reference to isospin symmetry breaking. We have not included explicitly contributions from OPE and one-photon exchange [51,52]. As noted in Ref. [11], this π-γ interaction is small and ambiguous, and requires regularization at short distances. So its main effect can be effectively shifted into a counter-term. While this can be improved, we will see below our final fitting results do not seem to require these long-range isospin breaking effects.
The potential v S 12 (k, K) is regularized via a Gaussian cutoff depending only on the momentum transfer k, which leads to a coordinate-space representation only mildly non-local, containing at most terms quadratic in the relative momentum operator. It reads (see where O l=1,...,6

12
have been defined above, referred to as T , τ z, σT , στ z, tT , tτ z, bT , bτ z. The four additional terms, denoted as p, pσ, pt, and ptτ , in the anti-commutator of Eq. (2.10) are p 2 -dependent. We consider, in combination with R L = (0.8, 1.0, 1.2) fm, R s = (0.6, 0.7, 0.8) fm, corresponding to typical momentum-space cutoffs Λ S = 2/R S from about 660 MeV down to 500 MeV. While the use of a Gaussian cutoff mixes up orders in the power counting-for example, the LO contact interactions proportional to C S and C T in Eq. (2.6) generate contributions at NLO and N3LO-such a choice nevertheless leads to smooth functions for the potential components v l S (r) and the resulting deuteron waves. Sharper cutoffs, like those ∝ exp [−(r/R) n ] with n = 4, as suggested in Ref [53], or n = 6, as in one of the earlier versions of the present model, generate wiggles in the deuteron waves at r ∼ R (as well as mixing of power-counting orders).

III. DATA ANALYSIS
Setting aside electromagnetic (EM) contributions (Coulomb and higher order ones) for the time being, the invariant on-shell scattering amplitude M for the N N system can be expressed in terms of five independent complex functionsthe Wolfenstein parametrization-as wherel,m,n are three orthonormal vectors along the directions of p + p, p − p, and p × p , and p , p are the final and initial relative momenta, respectively. The functions a, m, g, h, and c are taken to depend on the energy in the laboratory (lab) frame and the scattering angle θ in the center-of-mass (cm) frame. Any scattering observable can be constructed out of these amplitudes [30,31]. The N N amplitude is diagonal in pair spin S, and pair isospin and isospin projection T M T , and is expanded in partial waves as where L and J denote respectively the orbital and total angular momenta, the . . . , the S-matrix is simply given by in single channels with L = L = J, and by in coupled channels with S = 1 and L, L = J ∓ 1 ( J is the mixing angle). Hereafter, for notational simplicity we drop from the phase shifts unnecessary subscripts as well as the superscripts T M T , with T = 1 and M T = 1, 0, −1 for respectively pp, np, and nn. The S-matrix elements and phase shifts are obtained from solutions of the Schrödinger equation with suitable boundary conditions, as discussed Appendix C. In terms of the amplitudes M S M S M S , the functions a, m, g, h, and c then read 9) and this can be further simplified by noting that When EM interactions are included, the full scattering amplitudes M are conveniently separated into a part due to nuclear interactions and another one stemming from EM interactions, (3.10) The pp EM amplitudes contain Coulomb with leading relativistic corrections, vacuum polarization, and magnetic moments contributions, whereas the np ones contain magnetic moment contributions only (see Ref. [10] for a compendium of formulas and references to the original papers; for completeness, however, the determination of the pp phase shifts relative to EM functions and of the pp effective range expansion is summarized in Appendix D). Due to the finite range of the N N force, the nuclear part of the scattering amplitudes, M N , converges with a maximum total angular momentum of J = 15. In contrast, EM scattering amplitudes, M EM , require a summation of about thousand partial waves due to the long range and tensor character of the dipolar magnetic interactions. While these corrections are numerically tiny, they are nevertheless indispensable for an accurate description of the data [54]. We use the database developed in Granada and specified in detail in Ref. [10], where a selection of the large collection of np and pp scattering data taken from 1950 till 2013 was made. The adopted criterium was to represent the N N interaction with a general and flexible parametrization, based on a minimal set of theoretical assumptions so as to avoid any systematic bias in the selection process. The aim of the method, first suggested by Gross and Stadler [9], was to obtain a 3 σ self-consistent database. This entails removing 3 σ outliers and re-fitting iteratively until convergence. The procedure results in a database with important statistical features [12] and therefore amenable to statistical analysis, and leads to the identification of a consistent subset among the large body of 6713 np and pp experimental cross sections and polarization observables 3 . In the present study, in particular, we are concerned with a subset of this 3 σ-self-consistent database, namely data below 300 MeV lab energy. This database is organized in the following way: there are N sets of data, each one corresponding to a different experiment. Each data set contains measurements at fixed E lab and different scattering angles θ. However a few observables are measured at different E lab and fixed θ, like, for example, total cross sections since their measurement does not involve the scattering angle (θ = 0). An experiment may have a specified systematic error (normalized data), no systematic error (absolute data), or an arbitrarily large systematic error (floated data).
We briefly describe the fitting procedure. The total figure of merit is defined as the usual χ 2 function where χ 2 t refers to the corresponding contribution from each data set, which we explain next. In all cases, the χ 2 t for a data set is given by where o i and t i are the measured and calculated values of the observable at point i, δo i and δ sys are the statistical and systematic errors, respectively, and Z t is a scaling factor chosen to minimize the χ 2 t , (3.13) The last term in Eq. (3.12) is denoted χ 2 sys . For absolute data Z = 1 and χ 2 sys = 0, while for floated data use of Eq. (3.13) is made with δ sys = ∞ so that χ 2 sys = 0. Normalized data have in most cases Z = 1 such that χ 2 sys = 1 and the normalization is counted as an extra data point 4 . For some normalized data the systematic error can give a rather large χ 2 sys due to an underestimation of δ sys . In order to account for this, we float data that have χ 2 sys > 9 and no extra normalization data is counted. This is in line with the criterion used to build the pp and np database. Finally, the total χ 2 is the sum of all the χ 2 t for each pp and np data set. The minimization of the objective function χ 2 with respect to the LEC's in Eqs. (2.6) and (2.7) is carried out with the Practical Optimization Using no Derivatives (for Squares), POUNDerS [55]. This derivative-free algorithm is designed for minimizing sums of squares and uses interpolation techniques to construct residuals at each point. In the optimization procedure, we fit first phase shifts and then refine the fit by minimizing the χ 2 obtained from a direct comparison with the database. In fact, sizable changes in the total χ 2 are found when passing from phase shifts to observables, so this refining is absolutely necessary to claim reasonable fits to data. This is a general feature which is often found, and reflects the different weights in the χ 2 contributions of the two different fitting schemes. Indeed, the initial guiding fit to phase shifts chooses a prescribed energy grid arbitrarily, which does not correspond directly to measured energies, nor necessarily samples faithfully the original information provided by the experimental data. Moreover, there are different PWA's which describe the same data but yield different phase shifts with significantly larger discrepancies than reflected by the inferred statistical uncertainties [10][11][12].

IV. RESULTS
We report results for the potentials v 12 + v   Granada database of pp and np cross sections, polarization observables, and normalizations up to lab energies of 300 MeV, to the pp, np, and nn singlet scattering lengths, and to the deuteron binding energy. We list the number of pp and np data (including normalizations) and corresponding total χ 2 for the three models in Table III, where we also report for comparison the χ 2 corresponding to the AV18 [5] (of course, without a refit of it) and the same database.  Table IV. The values for the πN LEC's in the OPE and TPE terms of these models have already been given in Tables I and II. It is interesting to examine the extent to which these LEC's satisfy the requirement of naturalness. To this end, following Machleidt and Entem [23], we note that this criterion would imply that the LEC's of the charge-independent part v S,CI 12 of the contact potential have the following magnitudes where f π = 92.4 MeV and Λ χ = 1 GeV. A glance at Table IV indicates that the LEC's are generally natural, but for the following exceptions: C S,T in model c, C 7 in models b and c, and D 1 , D 10 , and D 12 in al three models considered. As already noted, however, the use of a (momentum-space) Gaussian cutoff mixes orders in the power expansion, since and, as an example, the spin-isospin independent central component of v S,CI 12 , after inclusion of this cutoff, is modified as suggesting that some of the LEC's multiplying terms linear and quadratic in k 2 may not be natural after all. In order to estimate the size of the (nominally) LO (Q 0 ) and NLO (Q 2 ) LEC's associated with the charge-dependent part v S,CD 12 of the contact potential, we note that the terms proportional to C IV 0 and C IT 0 in Eq. (2.7) should scale respectively as m 2 π and ∆m 2 π , where is related to the u-d quark mass difference-we assume that ∼ e = √ 4πα, e being the electric charge and α the fine structure constant-and ∆m 2 π is the squared-mass difference between the charged and neutral pions. Consequently, one would expect for the LO LEC's and for the NLO LEC's These expectations are not borne out by the actual values reported in Table IV. Particularly striking are the very large values obtained for the LEC's C IV 4 and C IT 4 associated with the spin-orbit term.  The S-wave, P-wave, and D-wave phase shits for np (in T = 0 and T = 1) and pp are displayed in Figs. 2-4 up to 300 MeV lab energies. The phases calculated with the full models a, b, and c including strong and electromagnetic interactions are represented by the band. The np phases are relative to spherical Bessel functions, while the pp phases are with respect to electromagnetic functions (see Appendix D). The cutoff sensitivity, as represented by the width of the shaded band, is very weak for pp, and generally remains modest for np, except for the T = 0 3 D 3 phase and 1 mixing angle, particularly for energies larger than 150 MeV. The calculated phases are compared to those obtained in partial-wave analyses (PWA's) by the Nijmegen [3,4], Granada [10], and Gross-Stadler [9] groups. Note that the recent Gross and Stadler's PWA was limited to np data only. We also should point out that, since the Nijmegen's PWA of the early nineties which was based on about 1780 pp and 2514 np data in the lab energy range 0-350 MeV, the N N elastic scattering database has increased very significantly. Indeed, in the same energy range the 2013 Granada database contains a total of 2972 pp and 4737 np data. Especially for the higher partial waves in the np sector and at the larger energies there are appreciable differences between these various PWA's. It is also interesting to observe that these differences are most significant for the T = 0 3 D 3 phase and 1 mixing angle, and therefore correlate with the cutoff sensitivity displayed in these cases by models a, b, and c.
The low-energy scattering parameters are listed in Table V, where they are compared to experimental results. The singlet and triplet np, and singlet pp and nn, scattering lengths are calculated with and without the inclusion of electromagnetic interactions. Without the latter, the effective range function is simply given by F (k 2 ) = k cot δ = −1/a + r k 2 /2 up to terms linear in k 2 . In the presence of electromagnetic interactions, a more which imply that charge symmetry and charge independence are broken respectively by and ∆a CIB = (a N pp + a N nn )/2 − a N np = 5.6 ± 0.6 fm , ∆r CIB = (r N pp + r N nn )/2 − r N np = 0.03 ± 0.13 fm . (4.10) The more significant values for ∆a CSB and ∆a CIB can be compared to those inferred from  Fig. 5 we show the 1 S 0 phase shifts for pp, np and nn calculated with and without the inclusion of electromagnetic interactions (only model b is considered). There is excellent agreement between these phases and those obtained in the the Granada, Gross and Stadler, and Nijmegen PWA's, when electromagnetic effects are fully accounted for. Particularly at low energies (see Fig. 6), the latter provide most of the splitting between the pp and np phases, with remaining differences originating from isospin symmetry breaking due to the OPE term in v L 12 and the central terms in v S,CD 12 , proportional to the LEC's C IT i and C IV i with i = 0-2. In the absence of electromagnetic interactions, the splitting between the pp and nn 1 S 0 phases is induced by the charge-symmetry breaking terms of v S,CD 12 proportional to the LEC's C IV i with i = 0-2; it is smaller than that between pp and np 1 S 0 phases. The effects of isospin symmetry breaking are also seen in the pp and np 3 P J phases with J = 0, 1, 2 in the upper right and lower panels of Fig. 5, especially at the higher energies. The calculated phases, which correspond again to model b, include electromagnetic effects, but the latter are negligible beyond 100 MeV. The splitting between the pp and np 3 P J phases is mostly due to the isotensor and isovector terms of v S,CD The static deuteron properties are shown in Table VI and compared to experimental values [60][61][62][63][64]. The binding energy E d is fitted exactly and includes the contributions (about 20 keV) of electromagnetic interactions, among which the largest is that due to the magnetic moment term. The asymptotic S-state normalization, A S , and the D/S ratio, η, are both ∼ 2 standard deviations from experiment for all models considered. The deuteron (matter) radius, r d , is exactly reproduced with model b, but is under-predicted (over-predicted) by about 1.4% (0.7%) with model a (model c). It is should be noted that this observable has negligible contributions due to two-body electromagnetic operators [65]. The magnetic moment, µ d , and quadrupole moment, Q d , experimental values are underestimated by all three models, but these observables are known to have significant corrections from (isoscalar) two-body terms in  nuclear electromagnetic charge and current [65]. Their inclusion would bring the calculated values considerably closer to, if not in agreement with, experiment. Finally, the S-and D-wave components of the deuteron wave function are displayed in Fig. 7, where they are compared to those of the Argonne v 18 (AV18) model. There is significant cutoff dependence as (R L , R S ) are reduced from the values (1.2, 0.8) fm of model a down to (0.8, 0.6) fm of model c. For r 1 fm, the S-wave becomes smaller (is pushed out), while the D-wave becomes larger (is pushed in) in going from model a to model c. The D-state percentage increases correspondingly (see Table VI).
We note in closing that in Appendix E we provide figures of the various components of potential models a, b, and c (their charge-independent parts only) as well as tables of numerical values for the pp and np S, P, D, F, and G phase shifts obtained with model b.

V. CONCLUSIONS
In the present study, we have constructed a coordinate-space nucleon-nucleon potential with an electromagnetic interaction component including first and second order Coulomb, Darwin-Foldy, vacuum polarization, and magnetic moment terms, and a strong interaction component characterized by long-and short-range parts. The long-range part includes OPE and TPE terms up to N2LO, derived in the static limit from leading and sub-leading πN and πN ∆ chiral Lagrangians. Its strength is fully determined by the nucleon and nucleon-to-∆ axial coupling constants g A and h A , the pion decay amplitude F π , and the sub-leading LEC's c 1 , c 2 , c 3 , c 4 , and b 3 + b 8 , constrained by reproducing πN scattering data (the values adopted for all these couplings are listed in Table I). In coordinate space, this long-range part is represented by charge-independent central, spin, and tensor components without and with the isospin dependence τ 1 · τ 2 (the so-called v 6 operator structure), and by charge-dependence-breaking central and tensor components induced by OPE and proportional to the isotensor operator T 12 . The short-range part is described by charge-independent contact interactions specified by a total of 24 LEC's (2 at LO, 7 at NLO, and 15 at N3LO) and by charge-dependent ones characterized by 10 LEC's (2 at LO and 8 at NLO), 5 of which multiply charge-symmetry breaking terms proportional to τ 1z + τ 2z and the remaining 5 multiply chargedependence breaking terms proportional to T 12 . In the NLO and N3LO contact interactions, Fierz transformations have been used in order to rearrange terms that in coordinate space would otherwise lead to powers of p-the relative momentum operator-higher than two. The resulting charge-independent (coordinate-space) potential contains, in addition to the v 6 operator structure, spin-orbit, L 2 , quadratic-spin-orbit, and p 2 components, while the chargedependent one retains central, tensor, and spin-orbit components.
The 34 LEC's in the short-range potential have been constrained by fitting 5291 pp and np scattering data (including normalizations) up to 300 MeV lab energies, as assembled in the Granada database, and the pp, np, and nn scattering lengths, and the deuteron binding energy. The global χ 2 (pp + np)/datum is 1.33 for the three different models we have investigated, each specified by a pair of (coordinate-space) cutoffs, respectively, R L and R S for the long-and short-range parts: (R L , R S ) = (1.2, 0.8) fm for model a, (1.0, 0.7) fm for model b, and (0.8, 0.6) fm for model c. These cutoffs are close to the 1/(2 m π ) ∼ 0.7 fm TPE range. The values of the LEC's corresponding to the three models are given in Table IV. They are generally of natural size, but for a few exceptions, most notably the LEC's C IV 4 and C IT 4 multiplying the charge-dependent spin-orbit terms, which lead to relatively large splitting between the pp and np 3 P 0 and 3 P 1 phase shifts-a splitting that is not consistent with that obtained in both the Nijmegen and Granada PWA's. It should also be noted that the degree of unnaturalness increases as the short-distance cutoffs are reduced.
Our results suggest that discrepancies between the phases calculated here and those from available PWA's in some of the partial waves, such as the 1 mixing angle, could hardly be resolved by carrying out the database selection using the present interaction. We should also note that the renowned Entem and Machleidt N3LO fit up to E lab = 290 MeV provides a χ 2 /datum of 1.1 for 2402 np data and 1.5 for 2057 pp, and hence a global χ 2 /datum of 1.3. In our case, we describe 2161 (2764) scattering data and 148 (218) normalizations for pp (np), which means that the average contribution to the χ 2 from each additional datum is homogeneous and of order one out of about 800 extra data. So, our fit is as good as the one of Entem and Machleidt with these additional data.
According to our findings the largest uncertainty in the chiral theory when fitting up to a maximum lab energy of 300 MeV is provided by the cutoff dependence. Under these circumstances it makes little sense to analyze further uncertainties, but it is nonetheless surprising that precisely the model implementing many QCD motivated theoretical constraints should end up magnifying the uncertainty to a larger extent than the spread historically found in all so far successful PWA's to pp and np scattering data. On the other hand, the reliability of the long distance chiral interaction does not depend on how the short distance unknown interaction is organized. This has been proven by the first chiral potential fits by the Nijmegen group from their pp [6] and np + pp [8] analyses and more recently verified with increased statistics by the Granada group [11]. This leaves open the possibility that better fits than those found here should be possible by properly altering the short distance structure. This point has recently been discussed in Ref. [67].
Of course, this cutoff uncertainty could be greatly reduced if the fitting energy range were to be lowered so as to ensure that differences between fitted data and fitting theory fulfill the normality requirement and, at the same time, statistical uncertainties remain at the same level as cutoff uncertainties. Following the recent suggestion [67], we find that this happens with the current form of the potential when E lab ≤ 125 MeV. In a companion paper we will analyze the statistical properties of the present fit and how there is a trade-off of different uncertainty sources.
We conclude by observing that, apart from the p 2 -dependent terms, the potential constructed here has the same operator structure of the AV18, and is of slightly better quality than the AV18 (the AV18 global χ 2 (pp + np)/datum on the same database up to 300 MeV lab energies is 1.46). It should be fairly straightforward to incorporate it in the few-nucleon calculations based on hyperspherical-harmonics expansion techniques favored by the Pisa group [66], or in the quantum Monte Carlo ones preferred by the ANL/ASU/JLab/LANL collaboration [18]. The Fortran computer program generating the potential will be made available upon request. The LO (OPE) terms corresponding to diagram (a) in Fig. 1 and x α = m πα r. The NLO terms corresponding to diagrams (b)-(d) read [24] v NLO where x = m π r (m π is the average pion mass) and K n are modified Bessel functions of the second kind. The NLO terms corresponding to diagrams (e)-(f) with a single ∆ intermediate state are given by where y = ∆M r (∆M is the ∆-nucleon mass difference) and the parametric integral over µ is carried out numerically. The NLO terms corresponding to diagram (g) with 2 ∆ intermediate states are Moving on to the loop corrections at N2LO, the terms corresponding to diagrams (h)-(k) are given by v N2LO c (r; / ∆) = 3 2 π 2 r 6 while those corresponding to diagrams (l)-(o) are given by Lastly, the contributions corresponding to diagram (p) read v N2LO c (r; 2∆) = − 2 81π 3 r 6 The radial functions of the charge-independent part of the potential v L 12 in Eq.
where r is the relative position and K −→ p = −i ∇ δ(r − r), the relative momentum operator. For the momentumspace operator structures present in Eqs. (2.6) and (2.7) one finds: Using the above expressions, the functions v l S (r) are obtained as v σ S (r) = C T C RS (r) + C 3 −C v bb S (r) = −D 9 1 v pσ S (r) = D 13 −C v pt S (r) = −D 14 C Note that in Eqs. (B8) and (B9) only the terms proportional to L 2 and (L · S) 2 are retained.
with the boundary conditions (again, reinstating superscripts and subscripts) w T SJ L L (r) r where L = J ∓ 1 is the orbital angular momentum of the incoming wave.
Appendix D: pp phase shifts and effective range expansion We discuss briefly the calculation of the pp phase shifts and effective range expansion with inclusion of the full electromagnetic potential v EM 12 [5]. Radial wave functions behave in the asymptotic region (r 30 fm) as where L = J for single channels or L = L = J ∓ 1 for coupled channels (the pair isospin and spin subscripts T and S have been dropped for simplicity), h where V C1 (V C2 ) and V V P are respectively the first-order (second-order) Coulomb and vacuum polarization terms. These terms include form factors to remove singularities in the r = 0 limit [5]. Note that the Darwin-Foldy and magnetic moment corrections are not included above, since at large r the former falls off exponentially and the latter behaves as 1/r 3 . Following Ref. [69] and treating the V C2 (r) and V V P (r) corrections in first order perturbation theory, one finds that F L (kr; η and G L (kr; η ) can be expressed as where the F L and G L are standard Coulomb functions, the function V (r) is proportional to V C2 (r) and V V P (r), and the phase shifts ρ L and τ L corresponding, respectively, to V C2 and V V P are given (in first order perturbation theory) by tan(ρ L + τ L ) ρ L + τ L = − ∞ 0 dr F L (kr; η ) V (r) F L (kr; η ) .
In the absence of V C2 and V V P , the solutions F L and G L reduce to the regular and irregular Coulomb functions. In the computer programs Eqs. (D5)-(D6) are used to construct the EM functions and Eq. (D8) to obtain the phase shifts ρ L and τ L . The effective range expansion in the 1 S 0 channel is obtained as [68][69][70] F EM (k 2 ) = − 1 a EM + 1 2 r EM k 2 + . . . , with F EM (k 2 ) = k C 2 0 (η ) where γ is Euler's constant, and the function I(r) entering the vacuum polarization potential V V P (r) is defined as in Ref. [69], The effective range function F EM (k 2 ) corresponding to model b is shown in Fig. 8. The numerical methods are stable down to lab energies of 1 keV.