Bayesian truncation errors in chiral effective field theory: nucleon-nucleon observables

TL;DR

This paper develops and applies a Bayesian model to quantify truncation errors in chiral effective field theory predictions for nucleon-nucleon observables, enabling robust statistical analysis and validation of EFT assumptions.

Contribution

It extends previous Bayesian truncation error models to a broader set of observables and parameters, providing a general statistical framework for EFT error quantification.

Findings

Bayesian model successfully quantifies EFT truncation errors.

Posterior distributions for EFT breakdown scale are consistent with theoretical expectations.

Model validation confirms applicability to regulator artifact-free potentials.

Abstract

Chiral effective field theory (EFT) predictions are necessarily truncated at some order in the EFT expansion, which induces an error that must be quantified for robust statistical comparisons to experiment. In previous work, a Bayesian model for truncation errors of perturbative expansions was adapted to EFTs. The model yields posterior probability distribution functions (pdfs) for these errors based on expectations of naturalness encoded in Bayesian priors and the observed order-by-order convergence pattern of the EFT. A first application was made to chiral EFT for neutron-proton scattering using the semi-local potentials of Epelbaum, Krebs, and Mei{\ss}ner (EKM). Here we extend this application to consider a larger set of regulator parameters, energies, and observables as a general example of a statistical approach to truncation errors. The Bayesian approach allows for statistical validations of the assumptions and enables the calculation of posterior pdfs for the EFT breakdown scale. The statistical model is validated for EKM potentials whose convergence behavior is not distorted by regulator artifacts. For these cases, the posterior for the breakdown scale is consistent with EKM assumptions.

Figures

Tables

Table 1. Table 1: Candidates for prior pdfs  Furnstahl  et al. ( 2015b ) .

Set	$\mathrm{pr}\left(c_n|\overline{c}\right)$	$\mathrm{pr}\left(\overline{c}\right)$
A	$\frac{1}{2\overline{c}}\theta \left(\overline{c}-\left|c_n\right|\right)$
B	$\frac{1}{2\overline{c}}\theta \left(\overline{c}-\left|c_n\right|\right)$	$\frac{1}{\sqrt{2\pi }\overline{c}\sigma }{e}^{-{\left(\ln \overline{c}\right)}^2/2{\sigma }^2}$
C	$\frac{1}{\sqrt{2\pi }\overline{c}}{e}^{-c_n^2/2{\overline{c}}^2}$

Table 2. Table 2: Comparison of notations for selected NN scattering observables.

Name	$X_{pqik}$	Others
Differential Cross Section	$I_{0000}$	$\sigma$, $d\sigma /d\Omega$
Vector Analyzing Power	$A_{00n0}$	$A_y$, $P_b$
Polarization Transfer	$D_{s^{\prime }0k0}$	$A$
	$D_{n0n0}$	$D$
Spin Correlation Parameters	$A_{00ss}$	$A_{xx}$
	$A_{00nn}$	$A_{yy}$

Equations

$$X=X_{\text{ref}}\sum _{n=0}^{\infty }{c}_n{Q}^n.$$

$${\Delta }_k\equiv \sum _{n=k+1}^{\infty }{c}_n{Q}^n.$$

$${\mathbf{c}}_k=(c_2,c_3,\dots ,c_k).$$

$$p={\int }_{-d_k^{(p)}}^{d_k^{(p)}}\mathrm{d}\Delta {\mathrm{pr}}_h(\Delta |{\mathbf{c}}_k).$$

$$⟨{\mathbf{c}}_k^2⟩=(k-1){\bar{c}}^2,⟨c_n⟩=0.$$

$${\mathrm{pr}}_h(\Delta |{\mathbf{c}}_k)=\frac{{\int }_0^{\infty }\mathrm{d}\bar{c}{\mathrm{pr}}_h(\Delta |\bar{c})\mathrm{pr}(\bar{c})\prod _{n=2}^k\mathrm{pr}(c_n|\bar{c})}{{\int }_0^{\infty }\mathrm{d}\bar{c}\mathrm{pr}(\bar{c})\prod _{n=2}^k\mathrm{pr}(c_n|\bar{c})},$$

$${\mathrm{pr}}_h(\Delta |\bar{c})\equiv \left[\prod _{i=k+1}^{k+h}{\int }_{-\infty }^{\infty }\mathrm{d}{c}_i\mathrm{pr}(c_i|\bar{c})\right]\delta \left(\Delta -{\Delta }_k^{(h)}\right).$$

$${\mathrm{pr}}_1(\Delta |{\mathbf{c}}_k)=\frac{{\int }_0^{\infty }\mathrm{d}\bar{c}\mathrm{pr}(c_{k+1}|\bar{c})\mathrm{pr}(\bar{c})\prod _{n=2}^k\mathrm{pr}(c_n|\bar{c})}{Q^{k+1}{\int }_0^{\infty }\mathrm{d}\bar{c}\mathrm{pr}(\bar{c})\prod _{n=2}^k\mathrm{pr}(c_n|\bar{c})},$$

$$Q=\frac{\text{max}\{p,m_{\pi }\}}{{\Lambda }_b}.$$

$$Q_{\text{interp}}(p)=\frac{m_{\pi }^n+p^n}{m_{\pi }^{n-1}+p^{n-1}}\frac{1}{{\Lambda }_b},$$

$$X_{\text{res}}\equiv X-X_{\text{NPWA}}.$$

$$\frac{\mathrm{d}\sigma }{\mathrm{d}\Omega }{X}_{pqik}=X_{\text{ref}}\sum _{n=0}^{\infty }{c}_n{Q}^n,$$

$$X_{pqik}=\sum _{n=0}^{\infty }{c}_n{Q}^n,$$

$$\mathrm{pr}(n|p,N)=\frac{N!}{n!(N-n)!}{p}^n(1-p{)}^{N-n}.$$

$$\mathrm{pr}(p|n,N)\propto \mathrm{pr}(n|p,N)\mathrm{pr}(p)\propto \mathrm{pr}(n|p,N),$$

$$X=X_{\text{ref}}\sum _{n=0}^{\infty }(c_n{\lambda }^n)\times {\left(\frac{Q}{\lambda }\right)}^n.$$

$$X\equiv X_{\text{ref}}\sum _{n=0}^{\infty }{b}_n{p}^n,$$

$$c_n={\Lambda }_b^n{b}_n.$$

$$\mathrm{pr}({\Lambda }_b|{\mathbf{b}}_k({\mathbf{\alpha }}_1),\dots ,{\mathbf{b}}_k({\mathbf{\alpha }}_m)),$$

$$\begin{array}{rl}\mathrm{pr}({\Lambda }_b|{\mathbf{b}}_k(& {\mathbf{\alpha }}_1),\dots ,{\mathbf{b}}_k({\mathbf{\alpha }}_m))\\ & =\frac{\mathrm{pr}({\mathbf{b}}_k({\mathbf{\alpha }}_1),\dots ,{\mathbf{b}}_k({\mathbf{\alpha }}_m)|{\Lambda }_b)\mathrm{pr}({\Lambda }_b)}{\mathrm{pr}({\mathbf{b}}_k({\mathbf{\alpha }}_1),\dots ,{\mathbf{b}}_k({\mathbf{\alpha }}_m))}.\end{array}$$

$$\mathrm{pr}({\mathbf{b}}_k({\mathbf{\alpha }}_1),\dots ,{\mathbf{b}}_k({\mathbf{\alpha }}_m)|{\Lambda }_b)=\prod _{i=1}^m\mathrm{pr}({\mathbf{b}}_k({\mathbf{\alpha }}_i)|{\Lambda }_b).$$

$$\begin{array}{rl}\mathrm{pr}({\Lambda }_b|{\mathbf{b}}_k(& {\mathbf{\alpha }}_1),\dots ,{\mathbf{b}}_k({\mathbf{\alpha }}_m))\\ & =\frac{\mathrm{pr}({\Lambda }_b){\prod }_{i=1}^m\mathrm{pr}({\mathbf{b}}_k({\mathbf{\alpha }}_i)|{\Lambda }_b)}{\mathrm{pr}({\mathbf{b}}_k({\mathbf{\alpha }}_1),\dots ,{\mathbf{b}}_k({\mathbf{\alpha }}_m))}.\end{array}$$

$$\mathrm{pr}({\mathbf{b}}_k|{\Lambda }_b)=\int \mathrm{d}{\mathbf{c}}_k\mathrm{pr}({\mathbf{b}}_k|{\mathbf{c}}_k,{\Lambda }_b)\mathrm{pr}({\mathbf{c}}_k|{\Lambda }_b).$$

$$\mathrm{pr}({\mathbf{b}}_k|{\Lambda }_b)=\int \mathrm{d}\bar{c}\mathrm{d}{\mathbf{c}}_k\mathrm{pr}({\mathbf{b}}_k|{\mathbf{c}}_k,{\Lambda }_b)\mathrm{pr}({\mathbf{c}}_k|\bar{c},{\Lambda }_b)\mathrm{pr}(\bar{c}|{\Lambda }_b).$$

$$\mathrm{pr}({\mathbf{c}}_k|\bar{c},{\Lambda }_b)=\prod _{n=2}^k\mathrm{pr}(c_n|\bar{c}).$$

$$\mathrm{pr}({\mathbf{b}}_k|{\mathbf{c}}_k,{\Lambda }_b)$$

$$\mathrm{pr}$$

$$\mathrm{pr}(b_n|c_n,{\Lambda }_b)=\delta \left(b_n-\frac{c_n}{{\Lambda }_b^n}\right),$$

$$\mathrm{pr}$$

Full text

I Introduction

The scope of ab initio nuclear structure and reactions has increased dramatically due to recent advances in many-body methods Barrett et al. (2013); Hagen et al. (2014); Carlson et al. (2015); Hergert et al. (2016); Barbieri and Carbone (2016); Lee (2016), continued growth in computational power, and new developments in chiral effective field theory (EFT) Epelbaum et al. (2009); Machleidt and Entem (2011); Gezerlis et al. (2014); Epelbaum et al. (2015a); Piarulli et al. (2015); Entem et al. (2017). To properly judge the successes and predictive power of ab initio nuclear theory, however, it is necessary that theory errors be understood. Thus, quantifying the theoretical uncertainties of nuclear calculations has now become a critical task for confronting experiment and theory and for extrapolating to unmeasured phenomena Dobaczewski et al. (2014); Ireland and Nazarewicz (2015); Carlsson et al. (2016).

Uncertainties in chiral EFT predictions arise from three sources Furnstahl et al. (2015a): uncertainty in the input data to which the EFT parameters are fit, errors in the Hamiltonian, and numerical approximations. Here we focus on quantifying the Hamiltonian truncation error as part of the larger BUQEYE program Furnstahl et al. (2015a) of quantifying all uncertainties for EFT predictions. Despite the promise of systematic expansions, uncertainties from truncation have been difficult to estimate and, when provided, generally lack a well-defined statistical interpretation.

In Ref. Furnstahl et al. (2015b), a Bayesian model for truncation errors originally applied to perturbative expansions in quantum chromodynamics Cacciari and Houdeau (2011); Bagnaschi and Jenniches (2014) was adapted to EFTs. The generic assumption is that the EFT provides us with a dimensionless expansion parameter $Q$, which is a ratio of scales, and an associated expansion for quantities $X$ (usually observables):

[TABLE]

Here, $X_{\text{ref}}$ is the natural size of $X$, which could be the leading-order estimate $X_{0}$, and the $c_{n}$s are dimensionless coefficients. For chiral EFT $c_{1}$ is zero by symmetry, and we have a double expansion in $Q=\{p,m_{\pi}\}/\Lambda_{b}$, where $p$ is the relative momentum of two scattering nucleons, $m_{\pi}$ is the pion mass, and $\Lambda_{b}$ is the EFT breakdown scale. The goal is to estimate the error incurred in the observable by truncating the expansion at order $k$. Note that this does not exclude an asymptotic expansion, but assumes that we truncate while the result is still improving.

In some cases the expansion in Eq. (1) may follow directly from a perturbative EFT expansion of a Lagrangian, i.e., through a sum of Feynman diagrams with powers of $Q$ coming from a simple power-law dependence on momentum or a mass (such as the pion mass in a chiral perturbation theory expansion). There will also be implicit $Q$ dependence, often in the form of logarithms, which vary much more slowly. But in other cases, such as EFT for more than one nucleon, the calculations are nonperturbative, and the dependence on momentum or energy will be complicated and nonlinear in general. Nevertheless, if the EFT is working we expect the calculation of $X$ to improve systematically as we go to higher orders.

Equation (1) can be interpreted as a summary of that expected systematic improvement. Namely, that the correction term with each successive order is on average a factor $Q$ smaller than the previous order. For this to be the case, we need the $c_{n}$ coefficients to be roughly the same size. Because the coefficients are unknown a priori, we treat them as drawn from a random distribution with a characteristic size. This is a realization of the underlying assumption that the naturalness of the low-energy constants (LECs) in the EFT Lagrangian propagates to the expansion for any observable. We have no general proof of this assumption, so we aim to validate it in each application.

In Ref. Furnstahl et al. (2015b) we made a first pass at formalizing and testing the assumptions behind the expansion in Eq. (1), building on an analogous Bayesian analysis applied to perturbative QCD calculations Cacciari and Houdeau (2011); Bagnaschi and Jenniches (2014). We considered various priors for the $c_{n}$s, made an application to a small subset of results from Epelbaum, Krebs, and Meißner (here EKM) for neutron-proton ($np$) scattering cross sections using their new semi-local potentials Epelbaum et al. (2015b, a), and tested the consistency of assumed expansion parameters, which are associated with the expected breakdown scales of the EFT implementation. Here we revisit the EKM application to further test and generalize those investigations, which will set the stage for extending our model of EFT truncation errors.

We seek to address the following questions:

Coefficients $c_{0}$–$c_{5}$ of the total cross section $\sigma$ given at four energies in Refs. Epelbaum et al. (2015b, a) were examined in Ref. Furnstahl et al. (2015b). Can we validate a posteriori our assumption that the observable coefficients follow some bounded random distribution about zero for all energies? 2. 2.

The truncation error model of Ref. Furnstahl et al. (2015b) has not yet been applied to other nucleon-nucleon (NN) observables calculated in chiral EFT, such as the differential cross section and various spin observables. How do the coefficient patterns compare for different NN scattering observables, considered both as functions of energy and scattering angle? Are the naturalness assumptions validated for these observables and for all values of the EKM regulator parameter $R$? 3. 3.

An appropriately assigned $(100*p)\%$ error band should capture the true value of an observable $(100*p)\%$ of the time. How can we utilize known order-by-order results to verify that the error band prescriptions work as advertised? Can information from different observables be treated as independent, and if not, how can we account for their relationships in our analysis? Is there a well-defined “correlation length” in energy or scattering angle beyond which expansion coefficients may be treated as independent of one another? 4. 4.

The identification of the expansion parameter $Q$, which in turn is based on identifying the scale $\Lambda_{b}$, is a key element in determining the convergence pattern. Is it consistent to take $\Lambda_{b}$ to be the same scale for every observable? To what extent can we extract $\Lambda_{b}$, given order-by-order expansions and our naturalness assumptions encoded in a Bayesian model?

In the present work we make progress on all these questions.

In Sec. II, we summarize and extend the relevant formulas from Ref. Furnstahl et al. (2015b). We refer the reader to that article for background on the use of Bayesian statistics in this context, derivations of the formulas we summarize here, and more general references. New results and analysis for the total cross section are given in Sec. III, and other observables are considered in Sec. IV. In Sec. V we perform Bayesian model checking Gelman et al. (2013) by applying a consistency check used in Ref. Furnstahl et al. (2015b) and by calculating posterior probability distribution functions (pdfs) for $\Lambda_{b}$. Section VI has our summary and outlook. For completeness and convenience, we show explicit formulas regarding the Bayesian model for truncation errors in Appendix A and summarize the notation and formulas used for NN observables in Appendix B. The Supplemental Material sup displays extra figures and data that helped inform our conclusions.

II Formulas

If the EFT expansion in Eq. (1) is truncated at order $k$, then the error induced is $X_{\text{ref}}\Delta_{k}$, where the scaled, dimensionless parameter that determines the truncation error is

[TABLE]

Generally it is only practical to approximate $\Delta_{k}^{(h)}$, the error due to the first $h$ omitted higher-order terms. For sufficiently small values of $Q$, the first omitted term $\Delta_{k}^{(1)}=c_{k+1}Q^{k+1}$ is a good estimate for $\Delta_{k}$, but we do not assume this in general.

We use the notation $\operatorname{pr}(x|I)$ to denote the probability density of $x$ given information $I$. Our pdf of interest is $\operatorname{pr}_{h}(\Delta|\mathbf{c}_{k})$: the probability distribution for $\Delta_{k}$ given the vector of relevant lower-order coefficients $\mathbf{c}_{k}$ that have been calculated, assuming that only $h$ higher-order terms contribute to the error and that $\Lambda_{b}$ is to be given from other considerations. The pdf $\operatorname{pr}_{h}(\Delta|\mathbf{c}_{k})$ is normalized in terms of the dummy variable $\Delta$, which is implied to be an estimate of $\Delta_{k}$ contingent on lower-order coefficients $\mathbf{c}_{k}$.

In contrast to Ref. Furnstahl et al. (2015b), here $c_{0}\notin\mathbf{c}_{k}$ because it does not provide insight into the convergence pattern of the observables; rather, the LO calculation provides scaling information. Again, $c_{1}\notin\mathbf{c}_{k}$ because $c_{1}=0$ in chiral EFT. Thus the relevant lower-order coefficients in the determination of $\Delta_{k}$ in chiral EFT are

[TABLE]

In the Bayesian framework, the posterior $\operatorname{pr}_{h}(\Delta|\mathbf{c}_{k})$ contains the complete information we claim to have about the dimensionless residual $\Delta_{k}$. In general, a posterior pdf can have complex structures such as multiple modes, heavy tails, large skewness, etc. Here we can capture most of the information with a small number of degree-of-belief (DoB) intervals.111These are also called “credibility” or “credible” intervals, or “Bayesian confidence intervals”. We use the highest posterior density (HPD) definition of DoB, which is the shortest interval that contains $(100*p)\%$ of the area Kruschke (2011); Gelman et al. (2013); Liu et al. (2015). This ensures that the probability density within the DoB is never lower than the density outside. The HPD definition is particularly well suited for skewed posteriors, as we will encounter in Sec. V. Because the $\operatorname{pr}_{h}(\Delta|\mathbf{c}_{k})$ that we consider here are unimodal and symmetric about $\Delta=0$, finding the DoB interval reduces to the inversion problem for $d_{k}^{(p)}$, where

[TABLE]

Hence, one believes with $(100*p)\%$ certainty the true value of the observable $X$ lies within $\pm X_{\text{ref}}\,d_{k}^{(p)}$ of the NkLO prediction. In general, Eq. (4) must be inverted numerically, but simplified results for certain priors and approximations (e.g., that the first omitted term dominates) are possible Furnstahl et al. (2015b).

In Ref. Furnstahl et al. (2015b), a statistical model for $\Delta_{k}$ in terms of the order-by-order coefficients of the EFT expansion was developed. It was assumed that naturalness could be implemented by treating the $c_{n}$s as random variables drawn from a shared distribution centered at zero with a characteristic size or upper bound $\bar{c}$. The coefficients at each order are treated as independent of one another—the value of $\bar{c}$ is the only way that information propagates between orders. These relationships can be encapsulated in a Bayesian network Ben-Gal (2007); Puga et al. (2015), as shown in Fig. 1. The nodes of the graph are random variables and the arrows denote causal relationships between them.

While the topology of Fig. 1 outlines the logic of our model, prescriptions in the form of priors $\operatorname{pr}(c_{n}|\bar{c})$ and $\operatorname{pr}(\bar{c})$ must be given to make quantitative statistical inferences of $\Delta_{k}$. When all we know is that there is an upper bound to the coefficients, an application of maximum entropy Gull (1998); Sivia and Skilling (2006) dictates that the least-informative distribution $\operatorname{pr}(c_{n}|\bar{c})$ is uniform for $|c_{n}|<\bar{c}$ and zero otherwise. Such uniformity is additionally appealing because it can lead to simple, analytic results. This uniform prior was the initial choice of Ref. Cacciari and Houdeau (2011). We employ it in priors we denote as “set A” and “set B” (see Table 1). The analogous prior of “set C” in Table 1 corresponds to the ensemble naturalness assumption of Ref. Schindler and Phillips (2009). This Gaussian prior follows from the maximum-entropy principle assuming knowledge of testable information on the mean and standard deviation of the $c_{n}$s Schindler and Phillips (2009):

[TABLE]

In addition we require a prior for $\bar{c}$: $\operatorname{pr}(\bar{c})$. Sets A and C of Table 1 use a log-uniform prior for $\bar{c}$ to reflect unbiased expectations regarding the scale of $\bar{c}$ Jeffreys (1939) (this was the choice in Ref. Cacciari and Houdeau (2011) and Ref. Schindler and Phillips (2009)). Such a prior cannot be normalized for $\bar{c}$ in $(0,\infty)$ and is therefore termed an “improper prior”. Limiting the range of $\bar{c}$ through the use of $\theta$ functions permits an examination of the otherwise ill-defined limiting behavior. When marginalizing (i.e. integrating) over $\bar{c}$, we can express complete ignorance of the scale of $\bar{c}$ by considering the limit of infinite range (${\rm A}_{\epsilon}$ or ${\rm C}_{\epsilon}$, see Furnstahl et al. (2015b)), or render the prior more informative through the use of a finite range $[a,b]$ (${\rm A}_{a\text{-}b}$ or ${\rm C}_{a\text{-}b}$). Alternatively, set B employs a log-normal distribution about zero Bagnaschi et al. (2015); Bagnaschi and Jenniches (2014), which sets the scale of $\bar{c}$ with the hyperparameter $\sigma$.

The general result for $\operatorname{pr}_{h}(\Delta|\mathbf{c}_{k})$ implied by Fig. 1 was derived as222We have corrected and simplified here the corresponding equation from Ref. Furnstahl et al. (2015b).

[TABLE]

where

[TABLE]

If we assume that the first omitted term dominates the truncation error, then Eq. (7) is easily evaluated by the $\delta$ function for any prior, and Eq. (6) reduces to

[TABLE]

where $c_{k+1}=\Delta/Q^{k+1}$ as enforced by the $\delta$ function. Error bands made under this assumption are denoted by the prior with a superscript (1), e.g., ${\rm A}^{(1)}$. Further progress, with or without the first-omitted-term approximation, requires an explicit choice of priors. The relevant equations for this work, such as posteriors and DoB intervals, are contained Appendix A.

III Total NN cross section

Truncation error DoBs were estimated in Ref. Furnstahl et al. (2015b) for the $np$ total cross section at laboratory energies of 50, 96, 143, and 200 MeV from results given explicitly in Refs. Epelbaum et al. (2015b, a) using the new R = 0.9 fm EKM potential. In this work we extend these calculations to all of the new potentials and for all energies up to 350 MeV. The first step is to extract the $c_{n}$ coefficients, defined by Eq. (1), from the order-by-order calculations of the $np$ total cross section $X=\sigma(E_{\rm lab})$. Here we choose our reference scale to be the leading-order calculation $X_{\text{ref}}=\sigma_{0}$, so $c_{0}\equiv 1$ by construction. Other reasonable choices of $X_{\text{ref}}$, such as a higher-order result or the experimental value, do not substantially change the convergence pattern for this observable.

We also need to specify the high-momentum scale $\Lambda_{b}$. Here we assume that $\Lambda_{b}$ is a given quantity, and adopt the values assumed by EKM. Their choice of $\Lambda_{b}\approx 400\mbox{--}600$ MeV (the particular value depending on a regulator parameter $R$) was based on a rough analysis of residual error plots (“Lepage plots”), validated by the observation that their choices resulted in natural coefficients in the EFT series for $np$ scattering cross sections Epelbaum et al. (2015b, a). In Sec. V, we make a statistical analysis of whether the EKM choices of $\Lambda_{b}$ (or nearby values) lead to self-consistent convergence patterns for observables, and explore directly determining a posterior probability distribution for $\Lambda_{b}$.

Because we have a double expansion in $p/\Lambda_{b}$ and $m_{\pi}/\Lambda_{b}$, we must develop a prescription to define $Q$. In Ref. Furnstahl et al. (2015b) we took $Q$ to be

[TABLE]

We expect that at low momenta the expansion will be dominated by powers of $m_{\pi}/\Lambda_{b}$ and at momenta much higher than $m_{\pi}$ it will be dominated by powers of $p/\Lambda$, so the appropriate choice of $Q$ in each region follows correctly from Eq. (9). However, it is not clear how we should parameterize the crossover region. To avoid cusps at $p=m_{\pi}$, we choose to replace Eq. (9) by a smooth interpolation function for $Q$:

[TABLE]

where $n$ is a sufficiently high integer (we take $n=8$ here). But we need to examine the behavior at low energies to assess whether the implicit equal weighting of the expansions is justified.

Coefficients $c_{2}$–$c_{5}$ for the total cross section, calculated with the $R=0.9\,$fm potential and $\Lambda_{b}=600\,$MeV, are shown as functions of energy in the left panel of Fig. 2. These include the results for four individual energies from Furnstahl et al. (2015b), but now we can see the global pattern. Except for the ${\rm N}{}^{4}{\rm LO}$ coefficient around $E_{\rm lab}\approx 50$ MeV, the coefficients at any fixed energy follow a distribution with a characteristic size of about one. If $Q=p/\Lambda_{b}$ were used for all energies instead of Eq. (10), then the coefficients would grow very large as $E_{\rm lab}$ gets small (i.e., as $p\rightarrow 0$), as shown in the right panel of Fig. 2. The onset of this behavior in $E_{\rm lab}$ increases with chiral order and, for ${\rm N}{}^{4}{\rm LO}$, is the source of the large coefficient near $E_{\rm lab}\approx 50\,$MeV. This reflects the increasing sensitivity at large order to the relative contribution of the two expansions in the crossover region. We do not yet have a model to address this behavior. If we exclude the crossover region, the underlying assumption of the priors $\operatorname{pr}(c_{n}|\bar{c})$ in Table 1 that the coefficients at a given energy are distributed with a characteristic size $\bar{c}$ is validated.

The observable coefficients for the other EKM potentials are shown in Figs. 3 and 4. For each potential we have adopted the value of $\Lambda_{b}$ advocated by EKM: $\Lambda_{b}$ equal to 600 MeV for $R=0.8\,$fm, $0.9\,$fm, and $1.0\,$fm, 500 MeV for $R=1.1\,$fm, and 400 MeV for $R=1.2\,$fm. We return in Sec. V to consider different choices of $\Lambda_{b}$.

The general assumption made in constructing a posterior for the truncation error in Ref. Furnstahl et al. (2015b), that the coefficients have a characteristic magnitude or upper bound $\bar{c}$, is based on the expectation that a well-formulated EFT will have a certain uniformity in the convergence pattern of observables. That is, with each successive order there is a steady convergence implied by the value of the expansion parameter (as shown below, this corresponds with a steady improvement of the prediction for the cross section). For an integrated observable such as the total cross section, we expect this to be particularly manifested. This justifies the use of lower-order results to inform our expectations for higher-order contributions.

The pattern of coefficients for $R=0.9\,$fm shows this uniformity, which is mostly still present for $R=0.8\,$fm and $R=1.0\,$fm. In particular, we see evidence for a characteristic size for the $c_{n}$s of order unity (in practice about three). However, the uniformity deteriorates significantly as one progresses to $R=1.1$, and $1.2\,$fm. This is a consequence of the growing cutoff artifacts at larger values of $R$. As the artifacts become more prevalent, there is a decreased contribution from mid-range pion physics at ${\rm N}{}^{2}{\rm LO}$ and ${\rm N}{}^{4}{\rm LO}$, which is counteracted by an increase in the contact terms at ${\rm NLO}$ and particularly ${\rm N}{}^{3}{\rm LO}$. This reflects a partial integrating-out of pion physics, which takes us closer to a pionless EFT convergence pattern with the dominant contributions at even orders in the expansion.

In Ref. Furnstahl et al. (2015b), we analyzed results only for $R=0.9\,$fm and $R=1.2\,$fm (as reported by EKM) and only at four energies. From this limited sample we concluded that the distribution of $c_{n}$s at these energies was consistent with a common $\bar{c}$ for $R=0.9\,$fm, but this was not the case for $R=1.2\,$fm. In particular, the latter case had uniformly small coefficients for N2LO and N4LO, consistent with there being no new short-range contributions at those orders and the regulator greatly reducing the pion tensor-range contribution. Now looking globally at the $R=1.2\,$fm coefficients, we see that N2LO and N4LO stay small for the full range of $E_{\rm lab}$. If we focus on N2LO in each graph, we see that the trend of the coefficients with energy is quite similar as $R$ increases (softening the interaction), but the overall scale decreases monotonically. The situation with N4LO is similar.

Turning to the N3LO coefficients for successive values of $R$, we find a transition from negative and fairly large (order $-3$) at $R=0.8\,$fm to positive and fairly large (order $+3$) at $R=1.0\,$fm and above. $R=0.9\,$fm is the middle of this transition. This is not unnatural, but may reflect a tendency toward overfitting at N3LO (see Ref. Wesolowski and Furnstahl (2017)). Taking all the orders together, the coefficients imply that the convergence pattern for $R=1.1\,$fm and $1.2\,$fm, for which regulator artifacts are significant, is not consistent with our statistical model.

Next we estimate DoB intervals for EFT truncation errors using the extracted coefficients. We apply at each energy the formulas from Sec. II to the coefficients from that energy only. The order-by-order results for the total cross section with $R=0.9\,$fm and prior set ${\rm C}_{\epsilon}$ (using Eq. (43) or (44)) are shown in Fig. 5. To amplify the patterns, in Fig. 6 and below we plot the residuals with respect to the Nijmegen partial-wave analysis (NPWA) Stoks et al. (1993a), where the residual for a calculated observable $X$ is defined as

[TABLE]

All plots of observables and residuals are shown with solid lines for the calculation at each order, with dark and light shaded bands denoting the 68% and 95% DoB interval for the truncation error at each kinematic point. Note that the errors are not Gaussian, as the 95% bands are not twice the size of the corresponding 68% bands.

The order-by-order convergence of the calculations in Fig. 6 to the NPWA result is clear, but not surprising—the potential was fit to reproduce the NPWA in each partial wave. The pattern of DoB intervals shown in both Figs. 5 and 6 is mostly systematic: the widths tend to increase with $Q$, decrease with order, and overlap with preceding order DoBs. Although we used set C for these truncation error estimates on the cross section, the results using set A are similar. See the Supplemental Material sup for plots displaying various error band prescriptions. We return to quantify the effects of prior choice on the success rate of the error bands in Sec. V.

An exception to the systematic and intuitive DoB intervals is for the NLO calculation near 200 MeV, where the intervals vanish. This is because the prior set ${\rm C}_{\epsilon}$ makes no assumption on the minimum (or maximum) size of $\bar{c}$, so the only information for the DoB at NLO is the NLO coefficient, which vanishes in that energy range. In Ref. Furnstahl et al. (2015b) we included $c_{0}=1$ at NLO, which effectively set a lower limit of $\bar{c}=1$. Considering the NLO coefficients over the full energy range, as well as the other coefficients, it is clear that we should use a prior with a non-zero $\bar{c}_{<}$. A more informative, but not too restrictive, choice of $\bar{c}_{<}=0.25$ (and $\bar{c}_{>}=10$) is used in Fig. 7. The DoB intervals at low order are now more plausible while there is no significant difference at the two highest orders.

Of course, it is not enough that the DoB bands are plausible; they should be statistically valid. If our DoB intervals are consistent, we might expect the NPWA line in Fig. 5 or the zero line in Figs. 6 or 7 to lie outside the 68% region roughly 1/3 of the time and outside the 95% region roughly 1/20 of the time. With this in mind, a rough examination of Fig. 7 shows the bands are not ideal: the 95% DoBs appear too large for ${\rm NLO}$ and ${\rm N}{}^{2}{\rm LO}$, while they are too small for ${\rm N}{}^{3}{\rm LO}$; the 68% DoBs underestimate the error on the ${\rm NLO}$ plot; and the ${\rm N}{}^{3}{\rm LO}$ DoBs do not perform well at low energies. In Sec. V, we perform a systematic analysis using Bayesian model checking, where we evaluate how well the DoBs predict the subsequent order correction. Section VI concludes with a reassuring proof of concept (Fig. 27), which shows that, on average, our DoBs accurately assess the error of the order-by-order results when compared to NPWA data.

IV Other NN scattering observables

In this section we extend our analysis to other NN scattering observables. For convenience we have collected in the Appendix the relevant notation and formulas we have used, as well as a brief comparison to other notations in the literature. We focus on the differential cross section and a set of the most commonly considered spin observables, namely the analyzing power $A_{y}$, polarization transfer coefficients $A$ and $D$, and the spin correlation parameters $A_{xx}$ and $A_{yy}$. Each observable has been generated from LO through N4LO, primarily using the $R=0.9\,$fm potential, which we have seen demonstrates the best convergence pattern. We consider the observables both at fixed energy as a function of angle and at fixed angle as a function of energy.

In Fig. 8, residuals for the differential cross section as a function of scattering angle at fixed $E_{\rm lab}=96\,$MeV for the $R=0.9\,$fm potential are shown as a characteristic example of this observable; other energies display similar characteristics. The detailed order-by-order convergence pattern does not seem to depend on angle, suggesting that it is plausible to describe the convergence statistically.

These observations are supported by the plot of coefficients as a function of angle, shown in Fig. 9, for which each order takes a turn at being the largest in magnitude. For this extraction, the leading-order result $X_{0}$ was taken for $X_{\text{ref}}$; the results are not sensitive to this choice. The scale of the dimensionless coefficients is roughly uniform at angles less than 150∘ (about four), which is several times larger than the scale for the integrated cross section at this energy. The scale is significantly larger at back angles, where the momentum transfer becomes twice the relative nucleon momentum, which may require a reexamination of the expansion in this region. However, our overall conclusion is that the naturalness assumption, in the form of a characteristic size for the coefficient variations, still holds without integrating over all angles.

For any given spin observable $X_{pqik}$, we assume a natural expansion for the full quantity

[TABLE]

which is the probability for a particle to scatter into a solid angle $\mathinner{\mathrm{d}\Omega}$, given that the beam and target particles are polarized in the $i$ and $k$ directions and the scattered and recoil particle spins are in the $p$ and $q$ directions, respectively. The dimensionful scale in Eq. (12) is set by the size of the differential cross section, so the natural choice is $X_{\text{ref}}=\mathrm{d}\sigma/\mathrm{d}\Omega$. That means that the expansion for the spin observable itself is

[TABLE]

with no additional prefactor. Below and in the Supplemental Material sup we see that this scaling is consistent with natural ranges for the $c_{n}$s.

In Fig. 10a we show the extracted coefficients as a function of scattering angle for six $np$ scattering observables at $E_{\rm lab}=250\,$MeV, calculated using the potential with $R=0.9\,$fm. The corresponding DoB bands for the residuals at this energy, following the same prescription as applied to the cross sections, are shown in Fig 10b. As already noted, the LO coefficient does not inform our model for truncation errors and so is not shown and is not used for the truncation error posteriors. These figures serve as a representative example; figures showing coefficients and DoB bands for many additional energies are given in the Supplemental Material sup .

Except for constraints on some observables at special angles (e.g., $A_{y}$ at $\theta=0^{\circ},180^{\circ}$ and $A$ at $\theta=0^{\circ}$), the coefficients truly look like independent bounded random functions of the angle, which supports our proposition that a statistical treatment of their behavior is warranted. As with the cross sections, the DoB bands decrease in size systematically and the lower-order bands overlap the higher-order bands. It also appears that, in general, it is not necessary to look at observables integrated over all angles to see a natural EFT convergence pattern. We may therefore apply our statistical model to estimate truncation uncertainties for NN angular observables.

Figure 10c shows coefficients at a fixed angle of $\theta=120^{\circ}$ as a function of energy. The corresponding DoB intervals are shown in the Supplemental Material sup . Enhanced ${\rm N}{}^{4}{\rm LO}$ coefficients in the double-expansion crossover region, noted earlier for the total cross section, are visible in most of the observables. If this region is omitted, the behavior of the coefficients with energies varies with a $\bar{c}$ scale of about 2, with no other systematic patterns apparent.

Thus the observed convergence patterns of NN scattering observables for the EKM interaction with $R=0.9\,$fm, considered as functions of energy or angle, satisfy the statistical model naturalness assumptions implied by Fig. 1. The DoB intervals above 50 MeV derived using prior set C exhibit reasonable patterns (as do those using set A; see Supplemental Material sup ), but do not by themselves validate the statistical model. For that purpose we turn to Bayesian model checking to assess the statistical consistency of all the EKM potentials as well as the sensitivity to the choice of prior sets from Table 1.

V Model Checking

The predictiveness of our statistical model for EFT truncation errors relies on how well our implementation of naturalness aligns with the true convergence pattern exhibited by the EFT. An EFT could fail to exhibit a natural convergence pattern because of regulator artifacts or a poorly chosen $\Lambda_{b}$. Our prior sets for $\operatorname{pr}(\bar{c})$ and $\operatorname{pr}(c_{n}|\bar{c})$, which encode our assumptions about the size of the higher-order coefficients, may also be called into question.

The efficacy of our approach for any given EFT or particular observables predicted by that EFT can be examined using Bayesian model checking Gelman et al. (2013). Here we make use of consistency checks to determine if the DoB intervals behave as advertised. We also investigate the possibility of determining $\Lambda_{b}$ solely from the convergence pattern and the assumption of naturalness.

V.1 Consistency checks

Once a posterior pdf for $\Delta_{k}$ is determined via Eq. (6), the probability that the truncation error is in a DoB interval follows directly from Eq. (4). If our statistical model for the error is valid, a $(100*p)\%$ DoB interval should on average contain the actual next order value of the observable $(100*p)\%$ of the time (we use the first-omitted-term approximation in this section). By applying this test for a range of $p$ values to a sufficiently large set of observables, we can test for inaccurate models or EFTs with irregular convergence patterns. Such a consistency check provides us with the statistical toolset to analyze the sensitivity to our choice of priors and the consistency of the breakdown scale $\Lambda_{b}$ taken from Ref. Epelbaum et al. (2015b).

The procedure for creating consistency plots333Calibration plots or curves are other common names for such tools. to implement model checking is as follows Furnstahl et al. (2015b):

Choose a set of independent observables for which the next-order calculation is available (not including LO).444In Ref. Furnstahl et al. (2015b) the LO to NLO success rates were included as part of the consistency checks. Because we want to test the convergence pattern only, the LO to NLO success rate is not relevant here, as in the previous sections where we omit $c_{0}$. 2. 2.

Select a grid of $(100*p)\%$ DoB intervals with $p$ ranging from 0 to 1. 3. 3.

Compute the $(100*p)\%$ DoB interval for each observable in the set, using the same priors throughout. 4. 4.

For each next-order calculation that is within the DoB interval of the previous order, count one success. 5. 5.

Take the number of successes $n$ and divide by the total number of observables $N$ to get the actual success rate. 6. 6.

Plot the success rates versus DoB interval percentage and compare to the ideal result given by a $45^{\circ}$ line.

Because we will have a finite number $N$ of observables, we expect fluctuations away from the ideal result for a true $(100*p)\%$ success rate, as given by the binomial posterior

[TABLE]

We apply Bayes theorem with a uniform prior on $p$ to convert to a posterior for $p$:

[TABLE]

and generalize Eq. (14) to continuous $n$ to calculate horizontal 68% and 95% confidence intervals for the DoB percentage, using the HPD prescription (see Sec. II). These become shaded bands in the consistency plots.555In Ref. Furnstahl et al. (2015b) a different procedure yielded bands for given $N$, $n$, and $p$ that are reflected about the 45 degree line from the ones here. Additionally, the bands in Furnstahl et al. (2015b) were calculated using equal-tailed credible intervals for DoBs rather than the HPD prescription. Both procedures approach symmetric bands for large $N$.

We can easily evaluate DoB intervals for choices of $\Lambda_{b}$ different from those identified by EKM, which we have adopted so far. We follow Refs. Bagnaschi and Jenniches (2014); Bagnaschi et al. (2015); Furnstahl et al. (2015b) in doing this by introducing a scaling factor $\lambda$ to generalize Eq. (1) as

[TABLE]

Varying $\lambda$ about unity shifts $\Lambda_{b}$; in the consistency plots here we consider 20% variations, namely $\lambda=0.8$ and $1.2$, with respect to the EKM choice with $\lambda=1.0$.

The logic of the remainder of this section is as follows. We begin with a reexamination of the consistency plots for the total cross section, as begun in Figs. 10 and 11 from Ref. Furnstahl et al. (2015b), exploring more energies and stability under prior choice (Figs. 11–13). Next we extend the previous analysis with results from the differential cross section and our selected spin observables (Figs. 14 and 15). Finally, consistency plots of EKM potentials with different regulators are examined (Figs. 16–18), including examples of potentials that fail our analysis (Figs. 19 and 20). For a more extensive survey of our results, see the Supplemental Material sup .

In Fig. 11 we show consistency plots for the total cross section calculated with the $R=0.9\,$fm EKM potential for prior set ${\rm C}_{0.25\text{-}10}^{(1)}$ (recall that the superscript indicates that the truncated error is assumed to be given by the first omitted term). Here, each line averages over the success rate of the ${\rm NLO}$, ${\rm N}{}^{2}{\rm LO}$, and ${\rm N}{}^{3}{\rm LO}$ error bands in predicting the corresponding next-order contributions at energies $E_{\rm lab}=20,40,\ldots,340$ MeV. The trends show that for $\lambda=1$, i.e. $\Lambda_{b}=600$ MeV, the predicted $(100*p)\%$ DoB aligns with the measured success rate to within the uncertainty predicted by Eq. (15). While using data at many $E_{\rm lab}$ values improves the statistics, the independence of the results may be questionable if calculated for too closely spaced kinematic variables. Dependent measurements would cause the gray error bands in Fig. 11 to be too restrictive, so the $\lambda=0.8$ and $\lambda=1.2$ lines may be consistent even though they are generally outside the 68% bands (cf. the leftmost plot in Fig. 16). In future work we will model the correlation length in energy using Gaussian processes Rasmussen and Williams (2006); Gelman et al. (2003); Jones et al. (1998) (GPs) to draw more robust conclusions about independence.

In Fig. 12 we decompose the $\lambda=1$ line of Fig. 11 into the contribution from each individual order, while Fig. 13 shows the same decomposition but using prior set ${\rm A}_{\epsilon}^{(1)}$. Given the slight changes between Figs. 12 and 13, as well as similar examples not shown, we conclude that prior choice has little effect on the predictions of EFTs with good convergence patterns. For such an EFT, we expect the predictions to improve with the order of the prediction, because the higher orders contain more information about the pattern of the observable coefficients. This is what we see, with the ${\rm N}{}^{3}{\rm LO}$ predictions being fully consistent within the gray bands.

Next we consider angle-dependent observables, which were not analyzed in Ref. Furnstahl et al. (2015b). Each observable is generated using the $R=0.9$ fm EKM potential with $\Lambda_{b}=600$ MeV ($\lambda=1$). Each is evaluated at $N=102$ kinematic points: 17 energies ($20,40,\ldots,340$ MeV, as for $\sigma$) with 6 angles ($40^{\circ},60^{\circ},\ldots,140^{\circ}$) for each energy. The prior set used in the analysis is ${\rm C}_{0.25\text{-}10}^{(1)}$ and the plots are decomposed order-by-order. Figure 14 shows a consistency plot for the differential cross section, while Fig. 15 shows, as an example, the individual spin observables $A_{y}$, $A_{xx}$, and $A_{yy}$. The ${\rm N}{}^{2}{\rm LO}$ DoBs consistently underestimate the size of the ${\rm N}{}^{3}{\rm LO}$ correction for the differential cross section, but the ${\rm NLO}$ and ${\rm N}{}^{3}{\rm LO}$ DoBs do fairly well. Some of the DoBs for spin observables, such as $A_{y}$ and $A_{xx}$, overestimate the higher-order corrections, while the $A_{yy}$ DoB performs well. The ${\rm NLO}$ and ${\rm N}{}^{2}{\rm LO}$ coefficients of $A_{y}$ and $A_{xx}$ are generally larger than ${\rm N}{}^{3}{\rm LO}$ and particularly the ${\rm N}{}^{4}{\rm LO}$ coefficients at $E_{\rm lab}\gtrsim 100$ MeV, while $A_{yy}$ tends to have coefficients that each take turns being the largest.

The spacing in angle and energy may be close enough that the calculations used for the consistency plots are significantly correlated, which will constrain the gray error bands unnecessarily due to the large number of non-independent points. The true impact of this correlation has not yet been quantified and is a topic for future investigation. From the aforementioned plots, we can conclude that although integrating over angles is not necessarily required to ensure a natural convergence pattern of coefficients, some observables do show notable patterns that adversely affect the predictive power of their respective DoBs.

Finally, we return to the topic of EKM potentials with varying regulators, first raised in Sec. III. Thus far we have mainly focused on the $R=0.9$ fm EKM potential due to its natural convergence pattern compared to the other potentials, as evidenced by Figs. 2–4. Now we relax this focus to gain insight into the effects that regulator choices and their consequent convergence patterns have on the reliability of the error bands generated by this analysis. We also test the proposed breakdown scale $\Lambda_{b}$ for each regulator by varying $\lambda$ defined in Eq. (16) about unity. In an attempt to ensure independent results for the chosen kinematic points, in Figs. 16–18 we use $E_{\rm lab}=96,143,200,300$ MeV and $\theta=60^{\circ},120^{\circ}$ (if applicable). The choices of separation length in $E_{\rm lab}$ and $\theta$ are based on a rough analysis of the coefficient curves, which suggests that energies spaced by about $70$–$80$ MeV and $\theta$ spaced by $30^{\circ}$–$40^{\circ}$ can be taken as independent for evaluating DoB successes.

We find in Figs. 16–18, which show results averaged over orders for $R=0.9\,$fm, $0.8\,$fm, and $1.0\,$fm, respectively, that our statistical model for truncation errors is generally successful for these parameters. For $R=0.9\,$fm, both $\sigma$ and $\mathrm{d}\sigma/\mathrm{d}\Omega$ show strong consistency with $\lambda=1$, meaning $\Lambda_{b}\approx 600$ MeV, but a wider range of $\Lambda_{b}$ is not ruled out. In contrast, the spin observables are more consistent with somewhat larger $\Lambda_{b}$, particularly if we accept the limits of the gray bands. The three sets of observables for $R=0.8\,$fm remain fairly consistent with a single choice for $\Lambda_{b}$ and overall this potential passes the test of a natural convergence pattern based on the expected level of consistency. Although the order-averaged consistency plots for the $R=1.0$ fm potential are reasonable, the order-by-order convergence pattern and plausibility of a single $\Lambda_{b}$ become suspect; see the Supplemental Material sup for more information.

The failure of our statistical model for truncation errors when applied to $R=1.1$ fm and $R=1.2$ fm, which was anticipated by the pattern of coefficients in Fig. 4, is best observed in the order-by-order consistency plots, where the impact of fluctuations in coefficient size becomes clear. For better statistics (larger $N$), we use $\theta=60^{\circ},120^{\circ}$ and $E_{\rm lab}=20,40,\ldots,340$ MeV; since the chosen angles are fairly representative, using more angles does not greatly affect the conclusions. Figures 19 and 20 explicitly show the unequal nature of the coefficient magnitude for these regulator values. Because much of the physics content at ${\rm N}{}^{2}{\rm LO}$ and ${\rm N}{}^{4}{\rm LO}$ is moved to ${\rm NLO}$ and ${\rm N}{}^{3}{\rm LO}$, the ${\rm N}{}^{2}{\rm LO}$ DoBs tend to underestimate the contribution due to ${\rm N}{}^{3}{\rm LO}$, while the ${\rm NLO}$ and ${\rm N}{}^{3}{\rm LO}$ DoBs overshoot the error estimates due to ${\rm N}{}^{2}{\rm LO}$ and ${\rm N}{}^{4}{\rm LO}$, respectively. Because the trade-off of large and small coefficients causes the error bands to be overestimated and then underestimated at alternating orders in the expansion, this effect could on average cancel out when comparing to the actual data, as that comparison highlights the size of all left-out higher-order terms.

V.2 Posterior for $\Lambda_{b}$

So far we have assumed that the EFT breakdown scale $\Lambda_{b}$ was a given quantity, and then calculated posteriors for EFT truncation errors contingent on the known coefficients $c_{n}$. We have also checked whether this posterior is statistically consistent with particular fixed choices for $\Lambda_{b}$. Here we explore whether we can extract a plausible range for $\Lambda_{b}$ by calculating a posterior pdf for $\Lambda_{b}$, contingent only on the order-by-order results. We combine results from different momenta and angles that are far enough apart that it is reasonable to assume the EFT calculations are uncorrelated, but also compare to much more closely spaced kinematics to improve the statistics. The eventual goal is to be able to use modeled correlations between observable calculations to calculate $\Lambda_{b}$ based on the calculations at many different momenta $p$ and angles.

We first rewrite Eq. (1) in terms of powers of $p$ instead of $Q$ (recall that this should not be interpreted as the explicit $p$ dependence of the observable):

[TABLE]

which defines the dimensionful coefficients $b_{n}$. The $b_{n}$ are trivially related to the $c_{n}$ from Eq. (1) by

[TABLE]

We proceed based on two independent assumptions: (1) the details of the chiral EFT description of low-energy QCD (e.g., renormalization scale and scheme) dictates a well-defined breakdown scale $\Lambda_{b}$, and (2) a well-formulated EFT implementation will lead to natural expansion coefficients for observables. Although any given $b_{n}$ can be extracted from order-by-order calculations without any reference to a breakdown scale or naturalness, our assumptions imply that the value of $b_{n}$ manifests the interplay between the underlying $\Lambda_{b}$ and the natural $c_{n}$ required by Eq. (18). This relationship is represented graphically as a Bayesian network in Fig. 21.

We note that $b_{0}$ will not give any information on the expansion parameter, because it will not modify the convergence pattern. Because it is known that $b_{1}=0$ in chiral EFT, it too is omitted from our analysis.

In general, we want to use order-by-order calculations of several observables at several kinematic points to inform our estimate of $\Lambda_{b}$. The full quantity of interest is therefore

[TABLE]

where the given information is $m$ sets of coefficients $\mathbf{b}_{k}(\bm{\alpha}_{i})\equiv(b_{2}(\bm{\alpha}_{i}),\cdots,b_{k}(\bm{\alpha}_{i}))$ labeled by $\bm{\alpha}_{i}=({E_{\rm lab}}_{\,i},\theta_{i},X_{i})$: the set of kinematic parameters and observable $X_{i}$ from which the $\mathbf{b}_{k}$ were calculated. Using Bayes theorem, we can express the posterior for $\Lambda_{b}$ as

[TABLE]

We have assumed statistical independence of coefficients at different orders, but it is possible for $b_{n}(\bm{\alpha}_{i})$ to be correlated with $b_{n}(\bm{\alpha}_{j})$. The coefficients can be correlated for multiple reasons: for a given observable, the kinematic parameters may be close to one another; two distinct observables could themselves be correlated; or a combination of both. Assuming statistical independence in $\bm{\alpha}_{i}$, we can factor the likelihood in Eq. (20) into

[TABLE]

Therefore, the posterior is given by

[TABLE]

The denominator of Eq. (22) is simply a normalization constant and the prior $\operatorname{pr}(\Lambda_{b})$ can be chosen later on, leaving only $\operatorname{pr}(\mathbf{b}_{k}(\bm{\alpha}_{i})|\Lambda_{b})$ to evaluate. For simplicity, we will refer to this as $\operatorname{pr}(\mathbf{b}_{k}|\Lambda_{b})$, noting that the likelihoods for all $\mathbf{b}_{k}(\bm{\alpha}_{i})$ simply need to be multiplied together to get the final posterior pdf in Eq. (22).

To express $\operatorname{pr}(\mathbf{b}_{k}|\Lambda_{b})$ in terms of the prior assumptions of naturalness, we first use marginalization Sivia and Skilling (2006) to introduce as auxiliary parameters the dimensionless coefficients $\mathbf{c}_{k}$ [see Eq. (3)]:

[TABLE]

Next, to express the prior pdf for the coefficients $c_{n}$, we integrate in the naturalness parameter $\bar{c}$:

[TABLE]

To simplify Eq. (24) we use independence as reflected in the causal relationship outlined in Fig. 21. Neither the $c_{n}$s nor $\bar{c}$ depend on $\Lambda_{b}$ if they are not mediated by $b_{n}$. We adopt a prior of independence between the $c_{n}$s as before. Thus

[TABLE]

The $b_{n}$s also only depend on their corresponding $c_{n}$ and $\Lambda_{b}$, and are independent of one another. This means that

[TABLE]

Therefore, Eq. (24) can be written as

[TABLE]

The pdf for $b_{n}$ contingent on $c_{n}$ and $\Lambda_{b}$ is simply

[TABLE]

which enables us to perform the $c_{n}$ integrations directly. Thus,

[TABLE]

where we have used $\prod_{n=2}^{k}\Lambda_{b}^{n}=\Lambda_{b}^{k(k+1)/2-1}$ and have set $c_{n}=b_{n}\Lambda_{b}^{n}$ from now on.

To evaluate Eq. (29), we must make choices for the priors, such as those from Table 1. Analytic expressions of Eq. (29) can be found for sets ${\rm A}_{\epsilon}$ and ${\rm C}_{\epsilon}$, which we will consider here. It is reasonable to assume no prior knowledge of the scale of $\bar{c}$, i.e. allow $\bar{c}_{<}\to 0$ and $\bar{c}_{>}\to\infty$, because the scale can vary wildly with a changing $\Lambda_{b}$. For set ${\rm C}_{\epsilon}$,

[TABLE]

The result for set ${\rm A}_{\epsilon}$ is similar to Eq. (30), with the replacement $\mathbf{c}_{k}^{2}\to\bar{c}_{(k)}^{2}$, where $\bar{c}_{(k)}=\max\{|c_{i}|:c_{i}\in\mathbf{c}_{k}\}$. These likelihoods are maximized for values of $\Lambda_{b}$ where the individual $c_{n}$s are about the same size.

The final step in specifying the posterior is to make a choice of $\operatorname{pr}(\Lambda_{b})$. Here we employ a non-informative log-uniform prior as we did for $\bar{c}$:

[TABLE]

which assumes we know only limits on the scale of $\Lambda_{b}$. Then for set ${\rm C}_{\epsilon}$,

[TABLE]

where the $\theta$ functions on $\Lambda_{b}$ are implicit. A more probable region in $\Lambda_{b}$ is singled out in Eq. (32) by the interplay of the $\operatorname{pr}(\mathbf{b}_{k}(\bm{\alpha}_{i})|\Lambda_{b})$ factors, which individually favor $\Lambda_{b}$s that make the order-by-order $c_{n}$s for each $\bm{\alpha}_{i}$ about the same size.

The lower limit $\Lambda_{<}$ in the prior of Eq. (31) requires comment: if $\Lambda_{<}$ is set less than the momentum scale $p_{i}$ corresponding to the lab energy where $\mathbf{b}_{k}(\bm{\alpha}_{i})$ is calculated, the expansion parameter $Q=p_{i}/\Lambda_{b}$ may be greater than one. If we have an EFT that converges according to our statistical model, $Q>1$ for the relevant kinematic points should be excluded by Eq. (22). If it instead favors values of $\Lambda_{b}$ for which $Q>1$, this would signal an inconsistency between the truncation error model and the EFT as implemented.

As already noted, we assume that because of their separations in energy or angle, the chosen sets of kinematic parameters can be treated as independent from one another and their probability densities multiplied. We make a similar assumption for the observables themselves, i.e., the set labeled $X_{pqik}$ includes $\mathbf{b}_{k}$ sets for each of the spin observables $A_{y}$, $A$, $D$, $A_{xx}$ and $A_{yy}$. The assumption of independence, particularly for observables at the same energy or angle, may be questioned. The exploration of methods to combine data from all kinematic parameters, such as through GPs Rasmussen and Williams (2006); Gelman et al. (2003); Jones et al. (1998), and assessments of observable independence, are currently in progress.

Given the above assumptions, we have applied Eq. (32) to various potentials, observable sets, and kinematic parameters; the resulting pdfs and DoB intervals, using the HPD prescription (see Sec. II), are presented in Figs. 22–26. In contrast to central credibility intervals, the HPD intervals ensure that massive extremes, such as the ${\rm N}{}^{3}{\rm LO}$ posterior for $\sigma$ near its lower boundary in Fig. 22, are not necessarily excluded from our DoB intervals Kruschke (2011). The posteriors of Fig. 22 mirror the conclusions drawn from Fig. 16—both $\sigma$ and $\mathrm{d}\sigma/\mathrm{d}\Omega$ predict $\Lambda_{b}\approx 600$ MeV, while the set of spin observables, taken together, prefers $\Lambda_{b}>600$ MeV. The relatively small amount of data used from $\sigma$ and $\mathrm{d}\sigma/\mathrm{d}\Omega$ do not allow for a very precise determination of $\Lambda_{b}$.

Figure 23 explores this result for spin observables by splitting out the posteriors for $A_{y}$, $A_{xx}$, and $A_{yy}$ separately. These posteriors can be qualitatively predicted from the order-by-order consistency plots for these observables given in Fig. 15. In general, the strength of the $\Lambda_{b}$ posterior at ${\rm N}{}^{4}{\rm LO}$ should be highly correlated with the pattern in the consistency plot at ${\rm N}{}^{3}{\rm LO}$ (e.g., do the points lie above or below the 45 degree line, which imply that $\lambda<1$ and $\lambda>1$ are more probable than $\lambda=1$, respectively). Similarly, the $\Lambda_{b}$ posterior at ${\rm N}{}^{3}{\rm LO}$ correlates with the pattern in the consistency plot at ${\rm N}{}^{2}{\rm LO}$. For $A_{y}$, this rule predicts that the ${\rm N}{}^{4}{\rm LO}$ posterior should have its strength concentrated well above 600 MeV, and that the ${\rm N}{}^{3}{\rm LO}$ posterior should be located to its left. For $A_{xx}$, the ${\rm N}{}^{4}{\rm LO}$ posterior should also be well above 600 MeV, but the ${\rm N}{}^{3}{\rm LO}$ posterior should be to its right. Finally, for $A_{yy}$, the consistency plots predict the ${\rm N}{}^{4}{\rm LO}$ posterior will be concentrated near 600 MeV (i.e., $\lambda=1$), with the ${\rm N}{}^{3}{\rm LO}$ posterior shifted somewhat to the right. All of these expectations are realized in Fig. 23.

One may wonder to what extent the $\Lambda_{b}$ posteriors are stable under different choices of kinematic parameter sets. Figure 24 shows the posteriors as in Fig. 22, but with a different (lower) range of energies. Note that $E_{\rm lab}=50$ MeV is near the crossover region $p\sim m_{\pi}$, where the interpretation of the expansion parameter is unclear. However, while there are systematic shifts, both sets of posteriors for $\sigma$ and $\mathrm{d}\sigma/\mathrm{d}\Omega$ are consistent with the EKM value of $\Lambda_{b}=600$ MeV, with the ensemble spin observable posteriors favoring significantly higher values. In all cases there are wide posteriors at ${\rm N}{}^{3}{\rm LO}$ and more stability at ${\rm N}{}^{4}{\rm LO}$. The only major shift in the median between these energy sets is for the ${\rm N}{}^{3}{\rm LO}$ cross section result. If we use larger sets of observables more closely spaced in both energy and angle (see Fig. 25), neglecting the possible danger from correlations, the $\Lambda_{b}$ posteriors become more narrow and more Gaussian, but are systematically in accord with Fig. 22.

Figure 26 applies Eq. (32) to the $R=1.2$ fm EKM potential, for which EKM assumed that $\Lambda_{b}=400$ MeV. Most of the posteriors imply even smaller values of $\Lambda_{b}$, except for the ${\rm N}{}^{4}{\rm LO}$ posterior for the spin observables, which is completely inconsistent. Following the earlier discussion of the lower limit on $\Lambda_{b}$ in the prior of Eq. (31), we see that the posteriors are not only indicating much lower values of $\Lambda_{b}$, but they also weight areas where the expansion parameter $Q>1$. The cross section and differential cross section posteriors exhibit this behavior most, while the spin observable posteriors are maximized in regions where $Q<1$. The weighting of $Q>1$ regions of the posterior is another indication that the EFT convergence for this regulator is not well-described by the statistical model.

VI Summary and outlook

In this work, we extend the analysis from Ref. Furnstahl et al. (2015b) that applies Bayesian statistics to the quantification of theoretical uncertainties in chiral EFT. Our approach makes testable predictions of DoB error bands based upon assumptions about the convergence pattern of EFT observables and an implementation of naturalness. In particular, we assume that the scaled observable coefficients $c_{n}$ defined in Eq. (1) are effectively random functions of natural size whose magnitude provides an estimate of the error incurred by truncating the EFT expansion.

We apply this model to a set of $np$ scattering observables predicted by the semi-local chiral EFT potentials of EKM Epelbaum et al. (2015b, a), who also proposed a non-statistical protocol for uncertainty quantification. The EKM error estimates in Epelbaum et al. (2015b) correspond most closely to the leading approximation of set ${\rm A}_{\epsilon}$ (see Table 1). In particular their error bands at NkLO are $k/(k+1)*100\%$ DoB intervals Furnstahl et al. (2015b); i.e., they do not correspond to the same DoB at each order. Additionally, if the known next-order result does not lie in that $k/(k+1)*100\%$ DoB interval, EKM extends the interval to the next-order result. Therefore, at some orders it is possible to interpret the EKM intervals according to our truncation error model using set ${\rm A}_{\epsilon}$, but not always. To calculate consistent statistical DoBs at each order, we follow the statistical model outlined in Fig. 1, which assumes a natural convergence pattern for the EFT.

We begin by arguing that appropriate physical choices of scale in Eq. (1) are $X_{\text{ref}}\approx X_{0}$ for $\sigma$ and $\mathrm{d}\sigma/\mathrm{d}\Omega$, while $X_{\text{ref}}\approx 1$ for any spin observable $X_{pqik}$. We then validate a posteriori for these choices of $X_{\text{ref}}$ the natural distribution of observable coefficients for the total cross section, differential cross section, and a selection of spin observables for the $R=0.9$ fm potential. The $R=0.8$ fm potential is also consistent with our statistical model. In contrast, the convergence patterns of the $R=1.0$ fm, 1.1 fm and 1.2 fm potentials become increasingly unsystematic, and hence they become less well described by this model. The expansion parameter for these softer potentials is dominated by regulator artifacts, which EKM account for by adopting a smaller value for $\Lambda_{b}$. This reflects the effective cutoff in momentum space instead of the intrinsic breakdown scale of chiral EFT Furnstahl et al. (2015a). But the order-by-order convergence pattern is also modified as the long-range pion contributions at odd orders in $Q$ are significantly shifted to short-range contributions at even orders. This is manifested in the coefficients extracted in Fig. 4 and the failed validation of DoBs in Figs. 19 and 20.

Our higher-order results are generally insensitive to the specific prior choice; we compare sets A and C in Figs. 12 and 13, and provide more examples in the Supplemental Material sup . The results for potentials with good convergence patterns (we focus on $R=0.9\,$fm) can be summarized in terms of our progress on the questions raised in Sec. I:

The observable coefficients of the total cross section vary smoothly with energy, typically changing sign once over the energy range from 0 to 350 MeV. There is no apparent order-by-order pattern at any given energy, which supports our model of a random distribution characterized by a size $\bar{c}$. The unnatural size of the ${\rm N}{}^{4}{\rm LO}$ $\sigma$ coefficient at $E_{\rm lab}=50$ MeV was noted in Ref. Furnstahl et al. (2015b). At low energies, higher-ordered coefficients become more sensitive to the value assigned to the expansion parameter $Q$ around the crossover region $p\sim m_{\pi}$, for which we do not have a model. Hence we cannot make strong statements about the coefficient spectrum and its implication for naturalness in that energy region. We plan to test alternative schemes for $Q$ in the crossover region and to validate the presence of the crossover in the $c_{n}$s using a change-point analysis of the correlations modeled using a GP model Gramacy and Lee (2008); Rasmussen and Williams (2006); Jones et al. (1998); Gelman et al. (2003). 2. 2.

The observable coefficients for both the differential cross section and the chosen spin observables vary smoothly in both $E_{\rm lab}$ and $\theta$ with characteristic sizes between about 1 and 5 for $R=0.9\,$fm, which validates the assumption that naturalness propagates to these observables for this potential. The functional dependences show no obvious patterns, supporting the model of effectively random functions. As with the cross section, the interpretation of the coefficients is ambiguous for low energies. 3. 3.

Because each $c_{n}$ is a smooth function when plotted against both $E_{\rm lab}$ and $\theta$, the values of the observable coefficients at one value of the kinematic parameters are correlated within some neighborhood (a correlation length) of $E_{\rm lab}$ and $\theta$. Through a rough estimation, we find that the correlation length in energy is about 80 MeV, while the correlation length in $\theta$ is approximately $40^{\circ}$. These values were estimated visually here, but in the future we will determine them directly using a GP model for the $c_{n}$s Rasmussen and Williams (2006); Jones et al. (1998); Gelman et al. (2003). This additional information will then be incorporated into our statistical model for truncation uncertainties. 4. 4.

The checks in Sec. V.1 show that taking $\Lambda_{b}$ to be the same scale for both $\sigma$ and $\mathrm{d}\sigma/\mathrm{d}\Omega$ is statistically consistent for the $R=0.9\,$fm EKM potential. While the spin observable $A_{yy}$ is also consistent with that same scale, the ensemble of spin observables are more consistent with a somewhat larger value. The $R=0.8$ fm EKM potential also shows promise, but extracting $\Lambda_{b}$ becomes questionable for $R=1.0$ fm and worse for larger $R$. The posteriors for $\Lambda_{b}$, shown in Sec. V.2 and the Supplemental Material sup , lead to the same conclusions—probable ranges of $\Lambda_{b}$ consistent with the values proposed by EKM can be extracted for $R=0.8$ and $0.9$ fm, identifying probable ranges for $R=1.0\,$fm is questionable at best, and the other EKM potentials ($R=1.1$ and $1.2\,$fm) are not well described by our statistical model. Conclusions about $\Lambda_{b}$ are not warranted for the poorly behaving potentials.

An overall validation of our truncation error model as applied to the EKM potential with $R=0.9\,$fm is provided in Fig. 27, which shows order-by-order consistency plots for observables compared with the NPWA data. These plots are made by modifying step 4. of the procedure laid out in Sec. V.1 to count a success when the actual NPWA result is within the DoB interval at that order (as opposed to comparison to the next-order calculation). In general, the $(100*p)\%$ error bands work as advertised, predicting the discrepancy with the NPWA data at the $(100*p)\%$ level to within expected fluctuations.

The general success of the model for chiral EFT truncation errors motivates additional applications, further development of the model (e.g., GP models), and its full integration into parameter estimation of LECs. We plan to apply our truncation error model to other chiral interactions that are available order-by-order, such as the recent potential of Entem, Machleidt, and Nosyk in Entem et al. (2017). Our error model and Bayesian model checking diagnostics can be applied not only for other chiral interactions but also for other EFTs in general. They also apply generically to any observable calculation that fulfills the expansion model in Eq. (1), including calculations in perturbation theory. A Bayesian-type Lepage plot analysis Lepage et al. (2002); Furnstahl et al. (2015a); Wesolowski et al. (2016) of the power-law behavior of residuals as a function of energy/momentum will complement the statistically motivated model checks of Sec. V. Now that we have a framework of testable assumptions for treating $\Lambda_{b}$ as a random variable in the posterior pdf calculations of Sec. V, that information can be used to marginalize over $\Lambda_{b}$ as an auxiliary parameter in truncation error estimates rather than using a fixed $\Lambda_{b}$ value. Work in these areas is in progress.

Acknowledgements.

We thank E. Epelbaum, H. Griesshammer, N. Klco, and D. Phillips for useful discussions. Useful feedback on the manuscript was provided by N. Klco and D. Phillips. This work was supported in part by the National Science Foundation under Grant Nos. PHY–1306250 and PHY–1614460 and the NUCLEI SciDAC Collaboration under Department of Energy Grant DE-SC0008533.

Appendix A Derivations of $\Delta_{k}$ Posteriors

Here we continue the analysis started in Sec. II by giving explicit forms of posteriors $\operatorname{pr}_{h}(\Delta|\mathbf{c}_{k})$ for various prior sets. The most non-informative case of set A follows if we take $\bar{c}_{<}=\epsilon$, $\bar{c}_{>}=1/\epsilon$ and then take the limit $\epsilon\rightarrow 0$ at the end. We designate this as set ${\rm A}_{\epsilon}$, and the results for this set were first worked out in Cacciari and Houdeau (2011). If we further adopt the first-omitted-term approximation, designated ${\rm A}_{\epsilon}^{(1)}$, we have analytic expressions for Eq. (6),

[TABLE]

and for $d_{k}^{(p)}$ from Eq. (4),

[TABLE]

where

[TABLE]

and $n_{c}$ is the number of relevant known coefficients—here, $n_{c}=k-1$ since $c_{0}$ and $c_{1}$ do not contribute to our analysis, but equations are given in a general form for the reader.

Relaxing the first-omitted-term approximation can pose a numerical challenge, since the integration volume grows quickly with increasing $h$. Luckily, by following Ref. Coello Pérez and Papenbrock (2015), whose results we reproduce in Eqs. (36)–(38), Eq. (7) can be reduced to one integral for the hard-wall (hw) prior $\operatorname{pr}(c_{n}|\bar{c})$ in sets A and B, and exactly evaluated for the Gaussian (G) prior in set C. For sets A and B,

[TABLE]

which, for $h\to\infty$, becomes the atomic function $h_{a}$ Victor et al. (2009). For set C,

[TABLE]

where

[TABLE]

Equation (37) is easily evaluated for all $h$; we use $h=10$ unless otherwise specified, at which point the posteriors have converged numerically (see Ref. Furnstahl et al. (2015b)).

With Eq. (37), we can exactly evaluate Eq. (6) for set C in terms of special functions. By inserting the priors and making the variable substitution $x=1/\bar{c}$,

[TABLE]

where, of course,

[TABLE]

Equation (39) can be evaluated in terms of the incomplete $\Gamma$ functions via

[TABLE]

using the definition

[TABLE]

Of all sets in Table 1, the posterior as given by set C, via Eq. (A), is the only one for which we have found a closed-form expression for all $h$ and ranges of $\bar{c}$.

For the non-informative set ${\rm C}_{\epsilon}$, where $\bar{c}_{<}\to 0$ and $\bar{c}_{>}\to\infty$, Eq. (A) simplifies to a $t$-distribution:

[TABLE]

Rather than integrating Eq. (43), $d_{k}^{(p)}$ can be found by numerically solving the transcendental equation

[TABLE]

where ${}_{2}F_{1}$ is the hypergeometric function.

Appendix B Details on NN observables

For convenience of the reader and because of the multitude of different conventions in the literature, we have gathered in this appendix the formulas used here in the calculation of NN observables MacGregor et al. (1960); Bystricky et al. (1978); La France and Winternitz (1980); Moravcsik et al. (1989); Stoks and de Swart (1990); Stoks et al. (1993b); Carlsson et al. (2016).

B.1 Kinematics

In the context of NN scattering, one particle (the beam), with kinetic energy $E_{\rm lab}$, is incident on a stationary particle (the target). For $np$ scattering, the lab system (l.s.) is often taken to be the rest frame of the initial proton. In the center-of-momentum (c.m.) system, each particle has a relative momentum of $p_{\rm rel}$. It is convenient to relate these quantities for each NN experiment:

[TABLE]

where relativistic kinematics is used Epelbaum et al. (2015b). Unless otherwise stated, $\theta$ is the c.m. polar scattering angle while $\phi$ denotes the azimuthal scattering angle. For our purposes, $\phi$ can be set to zero because all observables can be defined relative to the scattering plane.

The spin states of the initial and final states can be expressed in the uncoupled basis $|i\rangle_{\text{spin}}=|m_{1}^{\prime}m_{2}^{\prime}\rangle$ and $|f\rangle_{\text{spin}}=|m_{1}m_{2}\rangle$, respectively, where we have suppressed $s_{1}=s_{2}=1/2$. We can also use the coupled singlet-triplet basis, where $|i\rangle_{\text{spin}}=|s^{\prime}m^{\prime}\rangle$ and $|f\rangle_{\text{spin}}=|sm\rangle$.

B.2 Observables

Because nucleons have nonzero intrinsic spin, observables in general are dependent not only on kinematic variables ($E_{\rm lab},\theta,\phi$), but also on the relative orientation of the particles’ spin. A generic spin observable can be written as

[TABLE]

where $\mathrm{d}\sigma/\mathrm{d}\Omega$ is the (unpolarized) differential cross section, $M(k_{f},k_{i})$ is the spin-scattering matrix and $\sigma_{v}=\bm{\sigma}\cdot\bm{v}$. The subscripts $p$, $q$, $i$, and $k$ refer to the polarization directions of the scattered, recoil, beam, and target particles, respectively. If a final-state subscript is zero, its polarization is not analyzed. If a initial-state subscript is zero, the corresponding particle was unpolarized.

When an observable is considered in the c.m. system, the polarization of each particle is often decomposed in a common basis using the unit vectors $\bm{\ell}$, $\bm{m}$, $\bm{n}$ defined as

[TABLE]

and shown in Fig. 28. Here we consider pure experiments, where the spin projections are solely along the basis vectors. Hence, for a c.m. observable, the subscripts $p$, $q$, $i$, and $k$ are some combination of $\ell$, $m$, $n$, and 0.

It is often convenient to express spin observables in the l.s., where the scattered and recoil particles deflect at angles $\theta_{1}$ and $\theta_{2}$, respectively. Lab system observables often use three sets of bases to define spin observables, defined by the beam, scattered, and recoil particle directions. The beam (scattered, recoil) frame aligns $\bm{k}$ ($\bm{k}^{\prime}$, $\bm{k}^{\prime\prime}$) with the lab particle momentum and defines $\bm{n}$ ($=\bm{n}^{\prime}=\bm{n}^{\prime\prime}$) to be normal to the scattering plane, which leaves $\bm{s}$ ($\bm{s}^{\prime}$, $\bm{s}^{\prime\prime}$) in the scattering plane such that $\bm{s}=\bm{n}\times\bm{k}$ ($\bm{s}^{\prime}=\bm{n}\times\bm{k}^{\prime}$, $\bm{s}^{\prime\prime}=\bm{n}\times\bm{k}^{\prime\prime}$). The initial-state subscripts $i$ and $k$ are then chosen to be $k$, $s$, $n$, or [math]. Similarly, the scattered-state subscript $p$ is $k^{\prime}$, $s^{\prime}$, $n$ or [math], and the recoil-state subscript $q$ is $k^{\prime\prime}$, $s^{\prime\prime}$, $n$ or [math].

One added complication of calculating l.s. observables involves accounting for the relativistic spin rotation angles

[TABLE]

which rotate the primed and double primed vectors about $\bm{n}$, respectively. It is the rotated vectors, denoted with subscripts $R_{1}$ and $R_{2}$, that correspond to the l.s. subscripts for the scattered and recoil particles. In the nonrelativistic case, $\alpha=0$ and $\beta=\pi/2$, which implies that

[TABLE]

This too is illustrated in Fig. 28. Whether an observable is defined in the c.m. system or the l.s. should be clear from the chosen subscripts. All of the observables considered here use the l.s. notation.

Notational inconsistencies abound in the literature. While the subscripts of Eq. (48) completely determine a given spin observable, often the $X$ is changed to match historical usage. Other times, the subscript notation is abandoned completely for a nondescript letter. Table 2 attempts to reconcile some differences by matching a common name with Eq. (48) and other popular notations found in literature.

The spin-scattering matrix is the part of the scattering $S$ matrix that is due to interactions. In our convention, they are related via $M=(2\pi/ip)(S-1)$. To evaluate $M$, it is useful to write it in singlet-triplet space, and then make a partial-wave expansion

[TABLE]

where $\hat{J}\equiv\sqrt{2j+1}$ and $\hat{L}\equiv\sqrt{2\ell+1}$, the second line contains two Wigner $3j$ symbols, and $\mathbf{J}=\mathbf{L}+\mathbf{S}$ is the total angular momentum decomposed into orbital and intrinsic angular momentum. The nuclear potential conserves the total angular momentum $j$, but generally mixes the states of $\ell$ and $s$. Equation (B.2) becomes useful only if a small number of $j$ waves are needed to accurately determine $M_{m^{\prime}m}^{s^{\prime}s}$.

When a partial wave state is uncoupled, such as when $j=0$, then $S^{j}$ can be parameterized by a real phase $\bar{\delta}^{sj}_{\ell}$ such that $S^{j}\equiv e^{2i\bar{\delta}^{sj}_{\ell}}$. For $j>0$, $S^{j}$ is 4-dimensional in angular momentum space and can be written compactly using the triplet submatrix $S_{T}^{j}$ and the singlet-triplet submatrix $S_{ST}^{j}$ via

[TABLE]

The triplet submatrix $S_{T}^{j}$ can be parameterized by introducing another real parameter, the mixing angle $\bar{\epsilon}_{j}$. Using the common notation that subscripts $+$ and $-$ refer to $\ell=j+1$ and $j-1$, respectively, then

[TABLE]

and similarly

[TABLE]

In the present work, $\bar{\gamma}_{j}=0$ for all $j$, leaving the singlet-triplet submatrix uncoupled and thus $s^{\prime}=s$.

Equations (58) and (59) employ the “Stapp”- or “bar”-phase shift parameterization. Another parameterization, with phases and mixing angle denoted here by $\delta^{sj}_{\ell}$ and $\epsilon_{j}$, was made by Blatt and Biedenharn Blatt and Biedenharn (1952):

[TABLE]

where

[TABLE]

The Blatt eigenphases are related to the Stapp phases via

[TABLE]

Given the partial-wave-projected potential $V_{\ell^{\prime}\ell}^{sj}(p^{\prime},p)$, it is convenient and numerically accurate to solve the Lippmann-Schwinger (LS) equation with standing wave boundary conditions,

[TABLE]

for the partial-wave-projected $R$ matrix (known as the $K$-matrix in other contexts). In the present work we use Gaussian quadrature to reduce the LS equation to a system of linear equations, from which $R_{\ell^{\prime}\ell}^{sj}$ is extracted Landau (1996). Having solved Eq. (B.2), the on-shell matrix $R^{sj}_{\ell^{\prime}\ell}=R_{\ell^{\prime}\ell}^{sj}(p_{\rm rel},p_{\rm rel};p_{\rm rel}^{2})$ then leads directly to the phases and mixing angle. For the uncoupled channels,

[TABLE]

For coupled channels, $R^{1j}_{\ell^{\prime}\ell}$ is 2-dimensional. The Blatt-Biedenharn phases $\delta^{sj}_{\ell}$ and $\epsilon_{j}$ are extractable via

[TABLE]

which can then be converted to the Stapp convention using Eqs. (62)–(64).