Neutron-proton scattering at next-to-next-to-leading order in Nuclear Lattice Effective Field Theory

TL;DR

This paper systematically studies neutron-proton scattering using Nuclear Lattice Effective Field Theory up to NNLO, demonstrating small lattice artifacts at fine lattice spacings and good agreement with experimental phase shifts.

Contribution

It provides a non-perturbative NNLO analysis of neutron-proton scattering in NLEFT with a radial Hamiltonian method, including uncertainty estimates and lattice artifact assessments.

Findings

NNLO results agree with Nijmegen partial-wave analysis at a = 0.98 fm

Lattice artifacts are small at finer lattice spacings

Physics below 100 MeV CM momentum is lattice-spacing independent

Abstract

We present a systematic study of neutron-proton scattering in Nuclear Lattice Effective Field Theory (NLEFT), in terms of the computationally efficient radial Hamiltonian method. Our leading-order (LO) interaction consists of smeared, local contact terms and static one-pion exchange. We show results for a fully non-perturbative analysis up to next-to-next-to-leading order (NNLO), followed by a perturbative treatment of contributions beyond LO. The latter analysis anticipates practical Monte Carlo simulations of heavier nuclei. We explore how our results depend on the lattice spacing a, and estimate sources of uncertainty in the determination of the low-energy constants of the next-to-leading-order (NLO) two-nucleon force. We give results for lattice spacings ranging from a = 1.97 fm down to a = 0.98 fm, and discuss the effects of lattice artifacts on the scattering observables. At a = 0.98 fm, lattice artifacts appear small, and our NNLO results agree well with the Nijmegen partial-wave analysis for S-wave and P-wave channels. We expect the peripheral partial waves to be equally well described once the lattice momenta in the pion-nucleon coupling are taken to coincide with the continuum dispersion relation, and higher-order (N3LO) contributions are included. We stress that for center-of-mass momenta below 100 MeV, the physics of the two-nucleon system is independent of the lattice spacing.

Tables

Table 1. Table 1: Hopping coefficients ω i subscript 𝜔 𝑖 \omega_{i} for the free nucleon action, for different levels of improvement.

	unimproved	$𝒪\left(a^2\right)$ improved	$𝒪\left(a^4\right)$ improved
${\omega }_0$	1	$5/4$	$49/36$
${\omega }_1$	1	$4/3$	$3/2$
${\omega }_2$	0	$1/12$	$3/20$
${\omega }_3$	0	0	$1/90$

Table 2. Table 2: Summary of lattice spacings a 𝑎 a (spatial) and a t subscript 𝑎 𝑡 a_{t} (temporal) and box dimensions L 𝐿 L . The physical spatial lattice volume V 𝑉 V is kept constant at ( L ​ a ) 3 ≃ ( 63 ​ fm ) 3 similar-to-or-equals superscript 𝐿 𝑎 3 superscript 63 fm 3 (La)^{3}\simeq(63~{}\mathrm{fm})^{3} .

$a^{-1}$ [MeV]	$a_t^{-1}$ [MeV]	$a$ [fm]	$L$	$La$ [fm]
100	150	1.97	32	63.14
120	216	1.64	38	62.48
150	337.5	1.32	48	63.14
200	600	0.98	64	63.14

Table 3. Table 3: Summary of the fitting procedure, indicating which parameters are fitted to what scattering channel at each order in NLEFT, and the resulting χ 2 / N dof superscript 𝜒 2 subscript 𝑁 dof \chi^{2}/N_{\mathrm{dof}} (for a = 0.98 𝑎 0.98 a=0.98  fm).

order	fit channels	fit parameters	${\chi }^2/N_{\mathrm{dof}}$
LO	${}_{}{}^{1}S_0^{}$, ${}_{}{}^{3}S_1^{}$	$C_{{}_{}{}^{1}S_0^{}}$, $C_{{}_{}{}^{3}S_1^{}}$, $b_s$	$30.38$
NLO	${}_{}{}^{1}S_0^{}$	$C_{{}_{}{}^{1}S_0^{}}$, ${\tilde{C}}_1$	$1.77$
	${}_{}{}^{3}S_1^{}$, ${ϵ}_1$	$C_{{}_{}{}^{3}S_1^{}}$, ${\tilde{C}}_2$, ${\tilde{C}}_3$	$88.81$
	${}_{}{}^{1}P_1^{}$	${\tilde{C}}_4$	$11.94$
	${}_{}{}^{3}P_0^{}$, ${}_{}{}^{3}P_1^{}$, ${}_{}{}^{3}P_2^{}$	${\tilde{C}}_5$, ${\tilde{C}}_6$, ${\tilde{C}}_7$	$6.51$
NNLO	${}_{}{}^{1}S_0^{}$	$C_{{}_{}{}^{1}S_0^{}}$, ${\tilde{C}}_1$	$0.36$
	${}_{}{}^{3}S_1^{}$, ${ϵ}_1$	$C_{{}_{}{}^{3}S_1^{}}$, ${\tilde{C}}_2$, ${\tilde{C}}_3$	$28.81$
	${}_{}{}^{1}P_1^{}$	${\tilde{C}}_4$	$2.79$
	${}_{}{}^{3}P_0^{}$, ${}_{}{}^{3}P_1^{}$, ${}_{}{}^{3}P_2^{}$	${\tilde{C}}_5$, ${\tilde{C}}_6$, ${\tilde{C}}_7$	$25.59$

Table 4. Table 4: Fitted constants and low-energy parameters for a = 0.98 𝑎 0.98 a=0.98  fm. The LO constants C S 0 1 subscript 𝐶 superscript subscript 𝑆 0 1 C_{{}^{1}S_{0}} and C S 1 3 subscript 𝐶 superscript subscript 𝑆 1 3 C_{{}^{3}S_{1}} are given in units of [ 10 − 4 superscript 10 4 10^{-4} MeV -2 ], and the C i subscript 𝐶 𝑖 C_{i} of the NLO interaction in units of [ 10 − 8 superscript 10 8 10^{-8} MeV -4 ]. Due to the large lattice ( L = 64 𝐿 64 L=64 ) for a = 0.98 𝑎 0.98 a=0.98  fm, an uncertainty analysis using the variance-covariance matrix as in Table  5 was numerically unfeasible. Hence, an estimated uncertainty of 2 % percent 2 2\% has been assigned, which is consistent with the uncertainties for larger a 𝑎 a . For entries with a dagger ( † † \dagger ), the deuteron energy E d subscript 𝐸 𝑑 E_{d} has been included as an additional constraint.

	LO	NLO	NNLO
$C_{{}_{}{}^{1}S_0^{}}$	$-0.101\left(2\right)$	$-0.105\left(2\right)$	$-0.106\left(2\right)$
$C_{{}_{}{}^{3}S_1^{}}$	$-0.118\left(2\right)$	$-0.087\left(2\right)$	$-0.088\left(2\right)$
$b_s$	$0.399\left(8\right)$	$-$	$-$
$C_{q^2}$	$-$	$0.00440\left(8\right)$	$0.135\left(2\right)$
$C_{I^2,q^2}$	$-$	$0.0373\left(8\right)$	$0.0303\left(6\right)$
$C_{S^2,q^2}$	$-$	$-0.0292\left(6\right)$	$-0.0301\left(6\right)$
$C_{S^2,I^2,q^2}$	$-$	$-0.00190\left(4\right)$	$-0.0254\left(5\right)$
$C_{{\left(q\cdot S\right)}^2}$	$-$	$0.0378\left(8\right)$	$0.0360\left(7\right)$
$C_{I^2,{\left(q\cdot S\right)}^2}$	$-$	$0.00200\left(4\right)$	$0.0212\left(4\right)$
$C_{\left(q\times S\right)\cdot k}^{I=1}$	$-$	$0.0150\left(3\right)$	$0.0165\left(3\right)$
$E_d$ [MeV]	$2.16\left(4\right)$	$2.22{\left(4\right)}^{†}$	$2.22{\left(4\right)}^{†}$
$r_{{}_{}{}^{1}S_0^{}}$ [fm]	$2.12\left(4\right)$	$2.50\left(5\right)$	$2.63\left(5\right)$
$a_{{}_{}{}^{1}S_0^{}}$ [fm]	$-22.5\left(4\right)$	$-23.4\left(5\right)$	$-23.7\left(5\right)$
$r_{{}_{}{}^{3}S_1^{}}$ [fm]	$1.73\left(3\right)$	$1.70\left(3\right)$	$1.74\left(3\right)$
$a_{{}_{}{}^{3}S_1^{}}$ [fm]	$5.4\left(1\right)$	$5.4\left(1\right)$	$5.4\left(1\right)$

Table 5. Table 5: Fitted constants and low-energy S 𝑆 S -wave parameters for a = 0.98 𝑎 0.98 a=0.98  fm. The LO constants C S 0 1 subscript 𝐶 superscript subscript 𝑆 0 1 C_{{}^{1}S_{0}} and C S 1 3 subscript 𝐶 superscript subscript 𝑆 1 3 C_{{}^{3}S_{1}} are given in units of [ 10 − 4 superscript 10 4 10^{-4} MeV -2 ], and the C i subscript 𝐶 𝑖 C_{i} of the NLO interaction in units of [ 10 − 8 superscript 10 8 10^{-8} MeV -4 ]. The smearing parameter b s subscript 𝑏 𝑠 b_{s} of the LO contact interactions is determined by the LO fit, and thereafter kept fixed at NLO and NNLO. The values in parentheses are the uncertainties calculated using the variance-covariance matrix according to Eq. ( 115 ).

order	fit parameters	$a=1.97$ fm	$a=1.64$ fm	$a=1.32$ fm
LO	$C_{{}_{}{}^{1}S_0^{}}$	$-0.4676\left(2\right)$	$-0.3290\left(7\right)$	$-0.201\left(5\right)$
	$C_{{}_{}{}^{3}S_1^{}}$	$-0.6377\left(2\right)$	$-0.4482\left(2\right)$	$-0.265\left(5\right)$
	$b_s$	$0.0524\left(2\right)$	$0.0917\left(2\right)$	$0.173\left(6\right)$
NLO	$C_{{}_{}{}^{1}S_0^{}}$	$-0.5\left(1\right)$	$-0.35\left(2\right)$	$-0.220\left(2\right)$
	$C_{{}_{}{}^{3}S_1^{}}$	$-0.44\left(7\right)$	$-0.21\left(1\right)$	$-0.152\left(4\right)$
	$C_{q^2}$	$-0.05\left(3\right)$	$-0.032\left(9\right)$	$-0.006\left(1\right)$
	$C_{I^2,q^2}$	$0.08\left(2\right)$	$0.075\left(2\right)$	$0.052\left(1\right)$
	$C_{S^2,q^2}$	$-0.06\left(3\right)$	$-0.046\left(3\right)$	$-0.0341\left(7\right)$
	$C_{S^2,I^2,q^2}$	$0.03\left(2\right)$	$0.029\left(2\right)$	$0.0081\left(2\right)$
	$C_{{\left(q\cdot S\right)}^2}$	$0.11\left(2\right)$	$0.091\left(4\right)$	$0.0553\left(2\right)$
	$C_{I^2,{\left(q\cdot S\right)}^2}$	$-0.11\left(2\right)$	$-0.074\left(4\right)$	$-0.0240\left(8\right)$
	$C_{\left(q\times S\right)\cdot k}^{I=1}$	$0.037\left(8\right)$	$0.026\left(4\right)$	$0.019\left(2\right)$
NNLO	$C_{{}_{}{}^{1}S_0^{}}$	$-0.5\left(1\right)$	$-0.33\left(4\right)$	$-0.21\left(2\right)$
	$C_{{}_{}{}^{3}S_1^{}}$	$-0.5\left(1\right)$	$-0.22\left(1\right)$	$-0.15\left(2\right)$
	$C_{q^2}$	$0.08\left(3\right)$	$0.093\left(7\right)$	$0.118\left(7\right)$
	$C_{I^2,q^2}$	$0.07\left(2\right)$	$0.0668\left(4\right)$	$0.045\left(4\right)$
	$C_{S^2,q^2}$	$-0.06\left(3\right)$	$-0.05\left(2\right)$	$-0.036\left(7\right)$
	$C_{S^2,I^2,q^2}$	$0.01\left(2\right)$	$0.005\left(3\right)$	$-0.014\left(4\right)$
	$C_{{\left(q\cdot S\right)}^2}$	$0.10\left(3\right)$	$0.086\left(7\right)$	$0.056\left(4\right)$
	$C_{I^2,{\left(q\cdot S\right)}^2}$	$-0.10\left(3\right)$	$-0.055\left(4\right)$	$-0.006\left(4\right)$
	$C_{\left(q\times S\right)\cdot k}^{I=1}$	$0.031\left(8\right)$	$0.025\left(4\right)$	$0.018\left(2\right)$

Table 6. Table 6: Low-energy S 𝑆 S -wave parameters, as a function of the lattice spacing a 𝑎 a and the order of the NLEFT expansion. E d subscript 𝐸 𝑑 E_{d} is the deuteron binding energy, and the a i subscript 𝑎 𝑖 a_{i} and r i subscript 𝑟 𝑖 r_{i} denote the scattering lengths and effective ranges in channel i 𝑖 i . The experimental value of E d subscript 𝐸 𝑑 E_{d} is from Ref.  VanDerLeun:1982bhg , and the scattering lengths and effective ranges are from Ref.  Machleidt:2000ge . For entries marked with a dagger ( † † \dagger ), the empirical deuteron energy E d subscript 𝐸 𝑑 E_{d} has been included in the fit as an additional constraint.

order	$a$ [fm]	$E_d$ [MeV]	$r_{{}_{}{}^{1}S_0^{}}$ [fm]	$a_{{}_{}{}^{1}S_0^{}}$ [fm]	$r_{{}_{}{}^{3}S_1^{}}$ [fm]	$a_{{}_{}{}^{3}S_1^{}}$ [fm]
LO	$1.97$	$2.00\left(1\right)$	$2.041\left(1\right)$	$-22.4\left(4\right)$	$1.686\left(1\right)$	$5.46\left(1\right)$
	$1.64$	$2.07\left(1\right)$	$2.093\left(5\right)$	$-22.5\left(7\right)$	$1.6932\left(8\right)$	$5.45\left(1\right)$
	$1.32$	$2.12\left(2\right)$	$2.11\left(2\right)$	$-22.5\left(5\right)$	$1.71\left(1\right)$	$5.44\left(1\right)$
NLO	$1.97$	$2.2246{\left(3\right)}^{†}$	$2.4\left(6\right)$	$-23\left(4\right)$	$1.79\left(3\right)$	$5.31\left(2\right)$
	$1.64$	$2.2246{\left(1\right)}^{†}$	$2.3\left(1\right)$	$-23\left(2\right)$	$1.73\left(1\right)$	$5.33\left(1\right)$
	$1.32$	$2.2246{\left(1\right)}^{†}$	$2.47\left(3\right)$	$-23\left(1\right)$	$1.70\left(1\right)$	$5.336\left(9\right)$
NNLO	$1.97$	$2.2246{\left(3\right)}^{†}$	$2.6\left(6\right)$	$-24\left(4\right)$	$1.82\left(3\right)$	$5.35\left(2\right)$
	$1.64$	$2.2246{\left(1\right)}^{†}$	$2.5\left(3\right)$	$-23\left(2\right)$	$1.74\left(1\right)$	$5.36\left(1\right)$
	$1.32$	$2.22457{\left(7\right)}^{†}$	$2.6\left(2\right)$	$-23\left(1\right)$	$1.744\left(7\right)$	$5.382\left(5\right)$
experiment	$-$	$2.224575\left(9\right)$	$2.77\left(5\right)$	$-23.740\left(20\right)$	$1.753\left(8\right)$	$5.419\left(7\right)$

Table 7. Table 7: Summary of fit results (in units of a 𝑎 a ) for the perturbative NLO+NNLO analysis at a = 1.97 𝑎 1.97 a=1.97  fm. Fitted values of E d subscript 𝐸 𝑑 E_{d} are indicated by a dagger ( † † \dagger ). Note that the values of C S 0 1 subscript 𝐶 superscript subscript 𝑆 0 1 C_{{}^{1}S_{0}} , C S 1 3 subscript 𝐶 superscript subscript 𝑆 1 3 C_{{}^{3}S_{1}} and b s subscript 𝑏 𝑠 b_{s} are fixed by the LO fit.

	LO	NLO	NNLO
$C_{{}_{}{}^{1}S_0^{}}$	$-0.462\left(8\right)$	$-$	$-$
$C_{{}_{}{}^{3}S_1^{}}$	$-0.633\left(6\right)$	$-$	$-$
$b_s$	$0.054\left(3\right)$	$-$	$-$
$\Delta C$	$-$	$-0.2\left(2\right)$	$-0.0\left(2\right)$
$\Delta C_{I^2}$	$-$	$-0.02\left(9\right)$	$0.03\left(9\right)$
$C_{q^2}$	$-$	$0.03\left(5\right)$	$0.12\left(5\right)$
$C_{I^2,q^2}$	$-$	$0.04\left(2\right)$	$0.03\left(3\right)$
$C_{S^2,q^2}$	$-$	$-0.05\left(5\right)$	$-0.02\left(5\right)$
$C_{S^2,I^2,q^2}$	$-$	$0.00\left(2\right)$	$-0.01\left(2\right)$
$C_{{\left(q\cdot S\right)}^2}$	$-$	$0.06\left(2\right)$	$-0.05\left(2\right)$
$C_{I^2,{\left(q\cdot S\right)}^2}$	$-$	$-0.10\left(2\right)$	$-0.07\left(2\right)$
$C_{\left(q\times S\right)\cdot k}^{I=1}$	$-$	$0.039\left(5\right)$	$0.038\left(5\right)$
$E_d$ [MeV]	$-2.02\left(4\right)$	$-2.224{\left(3\right)}^{†}$	$-2.224{\left(3\right)}^{†}$

Table 8. Table 8: Summary of fit results (in units of a 𝑎 a ) for the perturbative NLO+NNLO analysis at a = 1.64 𝑎 1.64 a=1.64  fm. Notation as in Table  7 .

	LO	NLO	NNLO
$C_{{}_{}{}^{1}S_0^{}}$	$-0.47\left(1\right)$	$-$	$-$
$C_{{}_{}{}^{3}S_1^{}}$	$-0.64\left(1\right)$	$-$	$-$
$b_s$	$0.091\left(5\right)$	$-$	$-$
$\Delta C$	$-$	$-0.2\left(2\right)$	$0.3\left(3\right)$
$\Delta C_{I^2}$	$-$	$-0.00\left(9\right)$	$0.1\left(1\right)$
$C_{q^2}$	$-$	$0.04\left(6\right)$	$0.18\left(6\right)$
$C_{I^2,q^2}$	$-$	$0.08\left(3\right)$	$0.06\left(3\right)$
$C_{S^2,q^2}$	$-$	$-0.05\left(5\right)$	$0.00\left(6\right)$
$C_{S^2,I^2,q^2}$	$-$	$-0.01\left(3\right)$	$-0.01\left(3\right)$
$C_{{\left(q\cdot S\right)}^2}$	$-$	$0.06\left(3\right)$	$0.08\left(4\right)$
$C_{I^2,{\left(q\cdot S\right)}^2}$	$-$	$-0.07\left(3\right)$	$-0.06\left(4\right)$
$C_{\left(q\times S\right)\cdot k}^{I=1}$	$-$	$0.051\left(9\right)$	$0.05\left(1\right)$
$E_d$ [MeV]	$-2.13\left(4\right)$	$-2.224{\left(2\right)}^{†}$	$-2.224{\left(2\right)}^{†}$

Table 9. Table 9: Summary of fit results (in units of a 𝑎 a ) for the perturbative NLO+NNLO analysis at a = 1.32 𝑎 1.32 a=1.32  fm. Notation as in Table  7 .

	LO	NLO	NNLO
$C_{{}_{}{}^{1}S_0^{}}$	$-0.44\left(1\right)$	$-$	$-$
$C_{{}_{}{}^{3}S_1^{}}$	$-0.59\left(1\right)$	$-$	$-$
$b_s$	$0.18\left(1\right)$	$-$	$-$
$\Delta C$	$-$	$0.0\left(2\right)$	$0.4\left(2\right)$
$\Delta C_{I^2}$	$-$	$0.05\left(9\right)$	$0.30\left(9\right)$
$C_{q^2}$	$-$	$0.04\left(6\right)$	$0.62\left(6\right)$
$C_{I^2,q^2}$	$-$	$0.19\left(4\right)$	$0.07\left(3\right)$
$C_{S^2,q^2}$	$-$	$-0.03\left(5\right)$	$0.07\left(5\right)$
$C_{S^2,I^2,q^2}$	$-$	$-0.01\left(3\right)$	$-0.12\left(3\right)$
$C_{{\left(q\cdot S\right)}^2}$	$-$	$0.09\left(4\right)$	$0.02\left(4\right)$
$C_{I^2,{\left(q\cdot S\right)}^2}$	$-$	$-0.05\left(4\right)$	$0.12\left(4\right)$
$C_{\left(q\times S\right)\cdot k}^{I=1}$	$-$	$0.12\left(1\right)$	$0.11\left(1\right)$
$E_d$ [MeV]	$-2.14\left(3\right)$	$-2.224{\left(1\right)}^{†}$	$-2.224{\left(1\right)}^{†}$

Table 10. Table 10: Summary of fit results with perturbatively improved OPE (in units of a 𝑎 a ) for the perturbative NLO+NNLO analysis at a = 1.97 𝑎 1.97 a=1.97  fm. Notation as in Table  7 .

	LO	NLO	NNLO
$C_{{}_{}{}^{1}S_0^{}}$	$-0.462\left(8\right)$	$-$	$-$
$C_{{}_{}{}^{3}S_1^{}}$	$-0.633\left(6\right)$	$-$	$-$
$b_s$	$0.054\left(3\right)$	$-$	$-$
$\Delta C$	$-$	$-0.2\left(3\right)$	$-0.0\left(3\right)$
$\Delta C_{I^2}$	$-$	$-0.1\left(1\right)$	$0.03\left(9\right)$
$C_{q^2}$	$-$	$-0.03\left(7\right)$	$0.05\left(7\right)$
$C_{I^2,q^2}$	$-$	$0.09\left(3\right)$	$0.06\left(3\right)$
$C_{S^2,q^2}$	$-$	$-0.05\left(6\right)$	$0.00\left(6\right)$
$C_{S^2,I^2,q^2}$	$-$	$0.00\left(2\right)$	$-0.03\left(3\right)$
$C_{{\left(q\cdot S\right)}^2}$	$-$	$0.02\left(2\right)$	$-0.03\left(3\right)$
$C_{I^2,{\left(q\cdot S\right)}^2}$	$-$	$-0.07\left(2\right)$	$0.10\left(3\right)$
$C_{\left(q\times S\right)\cdot k}^{I=1}$	$-$	$0.014\left(7\right)$	$0.012\left(5\right)$
$E_d$ [MeV]	$-2.02\left(4\right)$	$-2.224{\left(3\right)}^{†}$	$-2.224{\left(3\right)}^{†}$

Table 11. Table 11: Nuclear binding energies with 2N forces up to NNLO in the NLEFT expansion for a = 1.97 𝑎 1.97 a=1.97  fm, data taken from Ref.  Lahde:2013uqa . The first parenthesis gives the estimated Monte Carlo error in the calculation of E B NNLO superscript subscript 𝐸 𝐵 NNLO E_{B}^{\mathrm{NNLO}} , and the second parenthesis the error due to variance-covariance matrix in Eq. ( 102 ). For reference, we also show the experimental binding energies.

	$E_B^{\mathrm{NNLO}}$(2N)	$E_B$(exp)
4He	$-25.60\left(6\right)\left(2\right)$	$-28.30$
8Be	$-48.6\left(1\right)\left(3\right)$	$-56.35$
12C	$-78.7\left(2\right)\left(5\right)$	$-92.16$
16O	$-121.4\left(5\right)\left(7\right)$	$-127.62$
20Ne	$-163.6\left(9\right)\left(9\right)$	$-160.64$
24Mg	$-208\left(2\right)\left(2\right)$	$-198.26$
28Si	$-275\left(3\right)\left(2\right)$	$-236.54$

Equations

$$M\equiv :\exp \left(-{\alpha }_{t}H\right):,$$

$$H\equiv H_{\mathrm{free}}+V_{\mathrm{LO}}+V_{\mathrm{NLO}}+\dots ,$$

$$E_{\lambda }=-\frac{1}{a_t}\log \lambda ,$$

$$|R⟩\equiv \sum _{{\vec{R}}^{\prime }}{Y}_{l,l_z}({\hat{R}}^{\prime }){\delta }_{R,R^{\prime }}|{\vec{R}}^{\prime }⟩,$$

$$M_{{\vec{R}}^{\prime },\vec{R}}\to {\tilde{M}}_{R^{\prime },R}.$$

$${\psi }_l(r)=A{j}_l(qr)+B{n}_l(qr),$$

$$H_{\mathrm{free}}$$

$${\omega }_k^{\mathrm{str}}\equiv {\omega }_k^{\mathcal{O}(a^4)}+\mathcal{N}\left({\omega }_k^{\mathcal{O}(a^4)}-{\omega }_k^{\mathcal{O}(a^2)}\right),$$

$${\mathcal{O}}_1^{(0)}\equiv \frac{1}{2}:\sum _{\vec{n}}\rho (\vec{n})\rho (\vec{n}):,$$

$${\mathcal{O}}_2^{(0)}\equiv \frac{1}{2}:\sum _{\vec{n}}\sum _I{\rho }_I(\vec{n}){\rho }_I(\vec{n}):,$$

$$\Bigg{[}\begin{array}[]{c}C\vspace{.1cm}\\ C_{I}\\ \end{array}\Bigg{]}=\frac{1}{4}\Bigg{[}\begin{array}[]{rr}3&1\vspace{.1cm}\\ 1&-1\\ \end{array}\Bigg{]}\Bigg{[}\begin{array}[]{c}C_{{}^{1}S_{0}}\vspace{.1cm}\\ C_{{}^{3}S_{1}}\\ \end{array}\Bigg{]},$$

$${\mathcal{O}}_1^{(0)}\to \frac{1}{2L^3}:\sum _{\vec{q}}f(\vec{q})\rho (\vec{q})\rho (-\vec{q}):,$$

$${\mathcal{O}}_2^{(0)}\to \frac{1}{2L^3}:\sum _{\vec{q}}f(\vec{q}){\rho }_I(\vec{q}){\rho }_I(-\vec{q}):,$$

$$f(\vec{q})\equiv f_0^{-1}\exp \left(-b_s\frac{\vec{q}{}^4}{4}\right),$$

$$f_0\equiv \frac{1}{L^3}\sum _{\vec{q}}\exp \left(-b_s\frac{\vec{q}{}^4}{4}\right),$$

$$\frac{\vec{q}{}^2}{2}$$

$$P^{(0,1)}$$

$$P^{(1,0)}$$

$${\mathcal{O}}_{(0,1)}^{(0)}$$

$${\mathcal{O}}_{(1,0)}^{(0)}$$

$${\mathcal{O}}_1^{(2)}$$

$${\mathcal{O}}_2^{(2)}$$

$${\mathcal{O}}_3^{(2)}$$

$${\mathcal{O}}_4^{(2)}$$

$${\mathcal{O}}_5^{(2)}$$

$${\mathcal{O}}_6^{(2)}$$

$${\mathcal{O}}_7^{(2)}$$

$$\nabla_{l}f(\vec{n})\equiv\frac{1}{2a}\big{[}f(\vec{n}+a\hat{e}_{l})-f(\vec{n}-a\hat{e}_{l})\big{]},$$

$$\nabla_{l}^{2}f(\vec{n})\equiv\nabla_{l}\big{[}\nabla_{l}f(\vec{n})\big{]}.$$

$${⟨}^1{S}_0|V_{\mathrm{con}}^{(2)}{|}^1{S}_0⟩$$

$${⟨}^1{P}_1|V_{\mathrm{con}}^{(2)}{|}^1{P}_1⟩$$

$${⟨}^3{P}_0|V_{\mathrm{con}}^{(2)}{|}^3{P}_0⟩$$

$${⟨}^3{P}_1|V_{\mathrm{con}}^{(2)}{|}^3{P}_1⟩$$

$$⟨{}^{3}S{D}_1|V_{\mathrm{con}}^{(2)}|{}^{3}S{D}_1⟩$$

Full text

I Introduction

Nuclear Lattice Effective Field Theory (NLEFT) has recently gained prominence as an ab initio method for the study of nuclear structure formation at low energies. The advent of NLEFT has largely been due to rapid developments in computational algorithms and resources, which have enabled the efficient combination of lattice Monte Carlo methods with the low-energy effective field theory of QCD, known as Chiral Perturbation Theory or Chiral Effective Field Theory. Such progress has greatly increased our ability to exploit the advantages of the EFT method in the realm of many-body nuclear physics, which remains a highly challenging area of study. Hence, impressive progress has been made within NLEFT in furthering our understanding of the spectra, structure and scattering of light- and medium-mass nuclei Epelbaum:2011md ; Epelbaum:2012qn ; Epelbaum:2012iu ; Epelbaum:2013paa ; Bour:2014bxa ; Elhatisari:2015iga ; Elhatisari:2016owd , see also Ref. Lee:2008fa for an early review.

Chiral EFT provides a model-independent approach to hadronic interactions at the energy scales of interest for nuclear physics. Based on the spontaneous and explicit chiral symmetry breaking of QCD, Chiral EFT provides a systematic treatment of such interactions in terms of a generic soft scale ($Q$) which is commonly taken to refer to the Goldstone boson mass (such as the pion mass $M_{\pi}$) or to external nucleon momenta. In the nuclear physics context, the EFT is used to work out the interaction potential between the nuclear constituents. These chiral potentials are then used in an appropriate framework to generate the bound and scattering states. For the case of the nucleon-nucleon (NN) interaction considered here, Chiral EFT also clarifies the observed hierarchy between many-body contributions to the nuclear force. This power counting can be expressed in terms of $Q/\Lambda_{\chi}$, where $\Lambda_{\chi}$ refers to the hard scale at which chiral symmetry is restored Weinberg:1990rz . The contributions to the NN force are then classified as leading order (LO) for $(Q/\Lambda_{\chi})^{0}$, followed by next-to-leading order (NLO) for $(Q/\Lambda_{\chi})^{2}$, and next-to-next-to-leading order (NNLO) for $(Q/\Lambda_{\chi})^{3}$ etc., in decreasing order of importance. For a recent review of Chiral EFT in nuclear physics, see Ref. Epelbaum:2008ga . It should also be noted that Chiral EFT provides a method to systematically estimate the uncertainty of a calculation at a given order in the EFT expansion, which is of special relevance is searches of physics beyond the Standard Model (BSM). With the advent of precision experiments searching for BSM physics, the importance of well-controlled error estimates for the nuclear contributions have become essential for the statistical interpretation of purported BSM signals and hence, ultimately, for any claim of detection, see e.g. Ref. Hoferichter:2016nvd (and references therein).

The fundamental problem of neutron-proton scattering in NLEFT was first studied at LO in Ref. Borasoy:2006qn , and later extended to NLO in Ref. Borasoy:2007vi , with phase shifts and mixing angles calculated on the lattice using the so-called spherical wall method Borasoy:2007vy . This method was used earlier in the context of variational calculations of resonant states in 4He Carlson:1984zz . We shall here revisit, in a systematical manner, the calculation of neutron-proton scattering observables, which also serve to determine the low-energy constants (LECs) of the NLO contact terms of NLEFT. Our work is based on an improvement of the spherical wall method known as the radial Hamiltonian formalism, which was proposed and pioneered in Ref. Elhatisari:2015iga in the study of alpha-alpha scattering on the lattice. In this formalism, the two-nucleon problem is formulated in terms of radial coordinates for each partial wave. Specifically, lattice points with the same radial coordinate are grouped together and weighted by the appropriate spherical harmonics, which eliminates the need to work with a computationally costly $L^{3}\times L^{3}$ problem (where $L$ is the linear dimension of the cubic lattice) without loss of precision. This approach can be further accelerated by binning lattice points with similar radial coordinates into segments of width $a_{R}$, as proposed in Ref. Elhatisari:2016hby . The advantages of the radial Hamiltonian method were already demonstrated in Ref. Lu:2015riz for the phase shifts and mixing angles of a system of two nucleons with a simplified model potential.

In the present work, we address the task of determining the LECs of the two-nucleon force at NLO and NNLO in NLEFT, by means of a chi-square minimization with respect to neutron-proton phase shifts and mixing angles. This procedure also allows us to provide quantitative estimates of the uncertainties of the NLO constants in NLEFT, along with estimates of their systematical errors and the impact of such errors on the binding energies of nuclei. It should be noted that the pioneering calculations of Refs. Borasoy:2006qn ; Borasoy:2007vi (and almost all calcuations of nuclear properties) were performed with a coarse lattice spacing of $a=1.97$ fm, which corresponds to a relatively low momentum cutoff of $\pi/a=314$ MeV ⋆⋆\star⋆⋆\star$\star$Note that such soft nucleon-nucleon interactions lead to better convergence properties in the calculations of many-nucleon systems and nuclear matter, see e.g. Bogner:2005sn .. Here, we now also study the effects of decreasing the lattice spacing to $a\simeq 1$ fm, which greatly decreases the impact of lattice artifacts and systematical errors, and discuss the possibility of further improving the lattice action to decrease remaining discretization effects. Note that a first study of discretization errors and lattice spacing variation at LO has been performed in Ref. Klein:2015vna . Finally, our study of lattice spacing variation requires that the two-pion exchange potential (TPEP) is explicitly accounted for. In prior work at $a=1.97$ fm, the TPEP at NLO and NNLO contributions were integrated out by means of a Taylor expansion in powers of $q^{2}/(4M_{\pi})$. Since we now use lattice spacings as small as $a\simeq 1$ fm, we need to include the full structure of the TPEP in our analysis.

Our paper is organized as follows: The lattice EFT formalism, the radial Hamiltonian method, and the NLEFT potentials up to NNLO are presented in Section II. In Section III, we give the results of a fully non-perturbative calculation of coupled-channel neutron-proton scattering up to NNLO, followed by a treatment where the NLO and NNLO contributions are computed perturbatively. We also study the lattice spacing dependence of the calculated phase shifts and mixing angles. Furthermore, we investigate how the uncertainty in the four-nucleon LECs propagates into the prediction of nuclear ground-state energies. In Section IV, we conclude with a brief discussion of planned N3LO calculations and other future directions.

II Lattice formalism

We begin with a detailed description of the NLEFT lattice Hamiltonian on which our calculations are based. We denote the (spatial) lattice spacing by $a$, the temporal lattice spacing by $a_{t}$, and we also define $\alpha_{t}\equiv a_{t}/a$. Our lattice is a periodic cube of volume $L^{3}$. For non-zero temporal lattice spacing, we define the transfer matrix as Borasoy:2006qn

[TABLE]

with the Hamiltonian

[TABLE]

where $H_{\mathrm{free}}$ is the free nucleon Hamiltonian and $V_{\mathrm{LO}}$, $V_{\mathrm{NLO}}$, etc. contain nucleon-nucleon interactions of progressively higher order in NLEFT. The colons in Eq. (1) denote normal ordering. The energy eigenvalues are given by

[TABLE]

where $\lambda$ denotes an eigenvalue of $M$.

Following Ref. Elhatisari:2015iga , we construct the transfer matrix in radial coordinates. Specifically, we group the lattice points $(n_{x},n_{y},n_{z})$ with the same radial coordinate, by weighting them with the spherical harmonics. Thus, instead of working with the full basis $|\vec{R}\rangle$, one obtains the reduced basis

[TABLE]

where $Y_{l,l_{z}}$ is the spherical harmonic for angular momentum quantum numbers $(l,l_{z})$ and $\delta$ denotes the Kronecker delta. We thus obtain the radial transfer matrix

[TABLE]

A similar approach with a refined grid for the radial lattice was performed in Elhatisari:2016hby .

To determine phase shifts and mixing angles, we apply the method proposed in Ref. Lu:2015riz , whereby these are extracted directly from the radial wave functions. Specifically, one defines three radii, $R_{\mathrm{in}}$, $R_{\mathrm{out}}$ and $R_{\mathrm{wall}}$. The NN interaction contributes in the range $0<r<R_{\mathrm{in}}$, while an infinite spherical wall barrier is applied for $r>R_{\mathrm{wall}}$. In the range $R_{\mathrm{in}}<r<R_{\mathrm{out}}$, the NN interaction vanishes and the wave function can be expanded as a linear combination of the spherical Bessel and Neumann functions, according do

[TABLE]

from which phase shifts and mixing angles can be extracted. For more details, see Ref. Lu:2015riz and the earlier work of Ref. Borasoy:2007vy . Typical values used later are $R_{\rm in}\simeq 24$ fm, $R_{\rm out}\simeq 28$ fm, and $R_{\rm wall}\simeq 30$ fm.

We shall now give a detailed description of the Hamiltonian $H$, and its various contributions. The free nucleon Hamiltonian is given by Borasoy:2006qn

[TABLE]

where the $\hat{e}_{l}$ with $l=1,2,3$ are unit vectors in the spatial directions, and $m_{N}$ is the nucleon mass. In Table 1, we give the hopping coefficients $\omega_{k}$ for lattice actions up to $\mathcal{O}(a^{4})$. Throughout our work, we use the so-called stretched action which is defined in terms of the $\mathcal{O}(a^{4})$- and $\mathcal{O}(a^{2})$-improved actions Lee:2008fa . This gives the stretched hopping coefficients

[TABLE]

where $\mathcal{N}=10$ is adopted in the present calculations.

In Chiral EFT, the NN force is decomposed into the long-range components arising from the exchange of pions, and short-range contributions described by contact interactions with increasing powers of momenta. Such two-nucleon contact operators introduce unknown coefficients which we determine by fitting the data on neutron-proton phase shifts and mixing angles. In what follows, we present our contact and pion exchange operators.

II.1 Contact interactions

We begin our treatment of the lattice Chiral EFT interaction by considering the various contact operators that appear up to NNLO in the chiral expansion. At LO, we consider the following operators

[TABLE]

and

[TABLE]

as the independent contact operators, with coefficients $C$ and $C_{I}$, respectively. Here, $\rho(\vec{n})$ and $\rho_{I}(\vec{n})$ are the local density and local isospin density operators on the lattice, which are defined in App. A. At LO, the coefficients $C$ and $C_{I}$ are determined by the spin-singlet (${}^{1}S_{0}$) and the spin-triplet (${}^{3}S_{1}$) $S$-wave channels, and can be parameterized as

[TABLE]

where $C_{{}^{1}S_{0}}$ and $C_{{}^{3}S_{1}}$ are determined by fitting scattering data in the ${}^{1}S_{0}$ and ${}^{3}S_{1}$ channels.

In Ref. Borasoy:2006qn , it was shown that an on-site interaction such as those shown in Eqs. (9) and (10) do not suffice to provide a favorable description of the $S$-wave phase shifts except at very low momenta. Hence, smeared contact operators were introduced according to

[TABLE]

and

[TABLE]

where the smearing factor $f(\vec{q}\,)$ is

[TABLE]

with $b_{s}$ a free parameter, and the normalization is given by

[TABLE]

with

[TABLE]

where the $q_{l}$ are lattice momentum components, and the $\mathcal{O}(a^{4})$-improved hopping coefficients $\omega_{i}$ are given in Table 1.

In the analysis of the Ref. Borasoy:2007vk , smeared contact operators were found to dramatically improve the convergence of the NLEFT expansion in the $S$-wave channels, at the price of introducing unwanted attractive forces in the $P$-wave channels. By means of the projection operators Epelbaum:2008vj ,

[TABLE]

for the $(S,I)=(0,1)$ and $(1,0)$ channels, good agreement at LO in the $P$-wave channels can be recovered (although a similar problem of unwanted forces in the $D$-wave channels persists). In the present work, we use the corresponding smeared LO contact operators

[TABLE]

for $(S,I)=(0,1)$, and

[TABLE]

for $(S,I)=(1,0)$, where $\rho_{S}(\vec{n})$ and $\rho_{S,I}(\vec{n})$ are local spin density and local spin-isospin density operators, defined in App. A.

According to chiral EFT power counting, there are seven independent contact operators with two derivatives at NLO. Here, we use the basis and lattice formulation of Ref. Borasoy:2007vi , which leads to the following NLO contact operators

[TABLE]

where $\Pi_{l}(\vec{n})$ and $\Pi_{l,S}(\vec{n})$ denote current density and spin-current density operators, the lattice definitions of which are given in App. A. Following the treatment of Ref. Borasoy:2007vi for the spin-orbit operator $\mathcal{O}_{7}^{(2)}$, we project onto $I=1$, giving

[TABLE]

which eliminates lattice artifacts in the $S=1$ even-parity channels. For the derivative operator $\nabla_{l}$ in the NLO contact terms, we use

[TABLE]

where $a$ is the spatial lattice spacing, and $\hat{e}_{l}$ is a unit vector in spatial direction $l$. For the double derivative operator $\nabla_{l}^{2}$, we take

[TABLE]

In the radial transfer matrix formalism, we project each of the NLO contact operators onto the NN partial waves under consideration, such that $V_{X}^{i}$ is the matrix element of operator $i$ in channel $X$. If we denote the complete set of NLO contact interactions by $V_{\mathrm{con}}^{(2)}$, we find

[TABLE]

for the uncoupled channels, and

[TABLE]

for the coupled ones. It is clear that only certain combinations of the contact operators contribute to each partial wave, which allows for a simplified fitting procedure. Specifically, we determine $C_{{}^{1}S_{0}}$ and $\widetilde{C}_{1}$ by fitting the ${}^{1}S_{0}$ channel, $\widetilde{C}_{4}$ by means of the ${}^{1}P_{1}$ channel, $\widetilde{C}_{5}$, $\widetilde{C}_{6}$ and $\widetilde{C}_{7}$ from a simultaneous fit to the ${}^{3}P_{0}$, ${}^{3}P_{1}$ and ${}^{3}P_{2}$-${}^{3}F_{2}$ channels, and finally $C_{{}^{3}S_{1}}$, $\widetilde{C}_{2}$ and $\widetilde{C}_{3}$ by fitting the ${}^{3}S_{1}$-${}^{3}D_{1}$ channel.

We note that the fitted LECs $\widetilde{C}_{i}$ are given in terms of those of the NLO operators in Eqs. (21) through (26) and (28) by the relation

[TABLE]

which can be inverted in order to find the original LECs $C_{i}$, once the $\widetilde{C}_{i}$ have been determined.

II.2 Long-range interactions

Next, we consider the long-range one-pion exchange (OPE) and two-pion exchange (TPE) contributions to the chiral EFT interaction. The latter contributes at NLO and NNLO. At LO, the OPE potential is given by Borasoy:2006qn ; Lee:2008fa

[TABLE]

where the pion propagator is

[TABLE]

with

[TABLE]

where the $k_{i}$ are lattice momentum components. For the denominator of Eq. (61), we take

[TABLE]

using the $\mathcal{O}(a^{4})$-improved hopping coefficients $\omega_{i}$ of Table 1. For the numerator of Eq. (61), we take

[TABLE]

which coincides with the choice of derivative operator in Eq. (29). We also include the isospin-breaking (IB) effects due to the pion mass differences. Specifically, we take

[TABLE]

for the isospin-triplet and isospin-singlet channels, respectively. This approach is consistent with the conventions of the Nijmegen partial wave analysis. For more details on the IB corrections to the NN interaction, see Refs. Epelbaum:2004fk ; Epelbaum:2008ga ; Machleidt:2011zz (and references therein).

The first contribution from the TPE potential appears at NLO in chiral EFT. We note that several prior continuum calculations including TPE exist. For instance, in Refs. Friar:1994zz ; Kaiser:1997mw , dimensional regularization (DR) was used to remove the divergence appearing in the loop integral, and a non-local momentum-dependent form factor was applied to suppress the high-momentum contributions when solving the Lippmann-Schwinger equation. In Ref. Epelbaum:2003gr , another regularization called spectral function regularization (SFR) was proposed. Compared to DR, the SFR method introduces an additional cutoff to remove the short-range components of the TPE potential. Recently, a new position-space regularization was proposed in Refs. Gezerlis:2013ipa ; Gezerlis:2014zia ; Epelbaum:2014efa . The study of effects in nuclear lattice EFT due to different choices of regularization of the TPE is beyond the scope of the current work. In this work, we use the DR expressions with discretized lattice momenta. We also note that the lattice spacing serves as a natural UV cut-off.

Thus far, nuclear lattice EFT calculations have been performed with a lattice spacing of $a=1.97$ fm, and hence the TPE potentials at NLO and NNLO have not been included explicitly, but rather been absorbed into the contact terms. Since we are here studying the effects of reducing the lattice spacing to $a\simeq 1$ fm, we shall for the first time include the full TPE structure. As for the smeared LO contact terms and the OPE potential, we define the lattice formulation of the TPE potential in momentum space, and Fourier transform the results to coordinate space. The TPE potential is of the form

[TABLE]

at NLO. The explicit expressions for the components of Eq. (66) are

[TABLE]

with

[TABLE]

and

[TABLE]

with

[TABLE]

and

[TABLE]

with

[TABLE]

where the function $F^{(2)}(\vec{q}\,)$ is given by

[TABLE]

and $L(\vec{q}\,)$ is the loop function

[TABLE]

in DR. In order to coincide with the definitions of the derivative operator (29) and the double-derivative operator (30), we take

[TABLE]

and

[TABLE]

which ensures that the divergences appearing in the loop diagrams can be absorbed by tuning the contact interaction LECs $C_{i}$.

Similarly, we parameterize the sub-leading (NNLO) contribution to the TPE as

[TABLE]

where

[TABLE]

with

[TABLE]

and

[TABLE]

with

[TABLE]

and

[TABLE]

with

[TABLE]

where the function $F^{(3)}(\vec{q}\,)$ is given by

[TABLE]

and $A(\vec{q}\,)$ is the loop function

[TABLE]

in DR. For the momenta $\vec{q}$, we again apply the conventions of Eqs. (63) and (76).

III Results

We now turn to a description of our calculational methods. We take $F_{\pi}=92.2$ MeV for the pion decay constant, and $g_{A}=1.29$ for the nucleon axial coupling constant to account for the Goldberger-Treiman discrepancy Epelbaum:2004fk . For the nucleon mass, we use $m_{N}=938.38$ MeV, and for the charged and neutral pion masses, we take $M_{\pi^{\pm}}=139.75$ MeV and $M_{\pi^{0}}=134.98$ MeV, respectively. We use the isospin-averaged pion mass

[TABLE]

in the TPEP expressions at NLO and NNLO. For the constants $c_{1}$, $c_{3}$ and $c_{4}$ that appear in the TPEP at NNLO, we use $c_{1}=-1.10(3)$ GeV*-1*, $c_{3}=-5.54(6)$ GeV*-1* and $c_{4}=4.17(4)$ GeV*-1* from the accurate Roy-Steiner analysis of pion-nucleon scattering adopted to the counting of the nucleon mass used here Hoferichter:2015tha . Also, as the uncertainties of these LECs are very small, we only consider the central values in the following.

We determine the optimal parameter values for the NLEFT action up to NNLO by performing a chi-square fit to neutron-proton phase shifts and mixing angles. For this purpose, we define the uncertainties of the empirical scattering observables (in each partial wave) according to Refs. Epelbaum:2014efa ; Epelbaum:2014sza , which gives

[TABLE]

where $\Delta_{i}^{\mathrm{PWA}}$ denotes the uncertainty of the PWA, while $\delta_{i}^{\mathrm{PWA}}$ signifies the phase shift (or mixing angle) in channel $i$ of the PWA (see also Ref. Stoks:1993tb ). Furthermore, $\delta_{i}^{\mathrm{NijmI}}$, $\delta_{i}^{\mathrm{NijmII}}$ and $\delta_{i}^{\mathrm{Reid93}}$ refer to the PWA results based on the Nijmegen I, Nijmegen II and Reid93 NN potentials, respectively. Hence, a measure of systematical error in the PWA is accounted for in our analysis. The $\chi^{2}$ function to be minimized is defined as

[TABLE]

where $i$ runs over all values of $p_{\mathrm{CM}}$ and channels included in the analysis. In Eq. (88), $\delta_{i}^{\mathrm{PWA}}$ is the phase shift (or mixing angle) at a given momentum $p_{\mathrm{CM}}$ from the Nijmegen PWA, $\delta_{i}^{\mathrm{cal}}$ is the corresponding calculated NLEFT value, and $\Delta_{i}$ is given by Eq. (III).

When fitting the phase shifts and mixing angles of the Nijmegen partial wave analysis, we note certain simplifying features. Specifically, at LO we determine $C_{{}^{1}S_{0}}$, $C_{{}^{3}S_{1}}$, and the smearing parameter $b_{s}$, by fitting the ${}^{1}S_{0}$ and ${}^{3}S_{1}$ phase shifts. At NLO and NNLO, we no longer update the value of $b_{s}$. At NLO, we determine $C_{{}^{1}S_{0}}$ and $\widetilde{C}_{1}$ by fitting the ${}^{1}S_{0}$ phase shift, $C_{{}^{3}S_{1}}$, $\widetilde{C}_{2}$ and $\widetilde{C}_{3}$ by fitting the ${}^{3}S_{1}$ phase shift and the mixing angle $\epsilon_{1}$, $\widetilde{C}_{4}$ by fitting the ${}^{1}P_{1}$ phase shift, and finally $\widetilde{C}_{5}$, $\widetilde{C}_{6}$ and $\widetilde{C}_{7}$ by fitting the the ${}^{3}P_{0}$, ${}^{3}P_{1}$ and ${}^{3}P_{2}$ phase shifts. The NNLO fits are similar, apart from the inclusion of the NNLO TPEP operators. We do not take the deuteron binding energy $E_{d}$ as an additional constraint in the LO fits, as we do not expect $E_{d}$ to be accurately reproduced in an LO calculation. At NLO and NNLO, the experimental value $E_{d}=2.224575(9)$ MeV is taken as an additional constraint. At LO, we fit up to center-of-mass momenta of $p_{\mathrm{CM}}^{\mathrm{max}}=100$ MeV, while at NLO and NNLO we fit up to $p_{\mathrm{CM}}^{\mathrm{max}}=150$ MeV. Our fitting procedure at each order in NLEFT is summarized in Table 3.

III.1 Phase shifts and mixing angles to NNLO

Prior NLEFT work has used a relatively coarse lattice spacing of $a=1.97$ fm, which corresponds to a momentum cutoff $\Lambda\sim\pi/a=314$ MeV. This relatively low cutoff may induce significant lattice artifacts, particularly at high momenta. With this in mind, we here aim to study the NN scattering problem for $a=(200~{}\mathrm{MeV})^{-1}=0.98$ fm, with a temporal lattice spacing of $a_{t}=(600~{}\mathrm{MeV})^{-1}$. The number of lattice points in each spatial dimension is $L=64$, thus the physical volume is $V=(La)^{3}\simeq(63~{}\mathrm{fm})^{3}$, which is expected to be large enough to accommodate the NN system without introducing significant finite volume effects for the energy region $p_{\mathrm{CM}}<200$ MeV studied here. Our lattice parameters are summarized in Table 2.

First, we consider the problem of neutron-proton scattering by treating all orders in NLEFT up to NNLO non-perturbatively, similar to what is done in the continuum. This means that we construct the transfer matrix according to

[TABLE]

where the potential terms are given by

[TABLE]

at LO,

[TABLE]

at NLO, and

[TABLE]

at NNLO. Our results for the smallest lattice spacing, $a=0.98$ fm, are shown in Fig. 1. Clearly, the description of the $S$-wave channels is quite good even at LO, particularly for ${}^{3}S_{1}$. Compared to LO, significant improvements occur at NLO and NNLO, in particular for the ${}^{1}P_{1}$, ${}^{3}P_{0}$ and ${}^{3}P_{2}$ channels, as well as for the mixing angle $\epsilon_{1}$. While the NLO contributions appear central for a good description of the $P$-waves and $\epsilon_{1}$, the TPE contributions at NNLO do not appear to produce a significant systematical effect, although we note that certain channels (such as ${}^{3}P_{2}$) show marked improvement at NNLO. While the results for the $D$-waves appear rather accurate, we note that the current way of smearing the LO contact interactions does produce unwanted additional forces in the $D$-wave channels, which should be dominated by OPE alone. The $D$-wave channels are expected to improve further upon addition of the N3LO contributions, which will be included in future work Du2017 .

In Table 3, we also give the value of $\chi^{2}/N_{\mathrm{dof}}$ for each of our fits ($a=0.98$ fm), where $N_{\mathrm{dof}}$ equals the number of fitted data points (phase shifts or mixing angles at a given momentum) minus the number of adjustable parameters. At LO with $a=0.98$ fm, we find $\chi^{2}/N_{\mathrm{dof}}\simeq 30$, which is reasonable given the rather stringent uncertainty criterion (III) of the PWA. This indicates that we have a satisfactory description of the ${}^{1}S_{0}$ and ${}^{3}S_{1}$ channels in the range $p_{\mathrm{CM}}<100$ MeV. At NLO, the main contribution to $\chi^{2}/N_{\mathrm{dof}}$ arises from $\epsilon_{1}$ with $p_{\mathrm{CM}}>100$ MeV, while at NNLO $\epsilon_{1}$ and the $P$-wave channels contribute roughly equally. These observations are consistent with the results shown in Fig. 1.

We also give the $S$-wave low-energy parameters for $a=0.98$ fm in Table 4, along with a summary of the fitted parameters. We find that the NLO and NNLO results clearly provide the closest agreement with the empirical scattering lengths and effective ranges, taken from Ref. Machleidt:2000ge . We note that $a_{{}^{3}S_{1}}$ and $r_{{}^{3}S_{1}}$ are both stable at various orders in NLEFT, and reasonably close to the empirical values. This is easily understood since the phase shift in the ${}^{3}S_{1}$ channel is accurately reproduced already at LO. For $a_{{}^{1}S_{0}}$ and $r_{{}^{1}S_{0}}$, a clear improvement is observed at NLO and NNLO compared to the results at LO. We also find that at NLO and NNLO, $E_{d}$ can be accommodated without sacrificing any accuracy in the other low-energy parameters. Finally, $C_{{}^{1}S_{0}}$ and $C_{{}^{3}S_{1}}$ for $a=0.98$ fm are in reasonably close agreement with the continuum results of Ref. Epelbaum:2014efa for a cutoff of $R=1.0$ fm, which suggests that lattice artifacts are under control.

III.2 Variation of the lattice spacing

Up to this point, we have mostly elaborated on our results for $a=0.98$ fm, which is the smallest lattice spacing we have considered. We shall next comment on our findings when the lattice spacing is varied in the range $1.97\geq a\geq 0.98$ fm, while the physical lattice volume is kept constant at $V=(La)^{3}\simeq(63~{}\mathrm{fm})^{3}$ (see Table 2 for a summary of lattice parameters). As we work within the transfer matrix formalism, the temporal lattice spacing $a_{t}$ should also be varied when $a$ is changed. Here, we choose $a_{t}$ such that $a_{t}/a^{2}$ is kept constant. This is motivated by the fact that the Hamiltonian scales with the lattice spacing as $H\sim 1/a^{2}$. For a pioneering LO calculation of the effects of varying $a$, see also Ref. Klein:2015vna .

In Table 5, we summarize the fitted constants of the NN interaction as a function of $a$, along with the $S$-wave low-energy parameters in Table 6. We note that the uncertainties of the fitted constants are obtained by an analysis of the variance-covariance matrix according to Eq. (115), while those of the $S$-wave parameters are obtained using Eq. (118). Our computed $S$-wave parameters appear very stable with respect to lattice spacing variation, which suggests that lattice spacing effects are small in the $S$-wave channels.

Our results for neutron-proton phase shifts and mixing angles for $a=1.32$ fm are shown in Fig. 2, for $a=1.64$ fm in Fig. 3, and finally for $a=1.97$ fm in Fig. 4. Together with the results for $a=0.98$ fm shown in Fig. 1, it is immediately apparent that lattice spacing effects are small for the $S$-waves in the range $0<p_{\mathrm{CM}}<200$ MeV, which is consistent with the behavior of the $S$-wave parameters. On the other hand, this situation is quite different for the $P$-waves and $D$-waves. For these higher partial waves, as well as for the mixing angles $\epsilon_{1}$ and $\epsilon_{2}$, the lattice spacing effects remain small only up to $p_{\mathrm{CM}}<100$ MeV. For $p_{\mathrm{CM}}>100$ MeV, the deviations from the Nijmegen PWA increase rapidly, but are nevertheless systematically reduced when $a$ is decreased.

To conclude, for the $S$-waves the lattice spacing effects remain small throughout the range of $p_{\mathrm{CM}}$ considered here, even for the (rather coarse) lattice spacing of $a=1.97$ fm. For the $P$-waves and $D$-waves, this situation holds only up to $p_{\mathrm{CM}}\simeq 100$ MeV. However, we note that $a=0.98$ fm suffices to give an accurate description for $p_{\mathrm{CM}}\simeq 200$ MeV, regardless of the channel under consideration. This suggests that the observed discrepancies could be eliminated by a combination of improved lattice momentum operators and N3LO effects, possibly taken together with a lattice spacing somewhat smaller than $a=1.97$ fm. We would like to stress that the phase shifts agree within uncertainties below 150 MeV (with a few exceptions) for the lattice spacings considered. This validates the statements made in Ref. Klein:2015vna about the lattice spacing independence of observables in the two-nucleon sector.

III.3 Perturbative treatment of higher orders

We have thus far demonstrated that non-perturbative fits to neutron-proton scattering data are feasible to any given order in NLEFT, provided that the requisite potential operators have been worked out. Nevertheless, for practical reasons (such as sign oscillations and increased computational complexity) the contributions of NLO and higher orders are usually treated perturbatively in Monte Carlo simulations of nuclear many-body systems. With this in mind, we show here how our analysis of phase shifts and mixing angles can be applied in a way consistent with current lattice Monte Carlo work.

Before discussing our results, we briefly summarize the differences between the perturbative and non-perturbative analyses. We again start with a LO fit, the parameters of which are fixed by fitting the ${}^{1}S_{0}$ and ${}^{3}S_{1}$ channels (but not $E_{d}$). As in the non-perturbative analysis, for the LO fits we consider data up to $p_{\mathrm{CM}}^{\mathrm{max}}=100$ MeV. For higher-order (NLO and NNLO) fits, we include data up to $p_{\mathrm{CM}}^{\mathrm{max}}=150$ MeV. Since higher orders in NLEFT are treated perturbatively, the transfer matrix is constructed in a different way than in Eq. (89). To be specific, in the perturbative analysis the transfer matrix is

[TABLE]

where

[TABLE]

and as in previous Monte Carlo studies of NLEFT, we introduce the additional operators Borasoy:2007vi

[TABLE]

which we classify as NLO perturbations and add to the NLO potential in Eq. (91) when $M^{\mathrm{pert}}$ is computed. This is done because the LO LECs are kept fixed and thus fitting these finite shifts is equivalent to a refit of the LO LECs, as it is done in the non-perturbative case. Additionally, $\Delta V$ and $\Delta V_{I^{2}}$ absorb part of the (sizable) short-distance contributions from TPE at NLO and NNLO. At NLO, we also studied an operator of the form $\sum_{i}\tau_{1,i}\tau_{2,i}q_{i}^{2}$ which accounts for rotational symmetry breaking effects on the lattice, but no significant effects were observed.

As for the non-perturbative case, we give results for a range of lattice spacings for the perturbative analysis. The fitted parameters for $a=1.97$ fm, $a=1.64$ fm and $a=1.32$ fm are given in Tables 7, 8 and 9, respectively. The corresponding phase shifts and mixing angles are shown in Figs. 5, 6 and 7. For each computed phase shift, we provide an estimated uncertainty according to

[TABLE]

where $\mathcal{E}_{ij}$ denotes the variance-covariance matrix of the fitted parameters, according to Eq. (114), and $J_{\delta}$ is the Jacobian vector of the phase shift (or mixing angle) in question. The last factor in Eq. (97) is the so-called Birge factor described in App. B, which approximately accounts for the systematical errors in the analysis.

At LO, we reproduce well the low-momentum region, and obtain a realistic deuteron binding energy. In particular, we note that the ${}^{3}S_{1}$ PWA data are almost perfectly reproduced. This is largely caused by the very accurate PWA data of this channel, which gives this channel a relatively high weight in the $\chi^{2}$ function. We note that this may potentially worsen the agreement in other channels, where a comparable accuracy of the PWA data is not available. Also, the expectation is that the $P$-waves should be well described at LO, since they are dominated by the OPEP contribution. The reason why this is not the case for our LO results is that, in the perturbative calculation, and in order to be consistent with the Monte Carlo simulations, we treat the momentum $\vec{q}\,{}^{2}$ in the denominator of the OPE as in Eq. (62), and factors of $\vec{q}$ as in Eq. (63). This choice considerably suppresses the OPEP contribution already at intermediate momenta, which worsens the description of the $P$-waves.

Moving to NLO, a significant improvement is found in some channels, particularly for ${}^{1}P_{1}$ and ${}^{3}P_{1}$, where the PWA is now well described up to $\sim 100$ MeV. On the other hand, we note that the ${}^{3}P_{0}$ and ${}^{3}P_{2}$ channels, as well as the $D$-waves, show little improvement. We attribute these features to the deficiencies in the OPE as mentioned above. The description of the $S$-wave channels is found to improve at intermediate momenta, which is mainly due to the NLO contact terms and to the parts of the NLO TPEP that contribute to the $S$-waves.

At NNLO, while no new unknown parameters contribute, the sub-leading TPEP enters as a prediction from $\pi N$ scattering in Chiral EFT. Thus, the NLO constants are refitted at NNLO in order to absorb the strong short-distance isoscalar contributions from the $\pi N$ LECs. The NLO and NNLO results appear in most cases virtually indistinguishable (as shown in Fig. 5) as far as the level of agreement with the PWA is concerned, except for the ${}^{1}D_{2}$ channel where the high-momentum tail is noticeably improved.

For our perturbative analysis, we have also compared the computed scattering observables at different orders in NLEFT with the continuum results of Ref. Epelbaum:2014efa . We find that our $S$-waves agree with the continuum results (within errors) up to at least $p_{\mathrm{CM}}\simeq 100$ MeV, and in some cases over the entire range of momenta considered. The $P$-waves show good agreement within errors only for some channels, and only for NLO/NNLO. As already mentioned, this is mainly due to the non-optimal description of OPE at LO. For the $D$-wave channels, only ${}^{3}D_{1}$ shows good agreement with the continuum calculations. For the ${}^{1}D_{2}$ channel, the LO and NLO results overshoot the continuum error band, while the NNLO result is in agreement due to the large uncertainty. For ${}^{3}D_{2}$, the NLO/NNLO terms do not contribute at all and hence cannot improve the result. Further, for ${}^{3}D_{3}$ the lattice calculations start to deviate from the PWA and the continuum results for $p_{\mathrm{CM}}>100$ MeV.

Finally, it is important to stress that for cms momenta below 150 MeV, the phase shifts agree within the uncertainties (with the exception of $\epsilon_{1}$, were deviations set in at about 110 MeV). This validates the statements made in Ref. Klein:2015vna about the lattice spacing independence of observables in the two-nucleon sector.

III.4 Further improvements

Next, we shall discuss two problems that require further study to resolve. First, while a clear improvement was observed in the non-perturbative case for the scattering observables (and fitted parameters) as $a$ was decreased, a similar improvement is not found in the perturbative analysis. At LO, the quality of the description improves in general with decreasing $a$, particularly for the ${}^{3}S_{1}$-${}^{3}D_{1}$ coupled channel. However, at NLO/NNLO the picture is more complicated. We note that the $P$-waves and remaining $D$-wave channels do improve, but the description of the ${}^{1}S_{0}$ and ${}^{3}S_{1}$-${}^{3}D_{1}$ channels may in fact deteriorate for smaller $a$. We attribute this effect to the increasing influence of the TPE potential. While the effect of TPE on the $S$-waves can be absorbed by smeared contact interactions as was done in the non-perturbative calculation, in the perturbative case we only have standard (without smearing) contact interactions available. This is sufficient for $a=1.97$ fm, as the TPEP contribution then closely resembles a contact interaction. A possible solution for smaller $a$ would be to include a smeared version of the NLO/NNLO contact interactions. Alternatively, one could use exact momentum operators for the NLO contact terms, which do have a higher influence at larger momenta. This was not necessary nor observable for $a=1.97$ fm, but may improve the ${}^{3}S_{1}$ channel once $a$ is decreased. Finally, we note that the choice of $c_{1}$, $c_{3}{}$ and $c_{4}$ may also have an effect, as it influences the strength of the different contribution to TPEP. However, to use the full power of chiral EFT, one should utilize the values determined from pion-nucleon scattering.

Second, we show preliminary results including a perturbative improvement of the OPE operator. In order to remedy the aforementioned discrepancies in the peripheral partial waves such that consistency with the Monte Carlo calculation is maintained, we introduce a new operator at NLO that accounts for the difference between OPEP with the momenta of Eqs. (62) and (63) and the “exact” lattice momentum $\vec{q}_{\mathrm{ex}}\equiv 2\pi\vec{k}/L$. This gives

[TABLE]

so that by adding $\Delta V_{\mathrm{OPE}}$ to $V_{\mathrm{OPE}}$, one recovers OPEP with the exact momentum. It should be noted that this differs slightly from treating OPEP at LO with the exact momentum, since $\Delta V_{\mathrm{OPE}}$ is treated as a perturbation, while $V_{\mathrm{OPE}}$ is implemented non-perturbatively. Also, $\vec{q}$ approaches $\vec{q}_{\mathrm{ex}}$ as $a\to 0$. This means that, simultaneously, $\Delta V_{\mathrm{OPE}}$ becomes less important, and $V_{\mathrm{OPE}}$ gives a better description of the $P$-waves, as we approach the continuum limit. This is consistent with Figs. 2-4 of the non-perturbative calculation, where the $P$-waves clearly improve as $a$ decreases.

Our perturbative results with $\Delta V_{\mathrm{OPE}}$ included are given in Fig. 8 and Table 10, where as expected one can observe a clear improvement in the description of the $P$-waves. The experimental results for the ${}^{1}P_{1}$, ${}^{3}P_{0}$ and ${}^{3}P_{1}$ channels are now well reproduced for the range of fitted momenta $p_{\mathrm{CM}}<150$ MeV. In general, we find that all the $P$-wave channels and the $\epsilon_{1}$ mixing angle appear much closer to the PWA at NLO with improved OPE, than without this correction. Additionally, we find that the $D$-waves (except for the ${}^{3}D_{3}$ channel) also improve significantly with respect to the LO result. In the case of ${}^{3}D_{3}$, the correction is too large and so the computed values fall below the PWA ones. Again, this improvement is mostly attributable to $\Delta V_{\mathrm{OPE}}$, although we recall that the leading (NLO) TPEP also contributes to the high-momentum tails in some of the $D$-wave channels.

III.5 Nuclear binding energies

In Monte Carlo simulations of NLEFT, the binding energies of nuclei receive perturbative energy shifts that depend on the NLO constants $C_{i}$ and their uncertainties, in addition to any inherent Monte Carlo uncertainties. For instance, in Ref. Lahde:2013uqa , only the Monte Carlo errors were taken into account, and the $C_{i}$ were assumed to be accurately known and uncorrelated. Since our analysis provides us with the complete variance-covariance matrix of the NLO parameters $C_{i}$, we are now in a position to estimate the uncertainties of the nuclear binding energies at NNLO, due to uncertainties and correlations of the $C_{i}$. From our present results, we observe larger correlations between $\Delta C$ and $C_{q}^{2}$, between $\Delta C_{I}^{2}$ and $C_{I^{2},q^{2}}$, between $\Delta C_{I}^{2}$ and $C_{S^{2},q^{2}}$, and also between $C_{(q\cdot S)^{2}}$ and $C_{I^{2},(q\cdot S)^{2}}$.

In order to obtain a first, rough estimate of the relative magnitude of Monte Carlo and fitting errors in calculations of nuclear binding energies $E_{B}$, we recall that these are calculated according to

[TABLE]

where summation over $i$ is assumed. In the Monte Carlo calculation, the LO binding energies are computed non-perturbatively, and the second term in Eq. (99) represents the perturbative shift due to the NLO constants $C_{i}$ in the 2NF, which we take from Ref. Lahde:2013uqa . We note that

[TABLE]

is a function of all the coupling constants up to NNLO, while the LO values

[TABLE]

equal the binding energies at $C_{i}=0$. In terms of the variance-covariance matrix from the perturbative analysis in Section III.3,

[TABLE]

gives us the uncertainties in the NNLO energy shifts due to the fitting errors of the $C_{i}$. The results so obtained are given in Table 11.

We note that the errors due to the uncertainties in the $C_{i}$ are of comparable magnitude to the Monte Carlo errors, even when $\mathcal{E}_{ij}$ has been evaluated without consideration of the systematical errors encoded by the Birge factor. This may suggest that the procedure of fixing the $C_{i}$ from two-nucleon data may, at present, be the main factor limiting the accuracy of NLEFT calculations beyond LO for heavier nuclei. This issue is currently under further investigation. It should also be noted that the quoted NLEFT binding energies in Table 11 are not expected to coincide with the empirical ones, as the $3N$ and higher-order contributions have been neglected (see Ref. Lahde:2013uqa for further discussion).

IV Summary

We have revisited the problem of neutron-proton scattering in NLEFT using the recently developed radial Hamiltonian method. For the first time, this has allowed us to perform a comprehensive and systematical analysis of neutron-proton phase shifts and mixing angles up to NNLO in the EFT expansion, and at several different lattice spacings in the range $1-2$ fm. We have also presented a comparison of fully non-perturbative NNLO calculations with a perturbative treatment of contributions beyond LO. Decreasing the lattice spacing to $a\sim 1$ fm necessitated the inclusion of TPEP at NLO and NNLO, and as a consequence the latter is now distinct from the NLO treatment, although no new adjustable two-nucleon parameters are introduced at NNLO.

By decreasing the lattice spacing $a$, we have found that a much improved description of neutron-proton scattering can be obtained for larger center-of-mass momenta. By considering the lattice spacings $a=1.97$ fm, $a=1.64$ fm, $a=1.32$ fm and $a=0.98$ fm, we found that $a=1.97$ fm provides a good description up to $p_{\mathrm{CM}}\simeq 100$ MeV, whereas the results for $a=0.98$ fm are reliable up to $p_{\mathrm{CM}}\simeq 200$ MeV. In general, the systematical errors are much reduced as smaller lattice spacings. Our results suggest that the range of applicability in most channels could be significantly extended by improving the lattice momenta in the OPE and NLO contact interactions. Most importantly, however, is the finding that for momenta $p_{\mathrm{CM}}\lesssim 100$ MeV, the physics of the two-nucleon system is independent of the lattice spacing $a$, when $a$ is varied in the range from 1 fm to 2 fm. Furthermore, we have also investigated the error propagation of the uncertainties of the four-nucleon LECs into the binding energies of alpha-type nuclei up to 28Si.

There are several directions in which the present work should be extended. The inclusion of N3LO contact terms and TPE contributions is underway and will appear in a separate publication, along with the inclusion of electromagnetic effects. Also, the preliminary error analysis presented here will be investigated further in a subsequent publication, especially for the propagation of the variances and covariances of the fitted NLO constants to the nuclear binding energies. Also, as pointed out in Ref. Elhatisari:2016owd , different smearing procedures and fitting also to scattering processes with more nucleons allows one to taylor interactions that might be preferable in larger nuclear systems.

Acknowledgements.

We are grateful to Serdar Elhatisari, Evgeny Epelbaum and Hermann Krebs for useful discussions. We acknowledge partial financial support from the Deutsche Forschungsgemeinschaft (Sino-German CRC 110), BMBF (Grant No. 05P12PDTEE), the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under contracts DE-FG02-03ER41260 and DE-AC05-06OR23177, the Magnus Ehrnrooth Foundation of the Finnish Society of Sciences and Letters, MINECO (Spain), and the ERDF (European Commission) grant FPA2013-40483. The work of UGM was also supported by the Chinese Academy of Sciences (CAS) President’s International Fellowship Initiative (PIFI) (Grant No. 2017VMA0025).

Appendix A Density and current operators

Here, we define the various nucleon density and current operators that we use throughout our discussion of the nucleon-nucleon interaction in the main text. Following Refs. Borasoy:2006qn ; Borasoy:2007vi , we define the local density operator

[TABLE]

the local isospin density operator

[TABLE]

the local spin density operator

[TABLE]

and the local isospin-spin density operator

[TABLE]

where $\sigma_{S}$ and $\tau_{I}$ denote the Pauli matrices for spin and isospin, respectively. Similarly, we define the current density operator

[TABLE]

the isospin-current density operator

[TABLE]

the spin-current density operator

[TABLE]

and the spin-isospin-current density operator

[TABLE]

and we recall that the current operators are used in the isospin-projected spin-orbit term (28), with $\nabla_{l}$ defined according to Eq. (29).

Appendix B Uncertainty analysis

From the definition of $\chi^{2}$ given in Eq. (88), we note that $\chi^{2}$ is a function of the LO and NLO coupling constants

[TABLE]

such that if $\chi^{2}$ is expanded around its minimum, one finds

[TABLE]

where the Hessian matrix is given by

[TABLE]

and $C_{i}^{\mathrm{min}}$ denotes the set of parameters that minimizes the $\chi^{2}$ function. Given that $\chi^{2}$ reaches its minimum value for $C_{i}=C_{i}^{\mathrm{min}}$, the terms with one derivative vanish. Keeping terms up to second order, we obtain the Hessian approximation to the error (or variance-covariance) matrix

[TABLE]

and the standard deviations

[TABLE]

of the fitted constants are obtained from the diagonal elements of the error matrix.

In the absence of systematical errors, we expect to find a normalized chi-square of $\tilde{\chi}^{2}\equiv\chi^{2}/N_{\mathrm{dof}}\approx 1$, where $N_{\mathrm{dof}}$ is the number of degrees of freedom (number of fitted data - number of free parameters) in the fit. However, in our analysis $\tilde{\chi}^{2}>1$ in most cases, particularly at LO and for larger values of the lattice spacing $a$. Such a systematical error suggests that the uncertainties computed from Eq. (115) are underestimated. Following Ref. Perez:2014yla , we therefore rescale the input errors by the Birge factor birge , according to

[TABLE]

which leads to the replacement

[TABLE]

such that $\chi^{2}/N_{\mathrm{dof}}\approx 1$ for $C_{i}=C_{i}^{\mathrm{min}}$. For a given observable $\mathcal{O}$, we assign an uncertainty according to

[TABLE]

where

[TABLE]

is the Jacobian vector of $\mathcal{O}$ with respect to the $C_{i}$.