Bayesian Analysis of $\chi$EFT at Leading Order in a Modified Weinberg Power Counting Approach

TL;DR

This paper applies Bayesian methods to analyze nucleon-nucleon interactions in chiral effective field theory at leading order, revealing insights into low-energy constants, cutoff dependence, and phase shifts, with implications for nuclear predictions.

Contribution

It introduces a Bayesian framework for analyzing $ ext{chi}$EFT interactions at leading order, addressing cutoff dependence, limit-cycle behavior, and spurious bound states.

Findings

No significant cutoff dependence in predictive distributions.

Multimodal posterior with an insignificant bound state mode.

Less constrained, more repulsive phase shifts compared to phase shift optimization.

Abstract

We present a Bayesian analysis of renormalization-group invariant nucleon-nucleon interactions at leading order in chiral effective field theory ($\chi$EFT) with momentum cutoffs in the range 400--4000 MeV. We use history matching to identify relevant regions in the parameter space of low-energy constants (LECs) and subsequently infer the posterior probability density of their values using Markov chain Monte Carlo. All posteriors are conditioned on experimental data for neutron-proton scattering observables and we estimate the $\chi$EFT truncation error in an uncorrelated limit. We do not detect any significant cutoff dependence in the posterior predictive distributions for two-nucleon observables. For all cutoff values we find a multimodal LEC posterior with an insignificant mode harboring a bound $^1{S}_0$ state. The $^3P_0$ and $^3P_2$ phase shifts emerging from the Bayesian analysis are less constrained and typically more repulsive compared to the results of a phase shift optimization. We expect that our inference will impact predictions for nuclei. This work demonstrates how to perform inference in the presence of limit-cycle-like behavior and spurious bound states, and lays the foundation for a Bayesian analysis of renormalization-group invariant $\chi$EFT interactions beyond leading order.

Figures

a](#S4.F6.sf1) we estimate that the EFT truncation error for $p\lesssim m_{\pi}$ is about $15\%$. The ppd for the spin-polarization observable $P_{b}$, which was not included in the MCMC sampling, reproduces experimental data rather poorly. This result is particularly striking in the ppd at $T_{\mathrm{lab}}=25.0$ MeV. The situation improves at higher energies, which is somewhat surprising for a low-energy EFT. As such, conditioning on low-energy $P_{b}$ data could improve the model. However, expanding the dataset in [Table 1](#S3.T1) requires an improved model for the EFT truncation error [[2

Tables

Table 1. Table 1: Total and differential n ​ p 𝑛 𝑝 np scattering observables used to condition the Bayesian inference in this work. The total number of data points is 1043. This is a subset of the data in the Granada database [ 65 , 66 ] , where the normalizations determined from the Granada analysis are included. The uncertainties in the normalizations are not taken into account. The definitions of the various scattering observables can be found in Ref. [ 56 ] .

Observable	$T_{\mathrm{lab}}$ [MeV]	# of data points
$\sigma$	${10}^{-6}$ – 99.0	324
$d\sigma /d\Omega$	2.72 – 99.0	698
${\sigma }_T$	3.65 – 17.1	13
${\sigma }_L$	4.98 – 66.0	8

Table 2. Table 2: n ​ p 𝑛 𝑝 np scattering phase shifts in degrees, used as observational data within the first wave of history matching. Values for T lab subscript 𝑇 lab T_{\mathrm{lab}} are given in MeV. The phase shift at T lab = 40 subscript 𝑇 lab 40 T_{\mathrm{lab}}=40  MeV was obtained through spline interpolation. All other phase shifts were taken from Refs.  [ 65 , 66 ] . The standard deviation for the theoretical error in Eq. ( 21 ) is presented in the last column as a relative error.

$T_{\mathrm{lab}}$	${\delta }_{{}_{}{}^{1}S_0^{}}$	${\delta }_{{}_{}{}^{3}S_1^{}}$	${\delta }_{{}_{}{}^{3}P_0^{}}$	${\delta }_{{}_{}{}^{3}P_2^{}}$	${\sigma }_{\mathrm{th}}/\delta$
1	62.1	147.6	0.181	0.022	0.3
5	63.7	118.0	1.66	0.258	0.3
10	60.0	102.3	3.75	0.727	0.3
25	51.1	80.2	8.48	2.63	0.3
40	44.8	67.6	10.7	4.70	0.3
100	27.2	42.7	9.14	11.06	0.8

Table 3. Table 3: n ​ p 𝑛 𝑝 np scattering observables used within the second wave of history matching. Note that σ 𝜎 \sigma , without any subscripts, is denoting the total n ​ p 𝑛 𝑝 np cross section, while σ th / y exp subscript 𝜎 th subscript 𝑦 exp \sigma_{\mathrm{th}}/y_{\mathrm{exp}} is the (relative) model error for each observable.

Observable	$T_{\mathrm{lab}}$ [MeV]	${\theta }_{c.m.}$ [deg]	$y_{\exp }$	${\sigma }_{\mathrm{th}}/y_{\exp }$
$\sigma$	0.12	-	12050.0 mb	0.3
$\sigma$	3.186	-	2206.0 mb	0.3
$\sigma$	12.995	-	749.0 mb	0.3
$\sigma$	28.0	-	321.5 mb	0.3
$\sigma$	40.0	-	217.8 mb	0.3
$\sigma$	97.2	-	76.0 mb	0.8
$d\sigma /d\Omega$	99.0	21.0	8.01 mb/sr	0.8
$d\sigma /d\Omega$	99.0	78.0	2.14 mb/sr	0.8
$d\sigma /d\Omega$	99.0	149.0	9.50 mb/sr	0.8
$P_b$	25.0	33.1	0.047	0.3
$P_b$	25.0	90.3	0.053	0.3
$P_b$	95.0	29.8	0.170	0.8
$P_b$	95.0	88.5	0.291	0.8

Table 4. Table 4: The model evidences and maximum a posteriori (MAP) estimates ( 𝜽 ∗ superscript 𝜽 \bm{\theta}^{*} ) of the LECs for all cutoffs Λ Λ \Lambda (in MeV). A posterior bimodality with comparable probability mass per mode exists for both Λ = 700 Λ 700 \Lambda=700 and Λ = 4000 Λ 4000 \Lambda=4000 MeV.

$\Lambda$	$\ln \left({𝒵}_K\right)$	MAP (${𝜽}^{\ast })$
400	-3280	(-0.1186, -0.1167, 5.386, 0.5962, 2.716)
450	-3312	(-0.1128, -0.09640, 5.115, 0.5118, 2.813)
500	-3325	(-0.1084, -0.07815, 4.713, 0.4560, 2.859)
550	-3305	(-0.1048, -0.06069, -1.534, 0.7440, 2.808)
600	-3296	(-0.1018, -0.04352, -0.8692, 0.9547, 2.790)
700	-3282	(-0.09706, -0.006363, -0.3682, -2.848, 2.776)
700	-3280	(-0.09707, -0.006305, -0.3803, 1.910, 2.761)
2000	-3278	(-0.07831, -0.1267, 0.01847, -0.02938, 2.745)
3000	-3286	(-0.07482, -0.03519, 0.02497, -0.01020, 2.763)
4000	-3283	(-0.07210, 0.09648, -1.140, -0.005023, 2.773)
4000	-3282	(-0.07210, 0.09641, 0.2060, -0.005050, 2.773)

Table 5. Table 5: Threshold LEC values for which a bound state appears in the partial waves relevant for the contact potential at LO. For the S 1 3 − 3 D 1 superscript 3 superscript subscript 𝑆 1 3 subscript 𝐷 1 {}^{3}S_{1}\mathrm{-}^{3}D_{1} channel and Λ ≥ 2000 Λ 2000 \Lambda\geq 2000 MeV it is the second bound state that appears at the tabulated value, and the first appears for C ~ S 1 3 > 0.4 subscript ~ 𝐶 superscript subscript 𝑆 1 3 0.4 \tilde{C}_{{}^{3}S_{1}}>0.4 . For the P 0 3 superscript subscript 𝑃 0 3 {}^{3}P_{0} channel and Λ = 4000 Λ 4000 \Lambda=4000 MeV it is the second bound state that appears at the tabulated value, and the first appears for C P 0 3 > 1 subscript 𝐶 superscript subscript 𝑃 0 3 1 C_{{}^{3}P_{0}}>1 .

$\Lambda$ [MeV]	${\tilde{C}}_{{}_{}{}^{1}S_0^{}}$ [${10}^4$ GeV-2]	${\tilde{C}}_{{}_{}{}^{3}S_1^{}}$ [${10}^4$ GeV-2]	$C_{{}_{}{}^{3}P_0^{}}$ [${10}^4$ GeV-4]	$C_{{}_{}{}^{3}P_2^{}}$ [${10}^4$ GeV-4]
400	$-0.1279$	$-0.0825$	$-1.6735$	$-2.3178$
450	$-0.1207$	$-0.0636$	$-1.1076$	$-1.6375$
500	$-0.1151$	$-0.0449$	$-0.7578$	$-1.2016$
550	$-0.1106$	$-0.0254$	$-0.5316$	$-0.9095$
600	$-0.1067$	$-0.0037$	$-0.3802$	$-0.7063$
700	$-0.1010$	$0.0515$	$-0.2015$	$-0.4528$
2000	$-0.0794$	$-0.0983$	$0.0268$	$-0.0263$
3000	$-0.0755$	$0.0010$	$0.0449$	$-0.0096$
4000	$-0.0735$	$0.3229$	$-0.0310$	$-0.0048$

Table 6. Table 6: Modes found in the initial MCMC sampling of the posterior pdfs. For each cutoff there are two starting regions ( s 𝑠 s ), indicated in the second column, one where S 0 1 superscript subscript 𝑆 0 1 {}^{1}S_{0} contains a bound state (b) and one where it does not (u).

$\Lambda$ [MeV]	$s$	$\ln \left({𝒵}_K\right)$	MAP (${𝜽}^{\ast })$
400	u	-3280	(-0.1186,-0.1168, 5.3536, 0.5890, 2.7193)
400	b	-3887	(-0.1594,-0.1214, 3.1047, 0.4464, 4.9501)
450	u	-3312	(-0.1129,-0.0964, 5.0422, 0.5193, 2.8031)
450	b	-3879	(-0.1448,-0.1000, 3.6050, 0.4650, 4.9108)
450	b	-3893	(-0.1445,-0.0995,-6.2068, 0.5853, 4.9336)
500	u	-3325	(-0.1084,-0.0781, 4.9732, 0.4489, 2.8520)
500	b	-3861	(-0.1340,-0.0807,-3.9206, 0.6015, 4.8125)
500	b	-3862	(-0.1341,-0.0809, 4.0452, 0.5022, 4.8321)
550	u	-3305	(-0.1048,-0.0607,-1.5204, 0.7656, 2.8144)
550	u	-3328	(-0.1048,-0.0607, 4.4330, 0.3724, 2.8492)
550	b	-3832	(-0.1259,-0.0629,-2.3130, 0.7294, 4.6720)
550	b	-3843	(-0.1260,-0.0633, 4.1715, 0.5688, 4.7062)
600	u	-3296	(-0.1018,-0.0436,-0.8677, 1.0069, 2.7920)
600	u	-3326	(-0.1017,-0.0438, 4.2580, 0.3304, 2.8839)
600	b	-3822	(-0.1196,-0.0453,-1.2360,-4.8906, 4.5958)
600	b	-3841	(-0.1197,-0.0455, 4.2284,-4.7557, 4.6951)
600	b	-3828	(-0.1196,-0.0455,-1.2646, 0.8954, 4.5454)
600	b	-3824	(-0.1197,-0.0456, 4.2902, 0.6672, 4.6155)
700	u	-3279	(-0.0971,-0.0063,-0.3798, 1.9723, 2.7659)
700	u	-3282	(-0.0971,-0.0064,-0.3771,-2.9677, 2.7711)
700	u	-3323	(-0.0970,-0.0066, 3.7883, 0.2950, 2.8984)
700	b	-3766	(-0.1105,-0.0078,-0.4676,-2.5175, 4.3600)
700	b	-3795	(-0.1105,-0.0086, 3.8219, 1.1784, 4.5288)
700	b	-3764	(-0.1105,-0.0080,-0.4887, 1.9523, 4.3593)
2000	u	-3278	(-0.0783,-0.1267, 0.0184,-0.0293, 2.7490)
2000	b	-3688	(-0.0812,-0.1266, 0.0165,-0.0285, 3.9715)
3000	u	-3286	(-0.0748,-0.0352, 0.0248,-0.0102, 2.7581)
3000	b	-3688	(-0.0767,-0.0350, 0.0222,-0.0100, 3.9804)
4000	u	-3279	(-0.0730, 0.0967,-1.0424,-0.0050, 2.7585)
4000	u	-3282	(-0.0730, 0.0964, 0.2322,-0.0050, 2.7627)
4000	b	-3683	(-0.0744, 0.0966, 0.0819,-0.0050, 3.9644)

Table 7. Table 7: Phase shifts and mixing angles in degrees for relevant partial waves.

$T_{\mathrm{lab}}$ [MeV]	${\delta }_{{}_{}{}^{1}S_0^{}}$	${\delta }_{{}_{}{}^{3}S_1^{}}$	${\delta }_{{}_{}{}^{3}D_1^{}}$	${\delta }_{{ϵ}_1}$	${\delta }_{{}_{}{}^{3}P_0^{}}$	${\delta }_{{}_{}{}^{3}P_2^{}}$	${\delta }_{{}_{}{}^{3}F_2^{}}$	${\delta }_{{ϵ}_2}$
1	$48.11269$	$144.56170$	$-0.00641$	$0.17917$	$0.13249$	$-0.01180$	$0.00001$	$-0.00154$
50	$55.43553$	$67.70410$	$-7.98257$	$4.45167$	$-2.94614$	$-3.27873$	$0.31327$	$-1.63134$
100	$47.35263$	$51.13999$	$-16.12281$	$6.19131$	$-18.21728$	$-9.11571$	$0.70523$	$-2.73674$
200	$32.27681$	$31.20600$	$-26.46882$	$8.71350$	$-51.24600$	$-21.64237$	$1.00719$	$-3.46681$
300	$19.60329$	$17.43370$	$-26.71548$	$9.20362$	$-82.25737$	$-29.45361$	$0.80583$	$-3.12625$

Equations

$$y_{\mathrm{th}}^{(n)}$$

$$V({\mathbf{p}}^{\prime },\mathbf{p})$$

$$\displaystyle+\tilde{C}_{{}^{1}S_{0}}+\tilde{C}_{{}^{3}S_{1}}+\left(C_{{}^{3}P_{0}}+C_{{}^{3}P_{2}}\right)p^{\prime}p\Bigg{]}.$$

$$f_{\Lambda }(p)=\exp \left[-\frac{p^4}{{\Lambda }^4}\right]\times \theta (\Lambda +(300\mathrm{MeV})-p),$$

$$\theta (x)=\{\begin{array}{l}1,x>0\\ 0,x\le 0\end{array}.$$

$$V({\mathbf{p}}^{\prime },\mathbf{p})\to f_{\Lambda }(p^{\prime })\sqrt{\frac{m_N}{E(p^{\prime })}}V({\mathbf{p}}^{\prime },\mathbf{p})\sqrt{\frac{m_N}{E(p)}}{f}_{\Lambda }(p),$$

$$T_{l^{\prime }l}^{sj}(p^{\prime },p)=V_{l^{\prime }l}^{sj}(p^{\prime },p)+$$

$$+\sum _{l^{\prime \prime }}{\int }_0^{\infty }dk{k}^2{V}_{l^{\prime }{l}^{\prime \prime }}^{sj}(p^{\prime },k)\frac{m_N}{p^2-k^2+iϵ}{T}_{l^{\prime \prime }l}^{sj}(k,p).$$

$$\mathrm{pr}\left(\mathbf{\theta }|D,I\right)=\frac{\mathrm{pr}\left(D|\mathbf{\theta },I\right)\times \mathrm{pr}\left(\mathbf{\theta }|I\right)}{\mathrm{pr}\left(D|I\right)}.$$

$$\mathrm{pr}\left(\mathbf{\theta }|D,I\right)\propto \mathrm{pr}\left(D|\mathbf{\theta },I\right)\times \mathrm{pr}\left(\mathbf{\theta }|I\right).$$

$$\mathbf{\theta }=\left(\mathbf{\alpha },\bar{c}\right),$$

$$\mathbf{\alpha }=\left({\tilde{C}}_{{}^1{S}_0},{\tilde{C}}_{{}^3{S}_1},C_{{}^3{P}_0},C_{{}^3{P}_2}\right)$$

$$y_{\text{exp}}=y_{\mathrm{th}}^{(n)}+\delta y_{\mathrm{th}}+\delta y_{\text{exp}},$$

$$y=y_{\mathrm{ref}}\sum _{k=0}^{\infty }{b}_k{\left(\frac{Q}{{\Lambda }_{\chi }}\right)}^k,$$

$$y_{\mathrm{th}}^{(0)}=y_{\mathrm{ref}}\Bigg{[}b_{0}+\left(\frac{Q}{\Lambda_{\chi}}\right)\mathcal{D}_{1}\Bigg{]},$$

$$y=y_{\text{th}}^{(0)}+y_{\text{ref}}\sum _{k=2}^{\infty }{c}_k{\left(\frac{Q}{{\Lambda }_{\chi }}\right)}^k.$$

$$\delta y_{\mathrm{th}}\equiv y_{\text{ref}}\sum _{k=2}^{\infty }{c}_k{\left(\frac{Q}{{\Lambda }_{\chi }}\right)}^k.$$

$$\mathrm{pr}\left(c_k|I\right)=\mathcal{N}\left(0,{\bar{c}}^2\right),$$

$$\mathrm{pr}\left(\delta y_{\mathrm{th}}|I\right)$$

$$\text{with }{\sigma }_{\mathrm{th}}^2$$

$$\mathrm{pr}\left(\delta y_{\exp }|I\right)=\mathcal{N}\left(0,{\sigma }_{\exp }^2\right).$$

$$\mathrm{pr}\left(y_{\exp }|\mathbf{\theta },{\sigma }_{\exp },I\right)=\mathcal{N}\left(y_{\mathrm{th}}^{(0)}(\mathbf{\alpha }),{\sigma }_{\mathrm{th}}^2(\bar{c})+{\sigma }_{\exp }^2\right),$$

$$\mathrm{pr}\left(\bar{c}|I\right)=\frac{{\beta }^{\alpha }}{\Gamma (\alpha )}(1/\bar{c}{)}^{\alpha +1}\exp \left(-\beta /\bar{c}\right),$$

$$I_M(\mathbf{\alpha })\equiv \underset{i}{\max }\sqrt{\frac{\text{absolutevalue}{y_{\mathrm{th},i}^{(0)}\left(\mathbf{\alpha }\right)-y_{\exp ,i}}^2}{{\sigma }_{\exp ,i}^2+{\sigma }_{\mathrm{th},i}^2}}.$$

$${\mathcal{Z}}_K=\mathrm{pr}{\left(D|I\right)}_K$$

$$\approx \mathrm{pr}\left({\mathbf{\theta }}^{\ast }|D,I\right)\frac{{\left(2\pi \right)}^{N/2}}{\sqrt{\det {\Sigma }^{-1}}},$$

$$\mathrm{pr}\left(y|D,I\right)=\int d\mathbf{\theta }\mathrm{pr}\left(y|\mathbf{\theta },D,I\right)\mathrm{pr}\left(\mathbf{\theta }|D,I\right).$$

$$\mathrm{pr}\left(y|\mathbf{\theta },D,I\right)=\mathcal{N}\left(y_{\mathrm{th}}^{(0)}(\mathbf{\alpha }),{\sigma }_{\mathrm{th}}^2(\bar{c})\right),$$

$$\mathrm{pr}\left(y|\mathbf{\theta },D,I\right)=\delta \left(y-y_{\mathrm{th}}^{(0)}\left(\mathbf{\alpha }\right)\right).$$

$$p^2=\frac{m_p^2{T}_{\mathrm{lab}}(2m_n+T_{\mathrm{lab}})}{(m_n+m_p{)}^2+2m_p{T}_{\mathrm{lab}}}.$$

$$⟨{\mathbf{p}}^{\prime }||l{m}_l,p⟩=\frac{\delta (p^{\prime }-p)}{p^2}{Y}_{m_l}^l({\hat{p}}^{\prime }),$$

$$\int d\hat{p}{Y}_{m_l^{\prime }}^{\ast l^{\prime }}(\hat{p})Y_{m_l}^l(\hat{p})={\delta }_{l^{\prime }l}{\delta }_{m_l^{\prime }{m}_l}.$$

Full text

I Introduction

The foundational principles of chiral effective field theory ($\chi$EFT) [1, 2, 3] promise a systematically improvable description of the strong force between nucleons that is consistent with quantum chromodynamics (QCD). However, establishing a power counting of chiral nuclear interactions that fulfills renormalization requirements presents several challenges—see, e.g., van Kolck [4] for an overview. Recently, Yang et al. [5] analyzed nuclear ground-state energies using renormalization-group (RG) invariant formulations of $\chi$EFT up to (perturbative) next-to-leading order (NLO) corrections. Apparently, the essential mechanism for nuclear-binding tends to fail already at leading order (LO) for selected $A>4$ nuclei when using available RG-invariant power counting schemes. In fact, similar results had already been found using $\chi$EFT [6] based on the canonical Weinberg power counting (WPC) [7, 8, 9, 10], and pionless EFT [11, 12, 13]. Yang et al. [5] put forward three possible reasons for these shortcomings at LO in $\chi$EFT: i) One (or more) scales critical to the physical description of finite nuclei might not be correctly captured by the contact terms at LO [14]; ii) The LO nucleon-nucleon ($NN$) interaction should possibly be complemented with other interaction terms such as sub-leading pion exchange [15] and many-nucleon interactions [16, 17, 18, 19]; iii) The description of the nuclear interaction might be finely tuned and therefore require careful calibration of the low energy constants (LECs). Indeed, Yang et al. [5] renormalized the relevant LECs by demanding exact reproduction of selected phase shifts at a single scattering energy. This procedure resulted in point estimates of the LECs without the possibility to analyze possible fine-tuning effects.

In this work we tackle overfitting and expose possible fine tuning using Bayesian methods. Specifically, we infer a posterior probability density function (pdf) for the values of the LECs conditioned on neutron-proton ($np$) scattering data. The advantages of a Bayesian approach are several. First, we obtain a probabilistic measure of our uncertainty about the values of the LECs and subsequent predictions, something that is not obtained when doing maximum likelihood estimation [6]. Second, with the Bayesian approach we can utilize the expected systematicity of $\chi$EFT as prior knowledge about the truncation error [20]. There exist several Bayesian studies of $\chi$EFT interactions [21, 22, 23, 24, 25] and quantified posterior predictive distributions (ppds) of nuclear properties [26, 27, 28]. However, so far all such studies are grounded in $\chi$EFT based on WPC.

In EFTs that are perturbative at all orders, power counting usually follows the momentum scaling of the different Feynman diagrams. This is known as naive dimensional analysis (NDA) [29]. However, in $\chi$EFT we must account for bound multi-nucleon states and therefore face non-perturbative physics with the consequence of infrared enhancement of purely nucleonic intermediate states. To deal with this, Weinberg [2, 3] suggested to apply NDA to the potential which is then iterated to all orders by solving, e.g., the Lippmann-Schwinger equation. This prescription assumes that the iteration to infinite order does not introduce additional divergences with the need for higher-order counterterms—an assumption that did not hold upon closer inspection. Indeed, taking the momentum cutoff $\Lambda$ of the regulator very large, i.e., far beyond the anticipated breakdown scale of $\chi$EFT, Nogga et al. [30] found that WPC does not yield RG-invariant amplitudes in attractive spin-triplet partial-waves. By now there exist several proposals on how to modify WPC, see, e.g., Refs. [31, 32, 33, 5, 18, 34, 35, 36, 37]. We refer to all such proposals as modified Weinberg power counting (MWPC). One can argue that WPC provides a consistent EFT framework as long as $\Lambda$ is kept in the vicinity of the breakdown scale $\Lambda_{\chi}$ [38, 39, 40, 41]. In that scheme, all orders are typically summed up non-perturbatively and, starting at third order, it can provide realistic predictions for selected nuclei [42] and nuclear matter [43]. However, this achievement does not imply that WPC provides the foundation for an EFT of QCD.

This paper is organized as follows: In Sec. II the LO potential and non-relativistic $np$ scattering are addressed along with relevant background regarding power counting and the renormalization of singular potentials. In Sec. III we discuss our Bayesian approach and in Sec. IV we review the numerical sampling of the posterior pdf for the LECs. Finally, in Sec. V ppds are presented and analyzed for some $np$-scattering observables and the deuteron ground-state energy. A conclusion and outlook is presented in Sec. VI.

II Theory

To keep this work self-contained we first discuss how counterterms (and associated LECs) are introduced to renormalize the singular nature of the one-pion exchange (OPE) potential at LO in $\chi$EFT. We also study limit-cycle-like behavior in some detail since they play an important role during the Bayesian inference.

II.1 Effective field theory expansion

The EFT approach promises an order-by-order improvable description of a nuclear observable $y$ residing in the low-energy domain below the relevant breakdown scale. In general, up to some finite order $n$ we can expect an expansion of the form [10, 44]

[TABLE]

and we specialize to $\chi$EFT assuming a breakdown scale $\Lambda_{\chi}=600$ MeV as per previous analyses [45, 21] and denote the momentum cutoff by $\Lambda$. We also assume a soft scale $Q=\max(p,m_{\pi})$ with $p$ denoting the external momentum and $m_{\pi}$ denoting the pion mass. The unspecified function $\mathcal{D}_{n+1}$ depends on ratios of the relevant low energy scales and, importantly, absorbs the residual cutoff dependence. In this work we set the reference scale $y_{\text{ref}}$ to the corresponding experimental value $y_{\text{exp}}$ of $y$. Pulling out $y_{\text{ref}}$ leads to dimensionless expansion coefficients $b_{k}$, and if the theory is renormalized order-by-order, these should not exhibit any cutoff dependence although they do depend on ratios of relevant scales. Also, along the lines of naturalness [46], we expect the $b_{k}$-values to be of order unity.

For a perturbative EFT—where NDA applied to the Feynman diagrams carries over to the amplitudes and hence to observables—the power counting in Eq. 1 follows straightforwardly from the Lagrangian. The relation between the Lagrangian and observables in a non-perturbative theory like $\chi$EFT is less direct. The non-perturbative effects in combination with the need to treat the problem numerically pose challenges to finding a consistent power counting. When the amplitudes cannot be obtained perturbatively, the LO $T$-matrix must scale as $T\sim Q^{-1}$ [47, 48, 44]. The summation in Eq. 1 can still start at $k=0$, since the $Q^{-1}$ dependence can be absorbed into $y_{\mathrm{ref}}$. Terms with $k=1,2,\ldots$ then correspond to higher-order corrections with respect to LO. We employ MWPC where amplitudes of the LO potential are computed non-perturbatively and sub-leading corrections should be accounted for using perturbation theory [4, 49, 50, 33, 51].

II.2 Two-nucleon potential and scattering amplitudes

The momentum-space and isospin-symmetric LO potential considered in this work is adopted from Ref. [5], and using our conventions it reads

[TABLE]

The first term is the OPE potential where $\bm{\tau}_{i}$ and $\bm{\sigma}_{i}$, for $i=1,2$, are the isospin and spin operators for the respective nucleon, $\bm{p}$ ($\bm{p}^{\prime}$) are the ingoing (outgoing) relative momenta with normalization $\innerproduct{\bm{p}^{\prime}}{\bm{p}}=\delta^{3}\left(\bm{p}^{\prime}-\bm{p}\right)$, $\bm{q}=\bm{p}^{\prime}-\bm{p}$ is the momentum transfer, and we have $p\equiv|\bm{p}|$ and $p^{\prime}\equiv|\bm{p}^{\prime}|$. The contact LECs: $\tilde{C}_{{}^{1}S_{0}},\tilde{C}_{{}^{3}S_{1}},C_{{}^{3}P_{0}},C_{{}^{3}P_{2}}$, carry implicit projection operators such that they act in the indicated partial-waves, ${}^{(2s+1)}l_{j}$, where $s,l,j$ denote the quantum numbers of the $NN$ spin, orbital angular momentum, and total angular momentum, respectively. We use the PDG [52] values for the axial coupling constant $g_{A}=1.275$, the pion decay constant $f_{\pi}=92.1$ MeV, and averaged pion mass $m_{\pi}=138.039$ MeV. Finally, for the partial-wave decomposed [53] potential, and operators, we use the notation $V_{l^{\prime}l}^{sj}(p^{\prime},p)\equiv\bra{l^{\prime}sj;p^{\prime}}V\ket{lsj;p}$. As discussed in Sec. II.3 and Ref. [5], this LO potential is restricted to low partial waves: ${}^{1}S_{0}$, ${}^{3}S_{1}\mathrm{-}^{3}D_{1}$, ${}^{1}P_{1}$, ${}^{3}P_{0}$, ${}^{3}P_{1}$, ${}^{3}P_{2}\mathrm{-}^{3}F_{2}$. Consequently, it is set to zero for all partial waves with $l>1$ that have no coupling to $l\leq 1$.

To render the integrals of the Lippmann-Schwinger equation finite, the relative momenta of in- and out-going nucleons are regulated using the function,

[TABLE]

which was also used by Yang et al. [5]111The exact form of the regulator was not given in Ref [5]. Thus, we define it here., where

[TABLE]

We straightforwardly account for the small effects of relativistic kinematics using the minimal relativity prescription[54, 9]. Combined with the momentum regulators $f_{\Lambda}$, the LO potential in Eq. 2 is thus modified according to

[TABLE]

where $E(p)=\sqrt{m_{N}^{2}+p^{2}}$ and $m_{N}=\frac{2m_{p}m_{n}}{m_{p}+m_{n}}$ is the nucleon mass with proton and neutron masses $m_{p}=938.272$ MeV and $m_{n}=939.918$ MeV, respectively [52].

We condition all inferences on $np$ scattering data and must therefore solve the corresponding Lippmann-Schwinger equation for the $T$-matrix

[TABLE]

We solve this numerically in momentum space using a standard matrix-inversion method by first converting to a real equation for the reaction matrix [55]. Furthermore, we sum the partial-wave amplitudes to construct the spin scattering matrix, $M^{s}_{m^{\prime}_{s}m_{s}}(p,\theta_{\mathrm{c.m.}})$, see Appendix A. Here, $m_{s}$, ($m^{\prime}_{s}$) is the in-(out-)going total spin projection and $\theta_{\mathrm{c.m.}}$ is the center-of-mass scattering angle. All scattering observables can be computed from this spin scattering matrix [56]. The types of observables that we condition our inferences on are listed in Table 1 and discussed in Sec. III.1.

II.3 Singular potentials and limit-cycle-like behavior

An attractive potential is singular near the origin if it behaves as $-\lambda/r^{n}$ with $n>2$ (or $n=2$ with sufficiently large $\lambda>0$) [57, 58]. Two particles interacting only via a singular potential will collapse towards the origin with increasing velocity. Akin to the infinities in quantum field theory, the singularity at the origin is cured via renormalization. The OPE potential is singular in attractive spin-triplet partial waves, e.g., ${}^{3}P_{0}$, ${}^{3}S_{1}\mathrm{-}^{3}D_{1}$, ${}^{3}D_{2}$ and ${}^{3}P_{2}\mathrm{-}^{3}F_{2}$ [30]. From an EFT perspective, the singular character of the OPE becomes physically meaningful only with the addition of counterterms that parameterize, and cure, these short-range pathologies [59, 60].

To avoid the introduction of infinitely many counterterms, the OPE potential in Eq. (2) is truncated to act only in partial waves with orbital angular momentum $l<l_{c}$. It is not obvious, however, where the limit $l_{c}$ should be set. Higher-order terms can be included using the distorted-wave Born approximation and there is evidence that this does not necessitate the introduction of additional counterterms, see Refs. [37, 41] for contrasting views. According to previous studies, OPE contributions to the $NN$ scattering amplitude in partial waves with orbital angular momentum $l>1$ can be treated perturbatively up to at least $p\sim 300$ MeV [35, 61, 4]. We therefore truncate the OPE potential such that it is nonzero only for channels that contain a partial wave component with $l\leq l_{c}=1$.

LECs associated with counterterms that renormalize a singular potential with $n>2$ usually exhibit a limit-cycle-like behavior [59, 60, 30], i.e., provided that there is only one counterterm per partial wave, the corresponding LEC will exhibit periodic discontinuities as a function of the regulator value. This was extensively investigated for one-pion exchange in Ref. [30]. We can reproduce those results exactly, but we also observe a slight shift in the location of the limit-cycles-like behavior when including the minimal relativity correction in Eq. (5). In Fig. 1 we show how the LECs in the ${}^{1}S_{0}$, ${}^{3}S_{1}$, ${}^{3}P_{0}$ and ${}^{3}P_{2}$ partial waves run with $\Lambda$ when renormalized to reproduce the Nijmegen partial-wave phase shifts [62] for the laboratory kinetic energy $T_{\mathrm{lab}}=50$ MeV of the projectile nucleon. There is no limit-cycle-like behavior observed in the ${}^{3}P_{2}$ wave for the cutoff region studied here.

When this renormalization procedure is used, spurious and deeply bound states appear in the singular and attractive partial-waves as the cutoff $\Lambda$ is increased. This is not a problem if the states are deeply bound and thus outside the applicable domain of the EFT. In practice, one can project these states out of the spectrum in a phase-equivalent fashion [30]. In Fig. 2 we show how the spurious states in ${}^{3}P_{0}$ and ${}^{3}S_{1}\mathrm{-}^{3}D_{1}$ appear at the threshold value of $\Lambda$ corresponding to where the limit-cycle-like behavior is observed in Fig. 1.

In this work we infer the LEC posteriors at fixed values of the cutoff $\Lambda$. There exist threshold LEC values for which the total potential becomes sufficiently attractive such that a spurious bound state appears. The ($\Lambda$-dependent) LEC values for which this happens are tabulated in Appendix B for the respective channels. Note that a spurious state is not deeply bound in the immediate vicinity of these threshold values and we stress that the behavior depicted in Fig. 2 is obtained when requiring the exact reproduction of a specific phase shift.

As the bound state pole moves near the threshold when the corresponding LEC is varied, the phase shift will change dramatically by Levinson’s theorem [63]. Rapidly varying phase shifts lead to particular challenges in the inference process since scattering observables that are included in the likelihood, and that depend sensitively on partial-wave amplitudes, might then constrain the LEC values to very narrow regions of parameter space. In Fig. 3 we exemplify this by showing phase shifts at $T_{\mathrm{lab}}=50$ MeV as a function of the LEC in the partial waves ${}^{3}P_{0}$ and ${}^{3}P_{2}$ for different cutoff values. The red, horizontal line indicates the empirical value of the phase shift according to the Nijmegen analysis [62]. Fig. 3 also shows that for certain cutoffs we can obtain a wide range of LEC values that reproduce the empirical phase shift with reasonable precision.

Figures 1 and 3 provide complementary views on limit-cycles-like behavior. For example, in Fig. 1 we see that $C_{{}^{3}P_{0}}$ changes from a large positive value to a large negative one as we approach $\Lambda\approx 780$ MeV from below. The same information is contained in Fig. 3 where the LO ${}^{3}P_{0}$ phase shift intersects with the Nijmegen result for a positive (negative) value of $C_{{}^{3}P_{0}}$ for $\Lambda=700$ $(1000)$ MeV. Enforcing the exact reproduction of a phase shift therefore implies that the spurious bound state will be deeply bound exactly when the limit-cycle-like behavior appears as shown in Fig. 1. For the ${}^{3}P_{2}$ wave, the phase shift only intersects with the Nijmegen result on one side of the discontinuity, at least for the low values of $\Lambda$ shown here. Consequently, there is no spurious bound state, or limit-cycle-like behavior appearing.

III Bayesian parameter estimation

Bayes’ theorem expresses the posterior pdf, $\mathrm{pr}\left(\bm{\theta}|D,I\right)$, for the relevant model parameters, $\bm{\theta}$, conditioned on experimental data, $D$, and other pertinent information, $I$, in terms of quantities that we can evaluate

[TABLE]

Here, $\mathrm{pr}\left(\bm{\theta}|I\right)$ is the prior, $\mathrm{pr}\left(D|\bm{\theta},I\right)$ is the likelihood, and $\mathrm{pr}\left(D|I\right)$ is the model evidence. The latter is independent of $\bm{\theta}$ and can be omitted for parameter estimation purposes leaving us with the simpler expression

[TABLE]

The set of model parameters in this study is collected in a vector

[TABLE]

where $\bm{\alpha}$ denotes the LECs of $\chi$EFT at LO and $\bar{c}$ governs the magnitude of the EFT truncation error, as described below in Eq. 17. The elements of the vector $\bm{\alpha}$ at LO are

[TABLE]

in units of $10^{4}$ GeV*-2* and $10^{4}$ GeV*-4* for the $S$- and $P$-wave LECs, respectively.

III.1 Likelihood

When relating theory $y^{(n)}_{\text{th}}$ and experiment $y_{\text{exp}}$ for some observable $y$ we must account for both the theoretical and experimental uncertainties. This can be achieved with a statistical model

[TABLE]

in which the uncertainties are modeled by random variables $\delta y_{\text{th}}$ and $\delta y_{\text{exp}}$, respectively, and where we suppressed all parameter dependencies for simplicity. One of the great benefits of working with EFTs is the expected systematicity of the truncation error. Assuming that QCD accounts for all relevant physics, and that our EFT is sound, i.e., it systematically approximates QCD in the low-energy regime, then we can express the value of some nuclear observable $y$ as

[TABLE]

in the notation of Eq. 1. Thus, the LO theoretical prediction can be expressed as

[TABLE]

where the second term contains the residual cutoff dependence through $\mathcal{D}_{1}$. For the present case we assume that $b_{1}=0$ by time-reversal and parity invariance [8, 9, 10]. With this, Eq. 12 can be written as

[TABLE]

The coefficients $c_{k}$ are related to $b_{k}$ and the power series expansion of $\mathcal{D}_{1}$. Since $b_{1}=0$, it implies that $c_{1}=0$ if the first term in the power-series expansion of $\mathcal{D}_{1}$ vanishes, which we assume it does by the same symmetry argument as for $b_{1}$. The coefficients $c_{k}$, like $b_{k}$, are also assumed to be of natural size, i.e., of order one. Finally, we note that $c_{k}$ might have a residual cutoff dependence inherited from $\mathcal{D}_{1}$.

We can now express the LO truncation error as

[TABLE]

Quantitatively, this is unknown to us since we have not computed anything beyond LO. However, we can evaluate the probability distribution of the random variable $\delta y_{\mathrm{th}}$ if we assume a distribution for the expansion coefficients $c_{k}$ guided by domain knowledge included in $I$. In this work we follow Ref. [20] and assume that

[TABLE]

where $\mathcal{N}(\mu,\sigma^{2})$ denotes a normal distribution with mean $\mu$ and variance $\sigma^{2}$. Assuming independent and identically distributed $c_{k}$ coefficients the distribution for $\delta y_{\mathrm{th}}$ at each kinematical point is given by [64]

[TABLE]

Hence, the variance of the truncation error is governed by $\bar{c}^{2}$, which is unknown a priori. However, operating with an EFT we expect $\bar{c}$ to be of order one. We will quantify this in Sec. III.2.

The experimental error for each datum is assumed to follow a normal distribution with variance $\sigma^{2}_{\mathrm{exp}}$, i.e.,

[TABLE]

Assuming that the EFT truncation error and the experimental errors are independent the likelihood for a single observation reads

[TABLE]

by Eq. 11, where we now explicitly indicate the relevant parameter dependencies. Following Refs. [24, 45] we further assume that EFT errors of scattering observables at different energies and angles are independent.

The experimental data set, $D=\{y_{\mathrm{exp},i}\}_{i=1}^{N}$, on which we condition the MCMC inference is listed in Table 1 and contains $N=1043$ scattering observables. This is not all the data in the Granada database [65, 66]. We attempted an inference using all $np$ observables with $T_{\mathrm{lab}}<100$ MeV, but a model check revealed that most calibration data were poorly reproduced. For example, low-energy total cross sections deviated significantly from experimental values, which in turn yielded significantly overbound deuteron states for most cutoffs. This is probably due to a misspecified EFT truncation error, leading to unnaturally large values for $\bar{c}$ and thus overestimated EFT errors for low-energy cross sections (that are expected to be reproduced relatively well at LO). There are several suggestions [22] for how to construct more sophisticated EFT error models that would accommodate the incorporation of more data in the likelihood. We postpone such developments until we have incorporated sub-leading orders in MWPC.

III.2 Priors

We treat $\bm{\alpha}$ and $\bar{c}$ as random variables and must therefore place priors on them. For the LECs we employ independent and normally distributed priors with a standard deviation of four times their naturalness estimates: $\frac{4\pi}{f^{2}_{\pi}}\approx 0.1\times 10^{4}\ \mathrm{GeV}^{-2}$ and $\frac{4\pi}{f^{2}_{\pi}\Lambda^{2}_{\chi}}\approx 0.4\times 10^{4}\ \mathrm{GeV}^{-4}$ for the $S$- and $P$-wave LECs, respectively [39, 9, 10]. Similarly, we expect $\bar{c}$ to be of natural size and always positive. As shown in Fig. 4, we employ an inverse gamma distribution

[TABLE]

with parameters $\alpha=3$ and $\beta=4.2$. This distribution has a mode around one, but also a heavy tail that allows for a significant variation.

IV Sampling the posterior at different cutoff values

The parameter inference proceeds in two steps. First, we employ a Bayes linear approach known as history matching [67, 68, 69, 28] to identify relevant domains for the LECs. Second, we employ an affine invariant MCMC algorithm called emcee [70] to numerically draw samples from the posterior pdf. The history-matching domains identified in the first step are used to initialize the MCMC sampling. Furthermore, to avoid detrimental influence of the limit-cycle-like behavior during sampling we select momentum cutoffs based on the analysis in Sec. II.3.

IV.1 History Matching

History matching is an iterative scheme that can be employed to identify and exclude regions of parameter space that produce model outputs inconsistent with observational data, taking relevant uncertainties into account. See, e.g., Ref. [28] for details of the algorithm and an application in nuclear physics. At the end, history matching returns the region(s) of parameter space that could not be ruled out, the so-called non-implausible domain.

In this work, however, we employ history matching as a precursor to a full Bayesian analysis. A similar use of history matching can be found in Ref. [71]. Our aim is then to identify regions of parameter space where we expect to find large contributions to the posterior probability mass. We will be less concerned with the observational property of the data that we utilize in this step (we will employ scattering phase shifts), or with the rigour of our uncertainty estimates, as the sole purpose is to pick promising starting points for the MCMC algorithm. For this reason we refer to excluded samples as inefficient starting points, whereas the final domain is at least non-inefficient. These designations replace the standard implausible and non-implausible labels that are found in the history-matching literature.

We employ an “inefficiency” measure $I_{M}(\bm{\alpha})$ that gauges the performance of the theoretical predictions for a selected subset of data given values for the LECs

[TABLE]

The index, $i$, enumerates the set of experimental observations that is included in the history matching. The measure $I_{M}(\bm{\alpha})$ is defined as the maximum over the observations $\{y_{\mathrm{exp},i}\}$ and is therefore governed by the datum that is reproduced the worst. The experimental and theoretical errors are incorporated via the variances $\sigma^{2}_{\mathrm{exp},i}$ and $\sigma^{2}_{\mathrm{th},i}$, respectively. A second measure, $I_{2M}(\bm{\alpha})$, is sometimes used to construct an additional constraint. It is analogously defined as the second maximum over the observations and obviously fulfills $I_{2M}(\bm{\alpha})\leq I_{M}(\bm{\alpha})$.

We perform two waves of history matching. In the first wave, the scattering phase shifts listed in Table 2 are considered as observational data. Since all LECs act in distinct partial waves, each LEC could be constrained independently within one-dimensional sub-waves. We estimate conservative phase-shift errors $\sigma_{\mathrm{th}}$ from Eq. 17 using fixed $\bar{c}\approx 6$. This choice implies a $30\%$ error for $T_{\mathrm{lab}}<40.6$ MeV (i.e., $p<m_{\pi}$) and $80\%$ for phase shifts with $T_{\mathrm{lab}}\approx 100$ MeV. For each sub-wave we use $10^{4}$ different LEC values in a space-filling Latin hypercube design [72] across a rather wide interval informed by the phase-shift analysis in Sec. II.3 and corresponding to the ranges shown on the $y$-axes in Fig. 5. Samples for which the inefficiency measures are larger than some cutoff values are deemed as poor candidates for initializing the MCMC algorithm. In this analysis we use a sequence of cutoffs, $I_{M}>3.0$ combined with $I_{2M}>2.5$, which is similar to Ref. [69] (therein based on Pukelsheim’s three sigma rule [73]). In the end, the relevant parameter volume is reduced by a large fraction. The ratio of the non-inefficient over initial volume was found to be as small as $10^{-5}$, although this ratio depends rather strongly on the value of the cutoff.

In the second wave, a set of 13 $np$ scattering observables is used to construct the inefficiency measures. In this wave all four LECs are active simultaneously. Here we employ $10^{5}$ samples, using the same space-filling design, in the non-inefficient domain resulting from the first wave. The larger number of samples, and the much reduced volume helps to provide sufficient resolution for detecting narrow domains. The set of observables and corresponding model uncertainties, estimated using the same prescription as in the first wave, can be found in Table 3.

The resulting, relevant ranges for each of the LECs are shown in Fig. 5. The domains are classified, and color-coded, according to if the LEC gives rise to a spurious (unphysical) bound state in the respective channel. The rapid variation of the phase shift by $180^{\circ}$ where a bound state appears—discussed and illustrated in connection to Fig. 3—gives rise to disjoint domains. Therefore, we can anticipate multimodal structures of the LEC posterior pdf.

Depending on the cutoff, we find spurious states in all channels. In the ${}^{1}S_{0}$ channel there is always a region with an unphysical bound state, but it becomes more narrow for increasing values of the cutoff. The physical deuteron state is always present in ${}^{3}S_{1}\mathrm{-}^{3}D_{1}$, and we observe additional (spurious) states in disjoint domains for several cutoffs. For $\Lambda\geq 3000$ MeV there are two disjoint domains that contain at least one spurious bound state each.

While there are no spurious bound states in ${}^{3}P_{0}$ for $\Lambda\leq 500$ MeV, a second domain with lower values of the LEC contains a spurious bound state for $\Lambda=550,600,700$ MeV. For larger cutoffs, $\Lambda\geq 2000$ MeV, there is only one narrow domain and it gives rise to a spurious bound state. In ${}^{3}P_{2}$ we observe larger domains and spurious bound states for $\Lambda\geq 2000$ MeV.

We can understand the $\Lambda$-dependence of the $P$-wave domains by studying Fig. 3 in light of the errors going into the inefficiency measure (21). In doing so we see that for ${}^{3}P_{0}$ and $\Lambda=550-700$ MeV the theoretical and empirical phase shifts match across wide regions of LEC values that qualitatively agree with the non-inefficient domains observed in Fig. 5. A similar argument holds for ${}^{3}P_{2}$ and greater values of $\Lambda$.

Unfortunately, history-matching fails to reduce the initial domain of $C_{{}^{3}P_{0}}$ and $C_{{}^{3}P_{2}}$ for some momentum cutoffs and thus provides limited information in these cases. In the following we perform a Bayesian analysis. The LEC posteriors presented in Section IV.2 are conditioned on observational data (see Table 1) and a more reliable error model. They provide more informative and credible inference compared to the domains shown in Fig. 5.

IV.2 Markov chain Monte Carlo sampling

We sample the posterior pdf in Eq. (8) using the affine-invariant ensamble sampler from the python package emcee [70] conditioned on the experimental data in Table 1. As noted in Sec. IV.1, there are disconnected non-inefficient regions for some LECs at several values of the cutoff. These regions might correspond to multimodal structures of the posteriors. In infinite time, the MCMC walkers will explore all parameter space, but in finite time they are likely trapped in local modes. To handle such convergence challenges we initiate sampling in several of the regions identified using history matching. Ideally, all combinations of such regions should be explored in order to locate the dominant posterior mode(s) with some certainty.

We initialize 50 MCMC walkers using a random subset of points from the non-inefficient domains for each cutoff and let them take 5000 steps each after a burn-in period of 1000 steps. This is typically sufficient for obtaining a good representation of each posterior mode of our five-dimensional posterior pdf. The initialization of walkers in the $P$-waves are chosen in the region of the domain without bound states, shown in orange in Fig. 5, except for the higher cutoffs in ${}^{3}P_{0}$. The same principle is applied in the ${}^{3}S_{1}$ partial wave, where we initialize in the region containing the least amount of spurious states. The posterior turns out to be relatively flat in the direction of $P$-wave LECs. Thus, the MCMC sampling is less sensitive to where the walkers are initialized in the rather large start domains. Moving forward, we focus on the possible multimodalities originating from the ${}^{1}S_{0}$ channel where we always found a parameter region giving rise to a shallow and unphysical bound state. The initialization of MCMC walkers in the ${}^{1}S_{0}$ channel is chosen either in the unbound (orange) or bound (blue) region. Proceeding like this for selected cutoff values between 400 MeV and 4000 MeV we indeed find multimodal structures in the posterior pdf, as summarized in Table 6 in Appendix C.

To analyze the probability mass gathered in each mode we calculate the model evidence $\mathcal{Z}_{K}$ for all domains $K$ enclosing a mode. We do this using the Laplace approximation, which is accurate if the probability distribution can be locally approximated with a multivariate Gaussian, i.e.,

[TABLE]

where $\Sigma^{-1}$ is the Hessian matrix of $\mathrm{pr}\left(\bm{\theta}|D,I\right)$ evaluated at the mode $\bm{\theta}^{*}$. Each mode listed in Table 6 in Appendix C is indeed found to be well approximated with a Gaussian. The numerical evaluation of Eq. 22 might be sensitive to the convergence of the MCMC chains. However, this sensitivity mainly impacts the evaluation of the determinant of the inverse Hessian which turns out to be a less significant contribution to the evidence compared to the density at the MAP point [74]. The $\ln(\mathcal{Z}_{K})$ values for all investigated modes are summarized in Table 6 of Appendix C. Fortunately, a well-defined “best” mode with highest evidence can be identified at all cutoffs except for $\Lambda=700$ and $\Lambda=4000$ MeV for which there are two $P$-wave modes with very similar evidences. All the significant modes identified in the evidence analysis are retained, and none of them contain an unphysical bound state in the ${}^{1}S_{0}$ channel.

To ensure that we obtain a high-quality representation of the relevant posteriors we perform the final sampling using 50 walkers taking 5000 steps each initialized in the vicinity of the largest $\ln(\mathcal{Z}_{K})$ mode for each cutoff. This lays the foundation for accurate inferences of the ppds presented in Sec. V. We use an autocorrelation analysis [24] and visual inspection of the traces to confirm MCMC convergence. At the cutoffs $\Lambda=700$ MeV and $\Lambda=4000$ MeV the two $P$-wave modes with similar evidences remain. The resulting values for the model evidences and MAP locations at several cutoffs are shown in Table 4.

The model evidences are of similar magnitude for all cutoffs. In fact, the biggest difference, on a logarithmic scale, is about $40$, i.e., the maximal ratio of evidences is $e^{40}$. The same difference can, for example, be obtained by shifting all predictions by about $5\%$, which is quite small compared to the EFT truncation error.

In Figs. 6a and 6b, we show the marginal posterior pdfs for the cutoff values $\Lambda=450$ MeV and $\Lambda=700$ MeV. The posterior pdfs for the remaining cutoff values are shown in Appendix C. The posteriors are consistent with our naturalness assumptions for both the LECs and the scale of the EFT truncation error as no tensions with our priors are observed.

The running of the inferred LECs can now be monitored. The cutoff dependence of the marginal posteriors for the LECs is shown in Fig. 7. We conclude that $\tilde{C}_{{}^{1}S_{0}}(\Lambda)$ and $\tilde{C}_{{}^{3}S_{1}}(\Lambda)$ have a similar behavior as the running of those couplings determined from the phase shift fits in Sec. II.3 and Fig. 1. This finding can likely be attributed to the fact that these LECs are quite well constrained by the employed low-energy $np$ scattering cross sections. However, the $P$-wave LECs $C_{{}^{3}P_{0}}$ and $C_{{}^{3}P_{2}}$ are not that well constrained by the $np$ data used in the inference and run differently with $\Lambda$ in the Bayesian analysis.

The inference of $\bar{c}$ is quite interesting in its own right. In Fig. 7 we show its marginal posterior and see that it is rather insensitive to cutoff variation, in particular for $\Lambda>700$. The breakdown scale $\Lambda_{\chi}$, $\bar{c}$ and the low-energy scale, $Q$, are connected through Eq. 17 and their posteriors have previously been investigated in WPC in Refs. [21, 23].

V Posterior predictive distributions

In this section we present ppds for scattering phase shifts as well as selected $np$ observables. The ppd for an observable $y$ is a marginalization of the model predictions over the posterior of the relevant parameters

[TABLE]

When sampling the model prediction, $\mathrm{pr}\left(y|\bm{\theta},D,I\right)$, we can choose to include the EFT truncation error, in which case it takes the form

[TABLE]

or we can exclude it, which corresponds to just propagating the parametric LEC uncertainty, by using a delta distribution

[TABLE]

The integral in Eq. 23 can be straightforwardly evaluated using the samples from the MCMC chains obtained in Sec. IV. In the following, ppds including (excluding) the EFT error are colored blue (purple).

V.1 Phase shifts

The Bayesian inference in this work differs from the point estimates obtained from a standard phase shift optimization. In Fig. 8 we compare the ppds from the present analysis with the optimized phase shifts from Yang et al. [5] for $\Lambda=450$ MeV and the Nijmegen partial-wave analysis [62]. The EFT error is not included in the ppd. One finds that the parametric uncertainty of the phase shifts, stemming from the LEC posterior, is rather small and the $95\%$ credible intervals shown in the figure are mostly visible in the ${}^{3}P_{2}$ partial wave. We find good agreement with the results by Yang et al. [5] except for the $P$-wave phase shifts which are more repulsive in the Bayesian analysis. Apparently, the $P$-wave LECs receive significant contributions when conditioning on scattering data. At $\Lambda=450$ MeV, The $P$-wave phase shifts ${}^{1}P_{1}$ and ${}^{3}P_{1}$, which are not part of the contact potential, are more attractive than the ones from the Nijmegen analysis, see Fig. 9. Note that the LO potential only includes low partial waves, for which the OPE is not perturbative. To prevent overfitting to higher order effects, which are better described by high partial waves, either the LO potential or the model for the LO EFT error likely needs to be revisited. We also foresee that the LO LECs will receive corrections at higher orders and this could significantly change the scattering amplitudes and the pole structure in a given channel. The $np$ phase shifts for the MAP estimate of the LECs for $\Lambda=450$ MeV are tabulated in Table 7 in Appendix D.

V.2 Observables in the $np$ sector

In Fig. 10 we compare experimental data for several $np$ scattering observables and the corresponding ppds, including the EFT truncation error, at LO and using the cutoff $\Lambda=450$ MeV. Total cross sections, $\sigma$, and differential cross sections, ${d\sigma}/{d\Omega}$, are reproduced quite well, although the ppd error bands are rather wide. Using Eq. 17 and $\bar{c}\approx 2.8$ from Fig. 6a we estimate that the EFT truncation error for $p\lesssim m_{\pi}$ is about $15\%$. The ppd for the spin-polarization observable $P_{b}$, which was not included in the MCMC sampling, reproduces experimental data rather poorly. This result is particularly striking in the ppd at $T_{\mathrm{lab}}=25.0$ MeV. The situation improves at higher energies, which is somewhat surprising for a low-energy EFT. As such, conditioning on low-energy $P_{b}$ data could improve the model. However, expanding the dataset in Table 1 requires an improved model for the EFT truncation error  [22] and likely access to order-by-order calculations beyond LO to guide the inference of $\bar{c}$.

In Fig. 11 we show the ppds for $\sigma$ and $P_{b}$ as a function of $\Lambda$, including the EFT error, for a handful of values of $T_{\mathrm{lab}}$ and $\theta_{\mathrm{c.m.}}$. Except for a slight variation of $P_{b}$ for $\Lambda<700$, no significant cutoff variation is observed. It can also be seen that the experimental value is well reproduced for $\sigma$ and more poorly for $P_{b}$. The LO EFT error for $\sigma$ appears to be on the conservative side and might shrink as we learn it better.

The ppd, without the EFT error, for the deuteron ground-state energy is shown in Fig. 12 and reveals an approximate 0.5 MeV underbinding, but no significant cutoff dependence. For some of the higher cutoffs additional spurious and deeply bound states appear at $E\approx-1$ GeV. Employing the same error model as for the scattering observables one can attempt to make an uncertainty estimate for the deuteron ground-state energy. Using the binding momentum $p_{gs}=\sqrt{m_{N}(-E_{gs})}\approx 50\ \mathrm{MeV}<m_{\pi}$ as a proxy for the relevant soft scale we arrive at an estimated EFT error of around $15\%$ of the predicted value. Including this error makes the lower end of the ppds touch the experimental value.

The ppds for the energies of the unphysical bound states in the channels ${}^{3}P_{0}$ and ${}^{3}P_{2}\mathrm{-}^{3}F_{2}$ are shown in Fig. 13, except for a very deep ($E\approx-1$ TeV) state in ${}^{3}P_{0}$ for $\Lambda=4000$ MeV. Some of the energies are around $25$ MeV, which is well within the validity of the EFT. Note also that the ppds in Fig. 13 show a strong dependence on the cutoff, which is generally the case for spurious states, which was also seen in Sec. II.3. In the ${}^{3}P_{2}\mathrm{-}^{3}F_{2}$ channel for the cutoff $\Lambda=700$ the deeply bound state at $E\approx-1$ GeV is produced by the very small tail of negative $C_{{}^{3}P_{2}}$-values shown in the marginal pdf in Fig. 7. We see spurious states in all singular channels for $\Lambda\geq 2000$ MeV.

VI Conclusions and outlook

In this work we have presented a detailed Bayesian study of a LO $\chi$EFT potential in MWPC. We have performed robust inferences of the LECs and the scale $\bar{c}$ of the EFT truncation error using history matching and MCMC sampling for several values of the momentum cutoff in the range from 400 to 4000 MeV. Below, we summarize the main takeaways.

• •

Multimodal LEC posteriors are induced via the existence of bound states, virtual states and resonances that move around near threshold for varying values of the LECs and cutoff $\Lambda$. This makes inference challenging. In particular, we could always find modes where the ${}^{1}S_{0}$ state was bound at an energy of approximately 1 MeV.

• •

We computed the model evidence for each mode in the Laplace approximation. At each cutoff we initialized a final round of MCMC sampling around the mode(s) with the largest evidence (see Table 4). Modes harboring a bound ${}^{1}S_{0}$ state were always found to be insignificant. However, pairs of modes in the ${}^{3}P_{2}\mathrm{-}^{3}F_{2}$ channel for $\Lambda=700$ MeV and in the ${}^{3}P_{0}$ channel for $\Lambda=4000$ MeV had similar evidences and could not be discerned. We expect that the poor constraints on the LECs might be improved with a better error model and with the inclusion of more types of scattering observables in the calibration data.

• •

Our inference shows that $\bar{c}\approx 2.8$ across all cutoffs investigated. Since the main posterior modes have similar evidences (see Table 4) this implies that the description of the low-energy scattering data is essentially equivalent for the different momentum cutoffs in the range 400 to 4000 MeV.

• •

Conditioning the LEC inference on scattering data leads to excessive repulsion in the ${}^{3}P_{0}$ and ${}^{3}P_{2}$ partial waves compared to the Nijmegen analysis [62] and to the phase-shift optimization performed by Yang et al. [5]. Despite the poor accuracy of $P$-wave phase shifts (see Fig. 8), we find that scattering observables are reproduced surprisingly well (see Figs. 10, 11 and 12) for a LO potential. We conclude that it might be problematic to use phase shifts for calibration of LO potentials in MWPC. These potentials, acting in the first few partial waves only, need to compensate for excluded higher-order contributions to reproduce observables. Consequently, the inference becomes crucially dependent on the specification of the EFT error model, of which we have limited information. This highlights the need for further studies at higher orders.

• •

We conclude that the ppds for various bound-state and scattering observables in the $np$ sector exhibit no significant cutoff dependence.

It remains to be seen what happens as sub-leading orders are included perturbatively and when the model for the EFT truncation error is refined. The inferred LO potentials, with broad distributions for the $P$-wave LECs, will likely lead to significant predictive uncertainties in nuclei. Therefore, the ppds for ground-state energies and radii should also be quantified to assess the quality and physical relevance of existing RG-invariant $\chi$EFT formulations.

Acknowledgements.

O.T. thanks C.-J. Yang for helpful discussions regarding the potential and conventions, as well as providing more detailed data from Ref. [5]. This work was supported by the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (Grant Agreement No. 758027), the Swedish Research Council (Grants No. 2017-04234, No. 2020-05127 and No. 2021-04507). The computations were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) partially funded by the Swedish Research Council through Grant Agreement No. 2018-05973.

Appendix A Neutron-proton scattering

The scattering geometry in relative-momentum coordinates is illustrated in Fig. 14.

For a stationary proton and an incoming neutron with kinetic energy $T_{\mathrm{lab}}$ in the laboratory frame of reference, the modulus of the relative momentum, $p$, is given by

[TABLE]

We partial-wave decompose the potential $V$ in Eq. 2 using the helicity formalism introduced in Ref. [53], and employ partial wave states $\ket{lm_{l},p}$ with the normalization

[TABLE]

where $Y^{l}_{m_{l}}(\hat{p})$ is the spherical harmonics, including the Condon-Shortley phase, normalized as

[TABLE]

We are interested in elastic scattering and express the $T$-matrix as a real-valued $R$-matrix defined as

[TABLE]

where $H_{0}$ is the free Hamiltonian and $E=p^{2}/m_{N}$ [75]. Consequently, we solve the Lippmann-Schwinger equation given by

[TABLE]

where $\mathcal{P}$ denotes the principal value. The momentum-space integral is discretized as described in Refs. [55, 76] using a Gauss-Legendre grid with $N_{p}=60$ points in the interval $[0,\Lambda+300]$ MeV, i.e. on the interval where the potential, and hence the part of the integral containing the potential is non-zero. The advantage of using a grid that does not extend to infinity is that the part of the integral that implements the principal value, and is outside the support of the potential, can be treated analytically, see, e.g., Ref. [77].

Following Refs. [75, 78] we obtain scattering phase shifts in the Blatt and Bidenharn convention from on-shell $R$-matrix elements in uncoupled channels via

[TABLE]

and for coupled channels via

[TABLE]

where $R^{sj}_{\pm\pm}\equiv R^{sj}_{l^{\prime}=j\pm 1,l=j\pm 1}(p,p)$. The more commonly employed Stapp- (bar-) phase shifts, which we use in this work, are given by the solution to the following equations [75]

[TABLE]

Note that numerical instabilities can occur when the mixing angle ($\bar{\epsilon}^{sj}$) in the Stapp-convention changes sign.

The spin-scattering matrix, $M$, is given by [79, 75, 80]

[TABLE]

The angles $\theta_{\mathrm{c.m.}}\in[0,\pi]$ and $\phi\in[0,2\pi]$ are the polar and azimuthal scattering angles, respectively. We set $\phi=0$ by cylindrical symmetry. Numerical instabilities in Eq. 34 can occur at $\theta_{\mathrm{c.m.}}=90^{\circ}$, but can be avoided by adding a small nugget.

Equipped with $M$, spin-scattering observables can be computed from traces in spin space. The spin averaged differential cross section is obtained as

[TABLE]

Other observables can be calculated in a similar fashion, or by decomposing the $M$-matrix into, e.g., Saclay-amplitudes as described in Ref. [56].

Appendix B Appearance of bound states

Table 5 lists the threshold LEC values for which $np$ bound states appear at different cutoff values, $\Lambda$, for the potential discussed in Eq. 5.

Appendix C Marginal posterior distributions

Table 6 lists the MAP LEC values for the modes found during the initial MCMC. Figures 15 and 16 show the marginal posterior pdfs for the parameters $\bm{\theta}$ at cutoff values $\Lambda=400,500,550,600,2000,3000,4000$ MeV.

Appendix D Phase shift benchmark

Table 7 lists selected phase shifts for the LO EFT potential defined in Sec. II with $\Lambda=450$ MeV and the MAP LEC values $\left(\tilde{C}_{{}^{1}S_{0}},\tilde{C}_{{}^{3}S_{1}},C_{{}^{3}P_{0}},C_{{}^{3}P_{2}}\right)$ = (-0.1128, -0.0964, 5.1151, 0.5118) in the units defined below Eq. 10.