{'fraction_non_alphanumeric': 0.03146252285191956, 'fraction_numerical': 0.01393053016453382, 'mean_word_length': 4.419159896968496, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Low-energy theorems for nucleon-nucleon scattering at unphysical pion masses

The longest-range part of the nuclear force due to the one-pion exchange governs the energy dependence of the scattering amplitude in the near-threshold region and imposes correlations between the coefficients in the effective range expansion. These correlations may be regarded as low-energy theorems and are known to hold to a high accuracy in the neutron-proton 3S1 partial wave. We generalize the low-energy theorems to the case of unphysical pion masses and provide results for the correlations between the coefficients in the effective range expansion in this partial wave for pion masses up to Mpi ~ 400 MeV. We discuss the implications of our findings for the available and upcoming lattice-QCD simulations of two-nucleon observables.


I. INTRODUCTION
The dependence of the nuclear forces and, more generally, nuclear observables upon the fundamental parameters of the Standard Model such as the light quark masses is one of the open problems in contemporary theoretical physics. Being a fascinating question on its own, it has attracted considerable attention on the theory side and has been shown to play important role for various phenomena in nuclear, astro-and particle physics. In particular, Bulgac et al. have speculated about an intriguing possibility of the appearance of a bound state in the 3 P 0 neutron-proton partial wave in the so-called chiral limit of Quantum Chromodynamics (QCD) corresponding to massless up and down quarks [1]. Such a scenario would clearly have far reaching consequences for nuclear physics. Another interesting observation was made in Ref. [2], where the proximity of the physical quark masses to the critical trajectory for an infrared renormalization group limit cycle of QCD was pointed out and conjectured to be responsible for the success of Efimov's program for describing the three-nucleon problem, see also Ref. [3] for a follow-up study along this line.
The quark mass dependence of the nuclear force plays an important role for constraining certain types of extensions of the Standard Model. In particular, many such theories allow for the parameters of the Standard Model such as the light quark masses to vary over time. The knowledge of the quark mass dependence of the nuclear force combined with the theory of Big Bang nucleosynthesis (BBN) and the observed nuclear abundances allow one to constrain a possible variation of quark masses at the time of BBN [4,5].
Recently, this topic has experienced renewed interest in connection with the anthropic considerations in nuclear physics [6]. In particular, the second 0 + state of 12 C, the so-called Hoyle state, is well-known to play a crucial role for a formation of life-essential elements such as 12 C and 16 O in hot, old stars due to its closeness to the 4 He-8 Be threshold [7]. Changing its excitation energy of = 379.47 (18) keV, measured with respect to the triple-α threshold, by more than ∼ 25% was shown to strongly decrease the production rate of either 12 C or 16 O [8,9]. The Hoyle state is, therefore, viewed as a promising candidate to address the fine-tuning problem of the fundamental constants of Nature in connection with the anthropic principle. The sensitivity of the energy difference to a variation of the light quark masses was recently analyzed within an ab initio framework of nuclear lattice simulations [10,11]. This theoretical approach makes use of a discretized formulation of chiral effective field theory (EFT) combined with the Monte Carlo method to perform Euclidean time evolution of an A-nucleon state, see Ref. [12] for a review article. Presently, by far the dominant source of the theoretical uncertainty in this calculation is related to the lack of knowledge of the quark mass dependence of the nuclear force or, more precisely, of the nucleon-nucleon (NN) S-wave scattering lengths.
While the ultimate answer to the question of quark mass dependence of hadronic observables is eventually to be arXiv:1504.07852v1 [nucl-th] 29 Apr 2015 provided by lattice QCD, one can gain useful insights into this topic within the framework of chiral EFT. During the past two decades, this theoretical approach was developed into a powerful tool to derive nuclear two-and manybody forces and the corresponding current operators in a systematically improvable and model-independent way, see Refs. [13,14] for recent review articles. It exploits the symmetry and symmetry-breaking pattern of QCD to formulate an appropriate EFT in terms of pions and nucleons (and possibly of the ∆(1232) isobar) which are the relevant degrees of freedom for low-energy nuclear physics in the non-strange sector. The resulting effective Lagrangian has been used to analyze hadronic observables in the Goldstone-boson and single-nucleon sectors as well as nuclear forces and the corresponding current operators by means of the chiral expansion, i.e. a perturbative expansion in powers of the pion mass M π and three-momenta of external particles. In particular, for the two-nucleon force, this expansion has been pushed to next-to-next-to-next-to-next-to-leading order [15,16] and demonstrated to provide very accurate description of low-energy NN scattering observables and the deuteron properties [15,17]. Notice that while the number of the relevant low-energy constants (LECs) in the effective Lagrangian grows with an increasing order of the calculation, the large amount of available experimental data on low-energy proton-proton and neutron-proton scattering observables make their determination unproblematic.
Given its reliance on the spontaneously broken approximate chiral symmetry of QCD, chiral EFT provides, at least in principle, a suitable theoretical framework to address quark mass dependence of the nuclear forces in a systematic way, see Refs. [2-5, 10, 11, 18-24] for recent studies along this line and Refs. [1,25] for related calculations at the more phenomenological level. Clearly, the major complication which prevents one from being able to perform accurate chiral EFT calculations of few-and many-nucleon systems at unphysical values of the quark masses is the lack of knowledge of the corresponding LECs. Quark mass dependence of the short-range interactions induced by integrating out the momentum scale associated with real pion production has been discussed in Ref. [26]. Further complications emerge due to the appearance of a finite cutoff in the calculations. This feature is unavoidable if the one-pion exchange potential is iterated in the non-relativistic Lippmann-Schwinger equation to all orders. Such an approach relies on implicit renormalization which makes it difficult in practice to control the quark mass dependence of short-range operators in a systematic way. Notice that the above-mentioned complications can be avoided by treating the exchange of pions in perturbation theory using e.g. the approach proposed by Kaplan, Savage and Wise (KSW) [27], see Ref. [19] for studies along this line and Ref. [23] for calculations using a closely related approach with dibaryon fields. The perturbative treatment of the one-pion exchange within this framework was, however, shown to be inadequate in spin-triplet channels of NN scattering [28,29].
Recently, an attempt was made to overcome the above mentioned difficulties associated with the unknown quark mass dependence of the short-range part of the nuclear force by employing a resonance saturation hypothesis for contact NN operators [5]. It is well known that the LECs accompanying the NN contact interactions can be understood at a semiquantitative level in terms of the exchange of heavy mesons in the sense of resonance saturation [30]. Using unitarized chiral perturbation theory and the information provided by lattice-QCD simulations to determine the pole positions of the relevant resonances as functions of the quark masses and assuming the validity of the resonance saturation picture at unphysical pion masses allows one to predict the quark mass dependence of the contact operators and to carry out chiral extrapolations in the NN sector. This strategy was followed in Ref. [5] up to next-to-next-to-leading order in chiral EFT within the Weinberg power counting scheme. It was, in particular, found that the deuteron is likely to become less bound at larger values of the light quark mass. These findings are in good agreement with the ones achieved in Ref. [25] using a more phenomenological approach by considering a set of Argonne models for NN interaction. A detailed comparison of these results with the available earlier chiral EFT calculations both within the Weinberg and the KSW power counting schemes can be found in Ref. [5]. We emphasize, however, that modelling the M π dependence of the contact interaction, as done in Ref. [5], provides a substantial source of uncertainty which was not estimated in that work.
Finally, chiral extrapolations of the NN S-wave scattering lengths and the deuteron binding energy were considered in Ref. [31] using the chiral EFT formulation of Ref. [32]. In this approach, the leading-order (LO) NN scattering amplitude is obtained by solving the three-dimensional integral equation introduced originally by Kadyshevsky [33]. This equation is an example of three-dimensional integral equations which satisfy relativistic elastic unitarity. Clearly, the Kadyshevsky equation turns into the usual Lippmann-Schwinger equation upon taking the non-relativistic limit of the two-nucleon propagator. An important feature of the Kadyshevsky equation for the LO NN scattering amplitude is its renormalizability. In particular, all ultraviolet divergences generated by iterations can be explicitly absorbed into redefinition of the NN derivative-less contact interaction. This feature allows one to remove the ultraviolet cutoff when calculating the amplitude and to avoid the above mentioned complications emerging in the non-relativistic framework with a finite cutoff. Assuming that the contributions of the M 2 π dependent short-range operators are suppressed relative to the derivative-less, M 2 π independent contact terms, the quark mass dependence of the S-wave phase shifts at LO emerges entirely from the propagator in the one-pion exchange potential. The resulting parameterfree predictions in the spin-triplet neutron-proton channel are in very good agreement with the calculations reported in Ref. [5]. Also the M π dependence of the 1 S 0 scattering length agrees qualitatively with the one reported in that paper.
While all above-mentioned studies [5,25,31], which are carried out using rather different theoretical approaches, predict less attraction in the deuteron channel for increasing values of the light quark masses, the available lattice-QCD simulations seem to indicate an opposite trend. The first pioneering calculation of NN scattering in the framework of lattice QCD has been carried out two decades ago by Fukujita et al. [34] within the quenched approximation and using pion masses of M π 550 MeV and heavier. Ten years later, the first fully dynamical lattice QCD calculation of the NN S-wave scattering length at pion masses of M π 350 MeV, 490 MeV and 590 MeV were reported in Ref. [35]. This study found smaller values of the 3 S 1 scattering length (corresponding to a stronger bound deuteron) but was inconclusive of its sign. Notice further that similarly to the calculation reported in Ref. [34], the authors of Ref. [35] did not address volume dependence. Recently, binding energies of light nuclei have been calculated in 3-flavor QCD at M π 390 MeV [36] and M π 800 MeV [37] as well as in 2 + 1-flavor QCD at M π 510 MeV [38], see also Ref. [39] for a quenched QCD calculation of the same observables. These fully dynamical lattice QCD calculations found a considerably stronger bound deuteron at these large pion masses, see section IV for more details. Notice further that the NPLQCD Collaboration also extracted the values of the effective range and the first shape coefficient at the pion mass of M π 800 MeV [40]. On the other hand, the HAL QCD Collaboration found no bound state in the NN 3 S 1 -3 D 1 channel by carrying out 3-flavor QCD simulations at pseudoscalar meson masses of 469 . . . 1171 MeV and employing a different approach by making use of the two-nucleon potential at the intermediate step of the calculation [41]. Very recently, Yamazaki et al. have performed 2 + 1-flavor QCD calculations of light nuclei at pion mass of M π 300 MeV [42]. They found the deuteron to be even stronger bound at M π = 300 MeV compared to their earlier study [39] at M π 510 MeV. We further emphasize that the preliminary results of the NPLQCD Collaboration at M π 430 MeV [43] also suggest the deuteron to be stronger bound as compared to the physical case. The current situation with the lattice-QCD results obtained by different groups is thus not completely clear although there seems to be a common trend towards stronger bound deuteron and other light nuclei at values of the light quark masses larger than the physical ones.
In this paper we address the quark mass dependence of low-energy NN scattering observables from the point of view of low-energy theorems (LETs). These theorems establish model-independent relations between the coefficients in the effective range expansion (ERE) which are determined by the long-range part of the nuclear force governed by the pion exchange. We generalize the LETs to unphysical pion masses and discuss their accuracy and the validity range with respect to the variation in M π . We apply the LETs to provide relations between the coefficients in the ERE in the 3 S 1 partial wave for pion masses up to M π 400 MeV. Our results open the way for nontrivial consistency checks of the existing and upcoming lattice-QCD calculations (see also Ref. [24] for a related discussion). As an example of a possible application, we employ a linear with respect to M π extrapolation of the 3 S 1 effective range suggested in Ref. [40] to predict the resulting M π dependence of the scattering length, shape parameters and the deuteron binding energy. The obtained results agree well with the lattice-QCD calculations reported in Refs. [36,38,43]. We also use the LETs to predict the values of the scattering length, effective range and the shape parameters using the deuteron binding energies calculated at M π = 300 MeV [42], M π = 390 MeV [36] and M π = 430 MeV [43] which can be tested on the lattice by calculating the phase shifts.
Our paper is organized as follows. In section II we discuss in detail the meaning of the LETs and provide results at the physical value of the pion mass in the 1 S 0 and 3 S 1 channels. A generalization of the LETs to unphysical values of the pion mass in the 3 S 1 partial wave is considered in section III. Implications of these findings for lattice-QCD results are addressed in section IV. The main results of our work are summarized in section V.

II. LOW-ENERGY THEOREMS FOR THE PHYSICAL PION MASS
Long-range interactions are responsible for the near-threshold left-hand singularities of the partial wave scattering amplitude and control its energy dependence [28,44]. In particular, they impose correlations between the coefficients in the effective range expansion which can be regarded as low-energy theorems. In the following, we discuss in some detail the meaning of the LETs following the lines of Refs. [45][46][47] and using the framework of non-relativistic quantum mechanics which is appropriate for analyzing low-energy NN scattering. For a related discussion on reconstructing the scattering amplitude in the physical region based on the discontinuities across the left-hand cuts and employing the unitarity constraints see Refs. [48,49] and references therein.
Consider two non-relativistic particles of mass m interacting via some short-range potential V . The corresponding S-matrix for an uncoupled channel with the orbital angular momentum l is parametrized in terms of a single phase shift δ l and can be written in terms of the T -matrix as where k denotes the scattering momentum in the center-of-mass system (CMS). In the complex energy plane, the partial-wave scattering amplitude and thus also the T -matrix possess a so-called unitarity cut, a kinematic singularity due to two-body unitarity. The unitarity cut starts from the branch point at the threshold (E = 0) and goes to positive infinity. Furthermore, there are singularities associated with the interaction mechanism and located at the negative real axis. In particular, in the case of the Yukawa potential (∼ exp(−M r)/r) corresponding to an exchange of a meson of mass M , the amplitude has a left-hand cut starting at k 2 = −M 2 /4. Bound and virtual states reside as poles at the negative real axis (k = i|k| and k = −i|k| for bound-and virtual-state poles, respectively) while resonances show up as poles at complex energies. The cut structure of the nucleon-nucleon scattering amplitude in the absence of long-range electromagnetic potentials is visualized in Fig. 1.
It is useful to express the T -matrix in terms of the so-called effective range function F l (k) ≡ k 2l+1 cotδ l (k) via Contrary to the scattering amplitude, the effective range function does not possess the kinematic unitarity cut and can be shown to be a real meromorphic (i.e. analytic except for poles) function of k 2 near the origin k = 0 for non-singular potentials of a finite range [50,51]. It can, therefore, be Taylor-expanded about the origin leading to the well-known effective range expansion (ERE) which has the form where a and r refer to the scattering length and the effective range, respectively, while v i denote the so-called shape parameters. The convergence radius of the ERE is bounded from above by the lowest-lying left-hand singularity associated with the potential. In particular, given that the longest-range part of the strong nuclear force is due to the one-pion-exchange potential (OPEP), the ERE for NN scattering is expected to converge for energies up to |E lab | ∼ M 2 π /(2m N ) = 10.5 MeV. Notice that the actual convergence range of the ERE might be smaller if the effective range function possesses poles corresponding to zeros of the scattering amplitude whose positions are determined by the strength of the interaction. Such a situation emerges, for example, when phase shifts change the sign. The problem of decreasing the applicability range of the ERE can be easily avoided with if the Taylor expansion is replaced by e.g. Pade approximants, see Ref. [52] for a related discussion. Furthermore, we emphasize that a generalization of the above considerations to coupled channels such as e.g. the 3 S 1 -3 D 1 channel in neutron-proton scattering is straightforward and amounts mainly to replacing the amplitude by a 2 × 2 matrix, see Ref. [53] for more details.
The framework of the ERE can be generalized to the case in which the potential is given by a sum of long-range (r L ∼ M −1 L ) and short-range (r S ∼ M −1 S M −1 L ) potentials, V L and V S , respectively. Following van Haeringen and Kok [54], one can define the modified effective range function In this equation, f L l (k) denotes the Jost function defined according to f L l (k) ≡ f L l (k, r) r=0 with f L l (k, r) being the Jost solution of the Schrödinger equation corresponding to the potential V L , i.e. the particular solution that fulfils lim r→∞ e −ikr f l (k, r) = 1 . (2.5) Further, δ L l (k) denotes the phase shift associated with the potential V L and the quantity R L l (k) can be computed from f L l (k, r) as follows: . (2.6) Here and in what follows, the superscript "L" indicates that the quantity of interest can be calculated solely from the long-range part of the potential V L . As proven in Ref. [54], the modified effective range function F M l (k 2 ) does not contain the left-hand singularities associated with the long-range potential and reduces, per construction, to the ordinary effective range function F l (k 2 ) in the case of V L = 0. It is a real meromorphic function in a much larger region set by r −1 S as compared to F l (k 2 ). In particular, for Yukawa-type potentials, the region in which the modified effective range function is meromorphic is set by |k| < M S /2. Similarly to the ERE, one can Taylor expand the function F M l (k 2 ) near the origin via This expansion is referred to as the modified effective range expansion (MERE).
The most frequently used application of the above framework concerns proton-proton scattering. In that case, the long-range interaction is due to the Coulomb potential V L (r) = α/r with α being the fine structure constant. The left-hand cut starts directly at the threshold so that the ERE has zero radius of convergence. Notice further that the Jost solution and, consequently, the function R L l (k) can be calculated analytically for the case of the Coulomb potential. For example, for l = 0 and the repulsive Coulomb potential, the modified effective range function takes the following well-known form: where the Coulomb phase shift is δ C ≡ arg Γ(1 + iη) and the quantity η is given by Further, the functions C 2 0 (η) (the Sommerfeld factor) and h(η) read Here, Ψ(z) ≡ Γ (z)/Γ(z) denotes the digamma function. For more details on the analytic properties of the scattering amplitude and related topics the reader is referred to the review article [55].
We are now in the position to clarify the meaning of the LETs. Assuming that the functions F l (k 2 ) and F M l (k 2 ) do not possess discrete poles in their meromorphic regions, the size of the coefficients in the ERE and MERE (except for the scattering length) is expected to be governed by the scales M L and M S associated with the lowest left-hand singularities, see [44] for a related discussion. If the long-range interaction is known, the quantities f L l (k), R L l (k) and δ L l (k) entering the right-hand side of Eq. (2.4) can be calculated explicitly. Using the MERE, truncated at a given order to approximate the function F M l (k 2 ), and inverting Eq. (2.4) then allows one to express the phase shift δ l (k) in terms of the known quantities f L l (k), R L l (k) and δ L l (k) and the first coefficients in the MERE which parametrize physics associated with the short-range interaction. For example, at LO, the modified effective range function is approximated as F M l (k 2 ) −1/a M . Thus, using a single piece of information about the short-range interaction in form of a M or, equivalently, the usual scattering length 1 a allows one to predict all coefficients in the ERE. These predictions are accurate up to corrections emerging from the second term in the MERE. The MERE for F M l (k 2 ) thus provides a systematically improvable expansion of the subthreshold parameters entering Eq. (2.3) in powers of M L /M S . The resulting correlations between the subthreshold parameters at each fixed order in this expansion are referred to as LETs. For an example of an exactly solvable model, which allows one to get analytical insights into the meaning of the LETs, the reader is referred to Ref. [45].
We are now in the position to address the LETs for S-wave neutron-proton scattering. Here and in what follows, we restrict ourselves to the case of exact isospin symmetry and neglect the long-range electromagnetic interactions between the neutron and proton. The longest-range part of the nuclear force is due to the OPEP which in momentum space has the form where q ≡ p − p is the momentum transfer of the nucleon while σ i (τ i ) denote the spin (isospin) Pauli matrices of the nucleon i. Further, g A and F π refer to the axial vector coupling of the nucleon and pion decay constant, respectively. Throughout this work, we adopt the values of F π = 92.4 MeV, g A = 1.267 for these constants at the physical value of the pion mass M phys π ≡ (2M π ± + M π 0 )/3 = 138.03 MeV. Upon performing Fourier transform, one obtains the well-known coordinate-space expression for the long-range part of the OPEP, namely where S 12 is the usual tensor operator, Notice that the above potential describing point-like nucleons interacting via one-pion exchange is singular and thus cannot be directly employed in the calculation of the LETs using the quantum mechanical framework outlined above. 2 One possible way to overcome this complication is to apply a properly regularized version of the OPEP. For example, one can employ a coordinate-space regularization procedure introduced in Refs. [57,58] or a closely related method used in Refs. [15,17]. These approaches are analogous to the Pauli-Villars regularization and maintain, per construction, the analytic structure of the amplitude in the sense that non-analyticities induced by the regulator appear at momenta k ≥ Λ. Here, Λ refers to the momentum scale associated with the regulator which in such an approach should be chosen of the order of the breakdown scale in the problem [45,59]. The properly regularized OPEP can then be used to compute the Jost solution by solving the Schrödinger equation and to work out the LETs for NN scattering.
Instead of using the approach described above, we employ here the modified version of Weinberg's chiral EFT proposed in Ref. [32] in order to work out the LETs for NN scattering. This momentum-space framework is more convenient for performing numerical calculations of the LETs and, as already mentioned in the introduction, has an advantage of being well suited for carrying out chiral extrapolations beyond the LETs [60].
Within the modified Weinberg EFT approach of Ref. [32], the LO amplitude is obtained by solving the Kadyshevsky equation [33] which, for the case of the fully off-shell kinematics, takes the form where G(p 0 , q) is the free Green function Correlations between the inverse scattering length a −1 , effective range r and the first three shape parameters v2, v3 and v4 induced by the one-pion exchange interaction in the 1 S0 (left panel) and 3 S1 (right panel) channels. Solid rectangles, diamonds, open triangles and circles correspond to the values of r, v2, v3 and v4, respectively, extracted from the Nijmegen partial wave analysis [53,61]. The results in the 3 S1 partial wave correspond to the Blatt-Biedenharn parametrization of the S-matrix [62].
Further, p ( p ) is the incoming (outgoing) three-momentum of the nucleon in the center-of-mass frame, p 0 = k 2 + m 2 N with m N denoting the nucleon mass and k being the corresponding three-momentum of an incoming (on-mass-shell) nucleon. The integral equation (2.13) is solved numerically in the partial wave basis.
In this work, we are interested in the 1 S 0 and the coupled 3 S 1 -3 D 1 channels of neutron-proton scattering. In both cases, the LO potential is given by the OPEP specified in Eq. (2.11) accompanied by a momentum-independent contact interaction, whose strength is adjusted to reproduce the S-wave scattering length. We regularize the integral equation by introducing an ultraviolet momentum cutoff Λ and take the limit Λ → ∞ when calculating the phase shifts. Notice that contrary to the non-relativistic framework based on the Lippmann-Schwinger equation, this is a legitimate procedure since the LO integral equation is renormalizable in the sense that all divergences emerging from its iterations are absorbable into redefinition of the LO contact interaction. Once the phase shifts are calculated, we extract the effective range and the shape parameters numerically by matching the truncated ERE to the effective range function. Given that the calculated S-wave scattering amplitudes fulfill elastic unitarity, correctly reproduce the first left-hand cut due to the OPEP and match the experimental value at the threshold (given by the scattering length), the resulting predictions for r and v i are equivalent to the LO LETs discussed above in the context of the MERE. We emphasize, however, that employing the cutoff-independent framework of Ref. [32] is not crucial for testing the LETs in NN scattering. A more conventional non-relativistic chiral EFT approach based on the Lippmann-Schwinger equation and utilizing a finite ultraviolet cutoff is equally well suited for this purpose.
The LO LETs for the 1 S 0 and the 3 S 1 channels have already been addressed in Ref. [32]. In Fig. 2, we show the LO LET predictions for the effective range and the first three shape parameters as functions of the inverse scattering length which is an input parameter in our calculations. Here and in what follows, the results in the 3 S 1 partial wave correspond to the Blatt-Biedenharn parametrization of the S-matrix [62]. To avoid the appearance of large/small numbers, we use dimensionless quantities by multiplying the coefficients in the ERE with the corresponding powers of the pion mass. As expected, the predicted dimensionless coefficients appear to be of a natural size for a rather large range of values of the scattering length. Notice further that the limit a → ±∞, which describes the situation in which there is a zero-energy bound state, does not lead to any peculiarities in r and v i . In both channels, one observes a clear tendency for r and v i to become unnaturally large in magnitude in the regions of 1/(aM π ) −0.5 and 1/(aM π ) 1. We found that the unnaturally large values of r and v i for 1/(aM π ) 1 are due to the appearance of a pole in the effective range function at relatively low positive energy. The pole position corresponds to the energy at which the corresponding phase shift crosses zero. For example, for 1/(aM π ) 1.4, the pole in the 3 S 1 effective range I: Various available results for the low-energy theorems for the neutron-proton 1 S0 and 3 S1 partial waves as described in the text. The results in the 3 S1 partial wave correspond to the Blatt-Biedenharn parametrization of the S-matrix [62]. function is located at k 60 MeV which is already within its meromorphic region. For larger values of 1/(aM π ), the pole moves towards the threshold so that the values of the r and v i are actually governed by the pole position rather than by M π . As already emphasized above, the appearance of a pole in the effective range function can be easily accommodated by replacing the Taylor expansion of F l (k 2 ) by e.g. Pade approximation. On the other hand, the case 1/(aM π ) −0.5 corresponds to the absence of a bound state and, given the natural value of the scattering length, it might already describe a perturbative regime of NN scattering. For weakly interacting systems, it is easy to see analytically (by making use of the Born approximation) that the coefficients in the ERE do not scale anymore with the range of the interaction. This is due to the appearance of an additional dimensionless small parameter associated with the weakness of the interaction.
The predictions for r and v i based on the LO LETs and corresponding to the experimentally observed values of the scattering lengths are compared in Fig. 2 with the empirical numbers from Ref. [53,61] based on the Nijmegen partial wave analysis (PWA). We also collect various available results for the LETs in the 1 S 0 and the 3 S 1 channels in Table  I. The results for the LO LETs shown in the table are taken from Ref. [32]. While the predicted values of r and v i do reproduce the empirical values in the 3 S 1 channel rather accurately, the agreement in the 1 S 0 partial wave is, at best, on a qualitative level. This pattern can be understood as follows. Generally, the accuracy of the LO LET is expected to be set by the ratio E L /E S , where E L (E S ) denotes the energy corresponding to the branch point of the left-hand cut due to the long-range interaction (branch point of the first left-hand cut in the scattering amplitude which is not correctly described by the employed approximation for the interaction). Therefore, given that we do not take into account the two-pion exchange potential in our calculations and thus cannot describe correctly the second left-hand cut, one may expect the LO LETs to be accurate at the level of ∼ 25%. In fact, the accuracy of the LO LETs in the 3 S 1 partial wave appears even to be somewhat higher. 3 This is presumably due to the fact, that we actually include the parts of the second and higher left hand cuts (see Fig. 1) which are associated with iterations of the OPEP. In fact, according to the chiral power counting, such iterative contributions to the scattering amplitude are expected to be more important than the irreducible two-and more-pion exchange potentials. On the other hand, the low accuracy of the LETs in the 1 S 0 channel is a consequence of the weakness of the OPEP projected on that partial wave. The strong tensor part of the OPEP, which accompanies the operator S 12 in Eq. (2.12), does not contribute to spin-singlet partial waves. The weakness of the OPEP in the 1 S 0 channel is evidenced by the fact, that the phase shift corresponding to the long-range OPEP reaches only about 12.6 • at maximum, which has to be compared with max(δ 1S0 (k)) ∼ 65 • in the real world.
We also list in Table I, for the sake of completeness, the predictions obtained within the KSW framework [27] which relies upon a perturbative treatment of the OPEP. These results are taken from Ref. [28] and show very large deviations from the empirical values. This shows that the perturbative treatment of the OPEP in these channels is not appropriate.
In Ref. [63], the subleading contact interaction C 2 (p 2 + p 2 ) has been taken into account in the 1 S 0 channel within the framework of Ref. [32], both perturbatively and nonperturbatively. Given that the long-range interaction included in that work is entirely given by the OPEP, the resulting predictions for the shape parameters, which are listed in Table  I for the case of the nonperturbative treatment of C 2 , are equivalent to the next-to-leading-order (NLO) LETs. Notice that the non-perturbative inclusion of the subleading contact interaction necessarily results in a residual subtraction scale dependence of the amplitude which is why the predictions for v 2,3,4 are shown within a range of values, see Ref. [63] for more details. While a significant improvement is observed in the description of all v i , the convergence is clearly rather slow in this channel. We find that the accuracy of the LETs in the 1 S 0 channel is insufficient for the purpose of providing constraints on the available lattice QCD results. For this reason, we will concentrate in the following entirely on the spin-triplet channel.
While the LETs in the 3 S 1 channel appear to be fairly accurate already at LO as far as one is concerned with predicting the values of r and v i from the scattering length, their accuracy is still insufficient for extracting the deuteron binding energy using the effective range as input, as we intend to do in the next sections. This is because of the function r(a −1 ) being rather flat in the relevant region of values of 1/(aM π ), which tends to magnify the uncertainty in the value of a −1 extracted from r. Notice that the observed nearly quadratic dependence of r on a −1 may be related to the universal behavior of this function found for van der Waals-like interactions [65], see also Ref. [66] for a related discussion. In addition, the deuteron binding energy depends, to a good accuracy, quadratically on a −1 , so that the error in the value of the binding energy is further magnified. For example, using the LO LETs together with the experimentally observed value of the effective range r = 1.75 fm as input, one extracts a = 7.2 fm and B d = 1.1 MeV which is two times smaller than the experimentally observed value of B d 2.22 MeV. To deal with this issue, we extend the LETs in the 3 S 1 channel to NLO by taking into account a subleading short-range interaction whose strength is adjusted to reproduce the experimentally observed value of the effective range. Notice that contrary to the 1 S 0 channel, where the long-range part of the OPEP is non-singular and a non-perturbative inclusion of the subleading contact interaction can be carried out semi-analytically in a close analogy to pionless EFT, see Ref. [63] for more details, we are unable to treat the subleading contact interactions non-perturbatively in spin-triplet channels without, at the same time, destroying the explicit renormalizability feature of the scattering amplitude. A perturbative inclusion of the NLO short-range terms will be discussed in detail in a subsequent publication [60]. Here, in order to avoid possible shortcomings due to the perturbative treatment of the subleading short-range terms, we follow a different approach and employ resonance saturation to model higher-order contact interactions by means of a heavy-meson-exchange. In particular, we use for the NLO potential where the heavy-meson mass M is set M = 700 MeV and the strength β is adjusted to reproduce the empirical value of the effective range. In agreement with the arguments given above, we have verified that the results for the LETs are insensitive to details of the short-range interaction. In particular, the deviation in the results caused by using the sigma-like scalar potential instead of the tensor one in Eq. (2.15) for the sub-leading short-range interaction is negligibly small. The sub-leading short-range interaction is only needed to describe the second term in the MERE ∼ r M which appears, in addition to a M , as an input parameter in our calculation. The employed form of the NLO potential ensures that the corresponding scattering amplitude is renormalizable so that the NLO calculations can be carried out in the same way as at LO. The predicted values of the shape parameters after fixing the values of C 0 and β are listed in Table I and, as expected, show a clear improvement as compared with the LO results. In fact, the NLO LETs appear to be accurate at the level of a few percent (except for v 2 ).

III. LOW-ENERGY THEOREMS FOR UNPHYSICAL PION MASSES
The LETs discussed in the previous section can be straightforwardly generalized to the case of non-physical quark (or pion) masses. The main effect from changing the pion mass, as far as the analytic structure of the amplitude is concerned, corresponds to shifts of the branch points of the left-hand cuts due to one-, two-and more-pion exchange, see Fig. 1. This effect originates from the explicit pion mass dependence in the pion propagator. Next, one needs to account for the change in the discontinuity across the left-hand cuts caused by the M π dependence of the strength of the OPEP, namely the ratio g A /F π . To account for this effect, we make use of the lattice-QCD results for the M π dependence of these quantities, see Fig. 3. In particular, we performed quadratic polynomial regression fits (as functions of M 2 π ) of the results from the 2 + 1 flavor simulations by the Budapest-Marseille-Wuppertal (BMW) collaboration for F π [67] and from the 2 and 2 + 1 + 1 flavor simulations reported in Ref. [69] for g A , see also Ref. [70] for recent lattice-QCD calculations of this quantity. In our fits, we included lattice-QCD data for pion masses up to M π = 500 MeV. Further, we enforced the fits to go through the experimental value for a given quantity at the physical value of the pion mass. The results of our fits together with the 67% confidence region are shown in Fig. 3. Notice that the observed tendency of lattice-QCD results for the axial charge of the nucleon to somewhat underestimate the experimental value may indicate that certain systematic corrections have not been properly taken into account in these studies. We expect this issue to be clarified in the near future, when more results near the physical point will become available and the statistical uncertainties of the simulations will be improved. Notice further that we did not include all available lattice QCD results in our fits since this will unlikely improve the accuracy of the interpolations. Especially for the pion decay constant, the uncertainty of the interpolation of the lattice-QCD results is very small and can be safely neglected at the accuracy level of our calculations.
Next, one also has to account for the pion mass dependence of the nucleon mass. Here, we follow the same strategy as for F π and g A and employ a quadratic interpolation of the recent lattice-QCD results for m N (M 2 π ) of the BMW collaboration [67], see Fig. 3. We emphasize that our fits of the M π dependence of F π , g A and m N are not intended to provide correct chiral extrapolations of these quantities to the chiral limit but serve exclusively to enable a smooth interpolation of the available lattice-QCD results.
Using the above results for the pion mass dependence of F π , g A and m N together with the explicit M π dependence in the pion propagator allows us to describe the discontinuity across the first left-hand cuts due to the OPEP for arbitrary values of the pion masses. Thus, we can immediately generalize the LO LETs by calculating r and v i as functions of the (inverse) scattering length at unphysical values of the pion mass. The resulting predictions are shown by the various thick lines in Fig. 4 for M π = 50, 100, 150, 200, 300 and 400 MeV. Notice that contrary to chiral EFT extrapolations, we do not make here any assumptions about the quark mass dependence of the included derivative-less short-range interaction ∝ C 0 , which is traded into the dependence of the calculated observables on the value of the scattering length. In particular, we perform calculations with different pion masses as if we lived in different worlds, in which the physical value of the pion mass were M π = 50, 100, . . . MeV. Thus for each given pion mass the value of C 0 is adjusted to reproduce the given value of the scattering length used as input. In this sense, our results are more general than the ones obtained for the 1 S 0 and 3 S 1 channels in Ref. [31] at LO within the modified Weinberg approach, where the dependence of the short-range interaction on M π was neglected in agreement with the estimation Yamazaki et al. NPLQCD NPLQCD, prelim.

FIG. 4:
Correlations between the inverse scattering length a −1 , effective range r and the first three shape parameters v2, v3 and v4 in the 3 S1 partial wave induced by the one-pion exchange interaction. Various thick lines show the predictions of the LO LETs while light-shaded bands (hardly visible for small Mπ) between thin lines depict the results of NLO LETs and reflect the estimated uncertainty due to unknown Mπ dependence of the subleading short-range interaction as explained in the text.
based on naive dimensional analysis (NDA).
In order to extend these results to NLO, we need to specify the M π dependence of the subleading short-range interaction which is modeled via resonance saturation, see Eq. (2.15). While this effect would be suppressed in the chiral EFT approach if one assumes the validity of NDA for short-range interactions, we decided to estimate it in order to stay as model independent as possible. Specifically, we assume that the M 2 π dependence of the strength β of the short-range interaction is within the envelope built by the lines which go through the physical point and describe a ±50% change in the value of β for M π = 500 MeV, i.e.: with δβ = 0.5 and ∆M 2 π ≡ (M 2 π − (M phys π ) 2 ) Mπ=500 MeV . Note that such an estimation is justified for the considered quantities such as g A , F π and m N , as can be seen from Fig. 3.
Our NLO LETs predictions for the effective range and the shape parameters viewed as functions of the inverse scattering length are visualized in Fig. 4 by the light-shaded bands. These bands correspond to the variation of β at a given value of M π according to Eq. (3.16). For pion masses below M π 200 MeV, one observes a rather small difference between the predictions based on LO and NLO LETs whose size can be viewed as an estimation of the accuracy of the LO LETs. Also, the employed variation of β is essentially invisible for such values of the pion mass. For heavier pions, both the differences between the LO and NLO LETs as well as the uncertainty in the results associated with the M π dependence of β start to increase. As shown in this figure, one cannot expect to have accurate predictions for pion masses above M π 400 MeV. Notice that decrease in the predictive power of the LETs is to be expected for heavier pions due to the decreasing separation between the soft and hard scales in the problem.   [40]. Right panel: Linear with M 2 π interpolation of the quantity Mπr in the 3 S1 partial wave according to Eq. (4.18). In both cases, solid squares refer to the value of Mπr at the physical point.

IV. IMPLICATIONS FOR LATTICE QCD CALCULATIONS
We are now in the position to confront the LETs with the available lattice QCD results in the NN sector. In table II, we list the published lattice-QCD results for the S-wave scattering parameters and energies of the bound states together with the experimental data. We do not show in the table the results from Ref. [35] where volume dependence was not addressed. Unfortunately, lattice calculations in the NN sector focused so far mainly on the binding energies and do not provide information on the scattering parameters. An exception is the work of Ref. [40], which provides, in addition to the binding energies, also the values of the scattering length, effective range and even the first shape parameter at the pion mass of M π 800 MeV. Clearly, such heavy pion masses are out of reach of the LETs. On the other hand, the authors of Ref. [40] conjectured that the quantity M π r may exhibit a nearly linear dependence on the pion mass. The suggested linear interpolation between the physical point and the lattice result has the form [40]: shape parameters and the deuteron binding energy. Our results for the deuteron binding energy B d , the ratio γ d /M π , where γ d = √ B d m N is the deuteron binding momentum, the ratio a/r and the first three shape parameters M 3 π v 2 , M 5 π v 3 and M 7 π v 4 are visualized in Fig. 6. In this figure, the dark shaded bands result from the variation of the constant β specified in Eq. (3.16) and reflect the uncertainty of the NLO LETs 4 . The light-shaded bands correspond to the resulting uncertainty which emerges from the theoretical uncertainty at NLO and the errors of the linear interpolation of M π r(M π ) (see the left panel of Fig. 5 and Eq. (4.17)) added in quadrature. Notice that we also show in Fig. 6 the preliminary lattice-QCD result of the NPLQCD collaboration at M π = 430 MeV [43]. Remarkably, the linear M π dependence of M π r suggested in Ref. [40] indeed appears to describe very well the common trend of the lattice-QCD results for the deuteron binding energy at intermediate pion masses. Also the NPLQCD Collaboration results of Ref. [40] for B d , γ d /M π , a/r and M 3 π v 2 at the pion mass of M π = 800 MeV can be well described by further extrapolating our results to heavier pion masses without introducing any strong curvature. Assuming the validity of Eq. (4.17) for pion masses below the physical one, we conclude that the deuteron becomes unbound for M π ∼ 50 MeV. It is also interesting to notice that the scattering length and the shape parameters show rather strong variations with the pion mass around and below the physical point. This nontrivial behavior is driven by the long-range physics associated with the pion exchange and is, in principle, testable in lattice QCD. The obtained results for the quantities γ d /M π and a/r, which probe the amount of fine tuning in the NN system, suggest that the physically realized value of the quark mass is close to the point, which separates the strong fine-tuning regime characterized by the rapidly growing scattering length from the regime featuring a fairly small amount of fine tuning with a/r = 2 . . . 3 within the large range of pion masses.
We emphasize that the observed agreement between the predicted M π dependence of the deuteron binding energy and lattice-QCD results is a rather nontrivial consequence of the assumed extrapolation of M π r. To illustrate this point, we consider an alternative scenario by assuming that M π r is a linear function of M 2 π rather than of M π . Notice that while a linear in M π dependence seems to be natural for NN observables, such a scenario cannot be excluded a priori. Using the lattice-QCD result of Ref. [40] for M π r at the pion mass of M π = 800 MeV, the resulting interpolation formula takes the form M π r ∼ = C ( 3 S1) + D ( 3 S1) M 2 π , where C ( 3 S1) = 1.149 +0.009  Fig. 7. The assumed quadratic dependence on M π is clearly in conflict with the lattice-QCD results and, therefore, appears to be highly unlikely. Interestingly, the shape parameters seem to be more robust with respect to the variation of the functional form of r(M π ) and show a qualitatively similar behavior in both scenarios.
Finally, as a last application of the LETs, we use the quoted values of the deuteron binding energy to predict the scattering length, effective range and the first three shape parameters at pion masses of M π = 300 MeV, M π = 390 MeV and M π = 430 MeV where lattice-QCD results are already available. Our results for the coefficients in the ERE are collected in Table III. In the case of the preliminary NPLQCD calculation at M π = 430 MeV, the values for the deuteron binding energy and the errors are extracted from a corresponding plot of Ref. [43]. These results may serve as useful consistency checks for the already published lattice-QCD results if the corresponding lattice data can be used to extract, in addition to the deuteron binding energy, also information about the scattering parameters.

V. SUMMARY AND CONCLUSIONS
The quark mass dependence of the low-energy NN scattering observables is investigated using the low-energy theorems which establish model independent relations between the coefficients in the effective range expansion governed by the long-range one-pion exchange potential. In order to clarify the meaning of the LETs, we employed the modified III: Predictions for the scattering length, effective range and shape parameters based on the lattice-QCD results for the deuteron binding energy. The first two errors correspond to the ones of the lattice-QCD results for the binding energy. The last quoted error is our estimation of the uncertainty of the NLO LETs associated with the unknown Mπ dependence of the subleading short-range interaction as described in the text. effective range expansion [54] to parametrize the short-range physics in a systematic way. Using the knowledge of the long-range interaction together with the MERE truncated at a given order allows one to calculate all coefficients of the ERE within a systematically improvable expansion in powers of the ratio of the long-and short-range scales. The explicit treatment of the left-hand cuts in the partial-wave scattering amplitudes due to the long-range pionic forces results in correlations between the coefficients in the ERE which are regarded as the LETs.
In this work we considered the LETs for NN scattering in the 3 S 1 and 1 S 0 channels. At leading order in the above mentioned expansion, the short-range physics is parameterized solely in terms of the modified scattering length. Therefore, in order to obtain LETs at LO, we considered two nucleons interacting via the OPEP supplemented by two contact interactions without derivatives, whose strengths were adjusted to reproduce the scattering length in each of these channels. The NN scattering amplitude was calculated by solving the Kadyshevsky equation which is exactly renormalizable at this order [32]. In the 3 S 1 channel, the predicted values of the effective range and the shape parameters agree rather well with the empirical values extracted from the Nijmegen PWA [64] already at LO. On the other hand, the OPEP projected onto the 1 S 0 partial wave is very weak, so that the LETs show much less predictive power in that channel. We, therefore, restricted ourselves to the spin-triplet channel in this paper. To further increase the accuracy of our predictions, we extended the calculations to NLO LETs by including the subleading short-range interaction and tuning its strength to reproduce the effective range. The resulting predictions for the shape parameters in the 3 S 1 partial wave are found to be accurate at the level of a few percent.
As a next step, we generalized the LETs in the 3 S 1 channel to study correlations between various low-energy observables at unphysical pion masses up to M π 400 MeV. To this aim, we made use of the available lattice-QCD results for the pion mass dependence of the pion decay constant, nucleon axial vector coupling and the nucleon mass. The predicted correlations between the low-energy NN parameters at unphysical pion masses open the way for nontrivial consistency checks of the ongoing and upcoming lattice QCD calculations in the NN sector and might also be useful for reducing the systematic uncertainty.
As an application, we extracted the low-energy parameters, namely the scattering length, effective range and first three shape parameters, from the lattice-QCD calculations of the deuteron binding energy at the pion masses M π = 300 MeV [42], 390 MeV [36] and 430 MeV [43]. Our predictions illustrate the importance of a simultaneous extraction of the scattering parameters in addition to the binding energies in lattice-QCD simulations of the NN system.
As another application, we applied the LETs to predict the M π behavior of the binding energy, scattering length and shape parameters based on the linear M π dependence of the effective range M π r suggested in Ref. [40]. We found that the predicted shape of the binding energy as a function of M π is in good agreement with the general trend of the lattice QCD data. We also found that the scattering length and the shape parameters show rather strong variations with the pion mass around and below the physical point. This nontrivial behavior is driven by the long-range physics associated with the pion exchange and is, in principle, testable in lattice QCD.
While our analysis establishes correlations between NN low-energy parameters, it does not provide extrapolation of lattice-QCD results to the physical point. The possibility of carrying out chiral extrapolations of NN scattering observables and the deuteron binding energy by extending the chiral EFT analysis of Ref. [31] beyond the LO will be discussed in a subsequent publication [60].