# Perturbative computations of neutron-proton scattering observables using
renormalization-group invariant $\chi$EFT up to N3LO

Oliver Thim<EMAIL_ADDRESS>Andreas Ekström  Christian Forssén
Department of Physics, Chalmers University of Technology, SE-412 96, Göteborg,
Sweden

###### Abstract

We predict neutron-proton scattering cross-sections and polarization
observables up to next-to-next-to-next-to leading order in a renormalization-
group invariant description of the strong nucleon-nucleon interaction. Low-
energy constants are calibrated to phase shifts, sub-leading corrections are
computed in distorted-wave perturbation theory, and we employ momentum-cutoff
values 500 and 2500 MeV. We find a steady order-by-order convergence and
realistic descriptions of scattering observables up to a laboratory scattering
energy of approximately 100 MeV. We also compare perturbative and non-
perturbative calculations for phase shifts and cross sections and quantify how
unitarity is gradually restored at higher orders. The perturbative approach
offers an important diagnostic tool for any power counting and our results
suggest that the breakdown scale in chiral effective field theory might be
significantly lower than estimates obtained in non-perturbative calculations.

##  I Introduction

Nuclear potentials used in ab initio [1] computations of atomic nuclei [2] are
almost exclusively derived using chiral effective field theory ($\chi$EFT) [3,
4, 5] based on Weinberg power counting (WPC) [6, 7]. Such potentials [8, 9,
10, 11, 12, 13, 14], now derived up to the fifth chiral order [15, 16, 17],
have furnished a wide range of structure and reaction predictions across the
nuclear chart [18, 19], but at the same time they grapple with the
renormalization challenge inherent to chiral nuclear forces [20]. Indeed,
numerical studies [21] of the nucleon-nucleon scattering amplitude have shown
that the contact operators, accounting for unresolved short-range physics,
already at leading order (LO) in WPC are not sufficient to renormalize the
singular nature [22] of the one pion-exchange potential. Consequently, LO
predictions based on WPC exhibit an unphysical dependence on the cutoff
$\Lambda$ that regularizes the amount of high-momentum (or short-range)
physics that is resolved.

Several PCs leading to renormalization-group (RG) invariant nucleon-nucleon
amplitudes have been proposed in the past two decades [23, 24, 25, 26, 27, 28,
29, 30, 31, 32, 33, 34, 35]. They can collectively be referred to as modified
Weinberg power countings (MWPCs). However, we typically know very little about
their predictive power for nuclei beyond the lightest-mass systems [36]. The
one exception is the recent study by Yang  _et al._ [37] that presented the
first ab initio predictions of binding energies in 4He, 6Li, and 16O using
$\chi$EFT potentials up to next-to-leading order (NLO) in several different
MWPCs. The calculations in that work revealed an $\alpha$-decay instability in
the ground states in 6Li and 16O. Subsequent analyses brought forward probable
causes for this instability as originating in ($i$) overfitting of the low-
energy constants (LECs) that parameterize the short-range interactions [38]
and ($ii$) underestimating the importance of few-nucleon forces [39] at LO in
MWPC.

The notable absence of MWPC-based predictions for heavier-mass nuclei is
likely due to a variety of factors. Firstly, potentials based on WPC are
easier to implement in present ab initio computer codes as one
straightforwardly sum leading and sub-leading corrections to the potential
before solving the Schrödinger equation, whereas in MWPC sub-leading
corrections should be added in perturbation theory [40]. Secondly, there
exists several widely available computer codes for evaluating matrix elements
of chiral nucleon-nucleon and three-nucleon potentials, as well as currents,
to very high orders in WPC. Finally, it is currently prohibitively costly to
converge ab initio predictions of nuclear properties at the large values of
the cutoff required for analyzing RG-invariance in MWPC.

In light of these facts we certainly see the utility of WPC, which might
provide a consistent EFT framework provided that renormalization is
interpreted in a fashion where the cutoff never exceeds the order of the
breakdown scale [41, 42, 43, 44]. However, the existence of MWPCs, where
renormalization does allow for the cutoff to be taken far beyond the breakdown
scale, calls for a continued effort. Given the fundamental importance of RG-
invariance it should be seriously explored whether MWPC approaches can furnish
a realistic and predictive framework for ab initio nuclear physics.

In this paper, we contribute to the meager list of quantitative predictions
grounded in RG-invariant formulations of $\chi$EFT. To the best of our
knowledge, and somewhat surprisingly, nucleon-nucleon scattering observables
have not been computed in MWPC beyond LO [41]. Here, we present predictions
for integrated and differential cross-sections, as well as polarization
observables, for elastic neutron-proton ($np$) scattering up to next-to-next-
to-next-to-leading order (N3LO) in the MWPC of Long and Yang [30, 45, 32],
where higher-order corrections to the potential are treated perturbatively
[21, 40]. This work serves as an important step in the development and
uncertainty quantification of any model of the nuclear interaction [46, 47,
48, 49, 50].

In Section II we review how to construct potentials in the PC of Long and
Yang, describe how to numerically compute the scattering amplitude in
distorted-wave perturbation theory, and explain how we calibrated LEC values.
In Section III we present results for scattering observables up to N3LO, and
we summarize and conclude in Section IV.

##  II Formalism

In $\chi$EFT, scattering amplitudes are expanded in a dimensionless ratio
$(Q/\Lambda_{b})^{\nu}$. Here, $\nu$ indicates the chiral order, $\Lambda_{b}$
is the underlying high-momentum scale of $\chi$EFT, and $Q$ denotes the
relevant low-energy scale. For nucleon-nucleon scattering, we assume
$Q\approx\text{max}(p,m_{\pi})$, where $p$ is the relative momentum in the
center of mass (c.m.) frame of the interacting nucleons, and the pion mass
$m_{\pi}$ is the relevant low-energy mass scale. In this work we adopt a
nomenclature where LO scales as $\left(Q/\Lambda_{b}\right)^{0}$ while sub-
leading orders are denoted by their relative scaling to LO. As such, NLO
scales as $\left(Q/\Lambda_{b}\right)^{1}$, next-to-next-to-leading order
(N2LO) as $\left(Q/\Lambda_{b}\right)^{2}$ and so on. In what follows, we
summarize relevant details regarding the MWPC that we use in this work, define
the potential terms $V^{(\nu)}$ entering at each chiral order, and explain how
we performed the perturbative calculations of scattering amplitudes.

###  II.1 The nucleon-nucleon interaction potential in the Long and Yang power
counting

We employ the MWPC of Long and Yang [30, 32, 51, 40], which adheres to the
following overarching principles:

  * •

The chiral order of a pion-exchange diagram, along with the necessary
counterterms for renormalizing pion loops, is determined by the naive
dimensional analysis (NDA) of its non-analytic part. This follows the same
principle as in Weinberg Power Counting (WPC).

  * •

Counterterms are promoted to lower chiral order only when needed to fulfill
the requirement of RG-invariance.

  * •

All corrections to the potential beyond LO are included perturbatively to
obtain RG-invariant amplitudes.

One-pion exchange (OPE) enters at LO in $\chi$EFT and must be treated non-
perturbatively, at least in the low partial waves where it is sufficiently
strong. The singular nature of OPE is increasingly alleviated by the
centrifugal barrier. Thus, at some point in the partial-wave expansion there
is sufficient angular momentum $\ell$ to furnish a perturbative treatment of
OPE [29, 52, 53] and consider it sub-leading.

At LO in the MWPC by Long and Yang, the OPE potential $V^{(0)}_{1\pi}$ is
considered non-perturbative in the ${}^{1}S_{0}$, ${}^{3}P_{0}$,
${}^{1}P_{1}$, ${}^{3}P_{1}$, ${}^{3}S_{1}\mathrm{-}^{3}D_{1}$ and
${}^{3}P_{2}\mathrm{-}^{3}F_{2}$ channels. OPE is attractive in ${}^{3}P_{0}$
and ${}^{3}P_{2}$. Renormalization requires promotion of counterterms to the
corresponding channels of the LO contact potential $V^{(0)}_{\mathrm{ct}}$
[21], thereby extending it beyond the canonical non-derivative ${}^{1}S_{0}$
and ${}^{3}S_{1}$ counterterms. At sub-leading orders ($\nu>0$), two pion-
exchange, $V_{2\pi}^{(\nu)}$, as well as higher-order contact potentials,
$V_{\text{ct}}^{(\nu)}$, enter perturbatively according to the principles
presented in the beginning of this subsection. The contributions to the
potential up to N3LO in the ${}^{1}S_{0}$, ${}^{3}P_{0}$, ${}^{1}P_{1}$,
${}^{3}P_{1}$, ${}^{3}S_{1}\mathrm{-}^{3}D_{1}$ and
${}^{3}P_{2}\mathrm{-}^{3}F_{2}$ channels are listed in the third column of
Table 1 labeled ”non-perturbative (at LO) channels”.

See Appendix A for detailed expressions of the potentials appearing in Table
1. Following Long and Yang, we do not consider any higher-order corrections to
OPE and employ potential expressions where pion loops are treated in
dimensional regularization. For the sub-leading two-pion exchange potential
$V^{(3)}_{2\pi}$ we use pion-nucleon LECs $c_{1},c_{3},c_{4}$ with central
values from the Roy-Steiner analysis in Ref. [54].

Table 1: Potential contributions at each in channels where OPE is treated non-perturbatively (column three) and perturbatively (column four). Detailed expressions for the potentials can be found in Appendix A. |  | non-perturbative (at LO) | purely perturbative
---|---|---|---
order | potential | channels | channels
LO | $V^{(0)}$ | $V^{(0)}_{1\pi}+V^{(0)}_{\mathrm{ct}}$ | 0
NLO | $V^{(1)}$ | $V^{(1)}_{\mathrm{ct}}$ | $V^{(0)}_{1\pi}$
N2LO | $V^{(2)}$ | $V^{(2)}_{2\pi}+V^{(2)}_{\mathrm{ct}}$ | 0
N3LO | $V^{(3)}$ | $V^{(3)}_{2\pi}+V^{(3)}_{\mathrm{ct}}$ | $V^{(2)}_{2\pi}$

Let us now turn to the channels with $\ell>1$ (and without any coupling to
$\ell\leq 1$). For these channels we consider OPE to be perturbative and
consequently set it to zero at LO. We follow Ref. [52] and suppress two-pion
exchanges by the same chiral power as OPE. Up to N3LO, there are no contact
potentials in the perturbative channels, and the contributions are listed in
the last column of Table 1. Other suggestions for the PC in perturbative
channels are discussed by, e.g., Pavón Valderrama  _et al._ [27].

###  II.2 A perturbative treatment of nucleon-nucleon scattering amplitudes

The perturbative computation of nucleon-nucleon scattering amplitudes proceeds
in two steps. First, we solve the Lippmann-Schwinger (LS) equation for the LO
amplitude in the ${}^{1}S_{0}$, ${}^{3}P_{0}$, ${}^{1}P_{1}$, ${}^{3}P_{1}$,
${}^{3}S_{1}\mathrm{-}^{3}D_{1}$ and ${}^{3}P_{2}\mathrm{-}^{3}F_{2}$
channels. Note that the LO potential is identically zero in all other
channels. Second, we perturbatively include higher-order potential corrections
to the amplitude, accounting for the distortion due to the non-perturbative LO
solution where necessary. In the following, we explain this procedure in
detail, see also Refs. [53, 32, 30].

The neutron-proton Hamiltonian in the center-of-mass (c.m.) frame can be
written

$H=\frac{\bm{p}^{2}}{m_{N}}+V_{\mathrm{I}}+V_{\mathrm{II}},$ (1)

where $\bm{p}$ denotes the c.m. momentum and $m_{N}=2m_{n}m_{p}/(m_{n}+m_{p})$
the nucleon mass. The projectile energy in the laboratory frame will be
denoted $T_{\mathrm{lab}}$. Furthermore, $V_{\mathrm{I}}$ denotes the LO
potential, and $V_{\mathrm{II}}$ denotes the sum of all sub-leading
potentials, which formally can be infinitely many. The PC helps us identify
important and less important contributions to the scattering amplitude $T$ and
therefore facilitates a meaningful truncation of $V_{\mathrm{II}}$. With the
notation for the chiral potentials $V^{(\nu)}$ introduced in Section II.1,
$V_{\mathrm{I}}$ and $V_{\mathrm{II}}$ read

$\displaystyle V_{\mathrm{I}}$ $\displaystyle=V^{(0)},$ (2) $\displaystyle
V_{\mathrm{II}}$ $\displaystyle=\sum_{\nu=1}^{\infty}V^{(\nu)}.$ (3)

The LO amplitude, $T^{(0)}$, is obtained (non-perturbatively) by solving the
LS-equation

$T^{(0)}=V^{(0)}+V^{(0)}G^{+}_{0}T^{(0)},$ (4)

where the free resolvent is given by

$G^{+}_{0}=\left(E-H_{0}+i\epsilon\right)^{-1},$ (5)

and $H_{0}=\bm{p}^{2}/m_{N}$. We use a notation where we suppress the explicit
dependence on the c.m. scattering energy, $E$, for the resolvents and
amplitudes.

In WPC, higher-order corrections are accounted for non-perturbatively by
solving the LS-equation for the sum $V_{\mathrm{I}}+V_{\mathrm{II}}$. In MWPC,
however, potentials beyond LO, i.e., the corrections ($V_{\mathrm{II}}$),
enter in perturbation theory to obtain RG invariant results [40]. Indeed,
higher-order corrections should be amenable to a perturbative treatment. If
not, they are non-perturbative in nature and belongs at LO.

Distorted-wave perturbation theory has been applied to compute scattering
amplitudes in several previous studies, see, e.g., Refs. [28, 53, 30, 32, 51,
55]. The perturbation series for the scattering amplitude can be derived and
expressed in various ways. The one that we find most instructive follows Refs.
[56, 57]. First, using the two-potential trick, the $T$-operator for the
Hamiltonian in Eq. 1 is written in the form

$T=T^{(0)}+\Omega^{\dagger}_{-}V_{\mathrm{II}}\sum_{n=0}^{\infty}\left(G^{+}_{1}V_{\mathrm{II}}\right)^{n}\Omega_{+},$
(6)

where the Møller wave operators are defined as

$\displaystyle\Omega_{+}$ $\displaystyle=\mathds{1}+G^{+}_{0}T^{(0)},$ (7)
$\displaystyle\Omega^{\dagger}_{-}$
$\displaystyle=\mathds{1}+T^{(0)}G^{+}_{0},$ (8)

and the full LO resolvent reads

$G^{+}_{1}=\Omega_{+}G^{+}_{0}.$ (9)

Inserting Eq. 3 in Eq. 6 gives for the full $T$-operator

$T=T^{(0)}+\Omega^{\dagger}_{-}\left[\sum_{\nu=1}^{\infty}V^{(\nu)}\right]\sum_{n=0}^{\infty}\left[G^{+}_{1}\left(\sum_{\nu^{\prime}=1}^{\infty}V^{(\nu^{\prime})}\right)\right]^{n}\Omega_{+}.$
(10)

Expanding both sums and organizing terms according to their chiral orders
$\nu$ yields the expressions for the first-, second-, and third-order
corrections to the LO amplitude as

$\displaystyle T^{(1)}$ $\displaystyle=\Omega^{\dagger}_{-}V^{(1)}\Omega_{+}$
(11) $\displaystyle T^{(2)}$
$\displaystyle=\Omega^{\dagger}_{-}\left(V^{(2)}+V^{(1)}G^{+}_{1}V^{(1)}\right)\Omega_{+}$
(12) $\displaystyle T^{(3)}$
$\displaystyle=\Omega^{\dagger}_{-}\Big{(}V^{(3)}+V^{(2)}G^{+}_{1}V^{(1)}+V^{(1)}G^{+}_{1}V^{(2)}+$
$\displaystyle+V^{(1)}G^{+}_{1}V^{(1)}G^{+}_{1}V^{(1)}\Big{)}\Omega_{+}.$ (13)

A diagrammatic representation of amplitudes up to NLO is presented in Fig. 1.
Note that the full amplitude at, e.g., third order (N3LO) is given by the sum
$T^{(0)}+T^{(1)}+T^{(2)}+T^{(3)}$. Clearly, the distorted-wave corrections in
Eqs. 11, 12 and 13 simplify dramatically when applied to the channels where
OPE is perturbative such that $T^{(0)}=0$, $\Omega_{+}=\mathds{1}$, and
$\Omega^{\dagger}_{-}=\mathds{1}$. In these channels we therefore recover
ordinary perturbation theory.

Figure 1: Diagrammatic representation of the LO neutron-proton amplitude
$T^{(0)}$ (hatched oval), obtained by solving the LS-equation, as well as the
first correction $T^{(1)}$ given in Eq. 11. The grey (black) solid blobs
represent the potentials $V^{(0)}$ ($V^{(1)}$).

The distorted-wave corrections to the amplitudes $T^{(\nu>0)}$ can
alternatively be obtained as solutions to a set of modified LS-type equations,
discussed in more detail in Refs. [58, 59], which read

$T^{(\nu)}=V^{(\nu)}+\sum_{i=1}^{\nu}V^{(i)}G^{+}_{0}T^{(\nu-i)}+V^{(0)}G^{+}_{0}T^{(\nu)}.$
(14)

We use this formulation to verify our numerical implementation of Eqs. 11, 12
and 13. We note that the alternative approach of modified LS-equations
requires a matrix inversion at each order, whereas the distorted-wave approach
requires matrix multiplications only. However, the number of matrix
multiplications increases rapidly as the chiral order is increased. For
example, at $\nu=10$, Eqs. 11, 12 and 13 require an order of magnitude more
matrix multiplications than the modified LS equations in Eq. 14. In this study
we only go to $\nu=3$ for which the number of matrix multiplications of the
two formulations are similar.

###  II.3 Numerical implementation

We project potentials and amplitudes to a partial-wave basis of states
$\ket{p,\ell,s,j}$ following the prescription in Ref. [60]111Note the mistake
in Eq. (4.22) pointed out in Ref. [4].. Here, $p=|\bm{p}|$, while $s,\ell,j$
denote the quantum numbers of the two-nucleon spin, orbital angular momentum,
and total angular momentum, respectively. Partial-wave matrix elements are
denoted by

$V^{js}_{\ell^{\prime}\ell}(p^{\prime},p)=\braket{p^{\prime},\ell^{\prime},s,j}{V}{p,\ell,s,j},$
(15)

where the conserved quantum numbers $s$ and $j$ are given as superscripts.

In the LS-equation, as well as in Eqs. 11, 12 and 13, infinite momentum
integrals appear and all potentials are regulated according to

$V^{js}_{\ell^{\prime}\ell}(p^{\prime},p)\to f_{\Lambda}(p^{\prime})\
V^{js}_{\ell^{\prime}\ell}(p^{\prime},p)f_{\Lambda}(p),$ (16)

where we choose a regulator function

$f_{\Lambda}(p)=\exp\left[-\frac{p^{6}}{\Lambda^{6}}\right]$ (17)

at all orders up to N3LO. In the calibration of the LECs, we use the cutoff
values $\Lambda=500$ MeV and $\Lambda=2500$ MeV.

Using Eqs. 7, 8 and 9, the terms in Eqs. 11, 12 and 13 can be expanded to sums
of products of the form $A_{1}G^{+}_{0}A_{2}$, of varying length. The
$A_{i}$’s are either $T^{(0)}$ or $V^{(\nu)}$ with $\nu=1,2,3$. For example,
the NLO correction in Eq. 11 reads

$\displaystyle T^{(1)}$
$\displaystyle=V^{(1)}+T^{(0)}G^{+}_{0}V^{(1)}+V^{(1)}G^{+}_{0}T^{(0)}$
$\displaystyle+T^{(0)}G^{+}_{0}V^{(1)}G^{+}_{0}T^{(0)}.$ (18)

Clearly, the fundamental matrix elements that need to be evaluated at sub-
leading orders are always of the form

$\bra{p^{\prime},\ell^{\prime}}A_{1}G^{+}_{0}A_{2}\ket{p,\ell},$ (19)

where we omit the $s$ and $j$ quantum numbers that are identical for the ket
and the bra. In Appendix B we show how to evaluate Eq. 19 using ordinary
matrix products and Gauss-Legendre quadrature. Longer products, e.g., of the
form $A_{1}G^{+}_{0}A_{2}G^{+}_{0}A_{3}$, are straightforwardly reduced to the
form in Eq. 19 by the associativity of matrix products. Knowing this, and the
distributive property with respect to addition, we can also reduce the
computational complexity of evaluating the perturbation series for $T$ by
computing and storing the composite operators $\Omega^{\dagger}_{-}$,
$\Omega_{+}$, and $G^{+}_{1}$.

For separable potentials of Yamaguchi type [61], both the distorted-wave
series and the LS equation can be solved analytically. We exploit this to
verify our numerical implementation and to inspect the stability of the
perturbative expansion. Numerical and analytical results for semi-realistic
and separable Yamaguchi potentials in the ${}^{1}S_{0}$ and
${}^{3}S_{1}\mathrm{-}^{3}D_{1}$ channels agree to at least single precision.

###  II.4 Calibrating the low-energy constants

Our focus in this work is to predict and analyze the description of $np$
scattering observables in MWPC and specifically the PC of Long and Yang. To
enable quantitative calculations, we calibrate the values of the unknown LECs
using the same approach as Long and Yang, i.e., by tuning the contact LECs to
achieve a good reproduction of the low-energy Nijmegen phase shifts [62] at
selected scattering energies.

Before discussing the details of the calibration, it is important to remember
that the order-by-order amplitudes

$T=T^{(0)}+T^{(1)}+T^{(2)}+\ldots$ (20)

are computed perturbatively and their sum is unitary only up to perturbative
corrections. To obtain real-valued phase shifts in the calibration of the LECs
we must compute phase shifts perturbatively by expanding the $np$ $S$-matrix
and matching to chiral orders, see Appendix C for details. If one instead
solves for the partial-wave $S$-matrix non-perturbatively from the order-by-
order sum of $T^{(\nu)}$ amplitudes, the corresponding phase shifts will have
a non-zero imaginary part that increases with scattering energy. Indeed,
Figure 2 shows phase shifts computed perturbatively and non-perturbatively in
the two channels ${}^{1}D_{2}$ and ${}^{3}D_{2}$. There are no LECs that need
to be calibrated in these channels at the orders considered in this work. The
imaginary part of the non-perturbative phase shift increases with scattering
energy. As that happens, the real part of the phase shift and the (real-
valued) perturbative phase shift differ progressively. This is consistent with
observations in Ref. [63].

Figure 2: $np$ scattering phase shifts in the ${}^{1}D_{2}$ (top row) and
${}^{3}D_{2}$ (bottom row) channels at NLO, N2LO, and N3LO using a momentum
cutoff $\Lambda=2500$ MeV. Phase shifts computed using the perturbative method
are shown with black solid lines. The red dashed and dot-dashed lines show the
real and imaginary parts, respectively, of the phase shift computed by summing
the $T$-matrix contribution and using the non-perturbative relation between
phase shifts and the $S$-matrix. The black dashed lines show phase shifts from
the Nijmegen analysis [62].

In the calibration of LECs, we do not account for uncertainties stemming from
the Nijmegen phase shifts or the truncation of the $\chi$EFT expansion. While
we are aware of the potential risk of overfitting in doing so, we opted for a
simple approach to establish a first quantitative potential and a baseline
understanding. The application of Bayesian inference methods [47, 48, 49] to
quantify the posterior probability distributions for the values of the LECs in
MWPC [38], though more robust, requires considerably more efforts. In this
work, we focus on studying the effectiveness of MWPC for realistic description
of elastic $np$ scattering.

The $T_{\mathrm{lab}}$ values of the Nijmegen phase shifts used as calibration
data are listed in Table 2 for each channel and order. The calibrated LECs up
to N3LO are compiled in Table 3 in Appendix A. We use a naming-convention
where capital letters $C,D,E,\ldots$ denote LECs with dimension MeV-2, MeV-4,
MeV${}^{-6},\ldots$, respectively.

Each LEC receives perturbative corrections at subsequent orders from where it
was first introduced. As an example, the LO LEC $C_{{}^{1}S_{0}}$ is expanded
into contributions

$C_{{}^{1}S_{0}}=C^{(0)}_{{}^{1}S_{0}}+C^{(1)}_{{}^{1}S_{0}}+C^{(2)}_{{}^{1}S_{0}}+\dots,$
(21)

where the superscript enumerates the perturbative correction and not the
chiral order. In the following we will exemplify the calibration procedure by
discussing in detail how we calibrated the LECs in the ${}^{1}S_{0}$ channel.

Table 2: Laboratory scattering energies $T_{\mathrm{lab}}$ (in MeV) of the Nijmegen phase shifts [62] used to calibrate the values of the LECs at each chiral order. In total, we employed 33 single-energy phase shifts—the same as the total number of contact LECs in the chiral expansion of Long and Yang up to N3LO. Channel | LO | NLO | N2LO | N3LO
---|---|---|---|---
${}^{1}S_{0}$ | 5 | 5, 25 | 5, 25, 50 | 5, 25, 50, 75
${}^{3}P_{0}$ | 25 | - | 25, 50 | 75, 100
${}^{1}P_{1}$ | - | - | 50 | 50
${}^{3}P_{1}$ | - | - | 50 | 50
${}^{3}S_{1}\mathrm{-}^{3}D_{1}$ | ${}^{3}S_{1}:30$ | - | ${}^{3}S_{1}:30,50.$ | ${}^{3}S_{1}:30,50.$
|  |  | $\epsilon_{1}:50$ | $\epsilon_{1}:50$
${}^{3}P_{2}\mathrm{-}^{3}F_{2}$ | ${}^{3}P_{2}:30$ | - | ${}^{3}P_{2}:30,50.$ | ${}^{3}P_{2}:30,50.$
|  |  | $\epsilon_{2}:50$ | $\epsilon_{2}:50$

At LO we calibrate the LEC $C^{(0)}_{{}^{1}S_{0}}$ such that the LO
${}^{1}S_{0}$ phase shift, $\delta^{(0)}$, reproduces the Nijmegen phase shift
at $T_{\mathrm{lab}}=5$ MeV. Two LECs are present in the ${}^{1}S_{0}$ channel
of the NLO potential: $D^{(0)}_{{}^{1}S_{0}}$ and $C^{(1)}_{{}^{1}S_{0}}$. The
latter is a perturbative correction to the LO LEC. These two LECs are
calibrated such that the LO phase shift plus the perturbative NLO correction,
i.e., $\delta^{(0)}+\delta^{(1)}$, reproduce the Nijmegen phase shifts at
$T_{\mathrm{lab}}=5$ and 25 MeV. The role of $C^{(1)}_{{}^{1}S_{0}}$ is to
ensure that the NLO correction vanishes for $T_{\mathrm{lab}}=5$ MeV. At N2LO
we have the LECs $C^{(2)}_{{}^{1}S_{0}},\ D^{(1)}_{{}^{1}S_{0}},\
E^{(0)}_{{}^{1}S_{0}}$ calibrated to phase shifts at energies
$T_{\mathrm{lab}}=5,25$ and 50 MeV. Finally, at N3LO the LECs
$C^{(3)}_{{}^{1}S_{0}},\ D^{(2)}_{{}^{1}S_{0}},\ E^{(1)}_{{}^{1}S_{0}},\
F^{(0)}_{{}^{1}S_{0}}$ are calibrated to reproduce the phase shifts at
$T_{\mathrm{lab}}=5,25,50$ and 75 MeV. An analogous scheme is employed for the
remaining partial waves and LECs. We calibrate all LECs for two different
momentum cutoffs: $\Lambda=500$ and 2500 MeV.

For the channels where OPE is perturbative there are no LECs present that need
to be calibrated. As a consistency check we compute and reproduce the
scattering phase shifts of Ref. [52]. Figure 3 shows our fit of the phase
shifts in the channels where OPE is non-perturbative. The bands indicate the
variation due to the two different cutoff values. There is an overall order-
by-order convergence in all channels up to around $T_{\mathrm{lab}}=100$ MeV
and we can reproduce the known results of [30, 32, 51]. The degree of cutoff
sensitivity varies notably among different channels. For instance, channels
like ${}^{1}P_{1}$ and ${}^{3}F_{2}$ show minimal sensitivity to the cutoff
value, while ${}^{3}P_{2}$ and $\epsilon_{1}$ demonstrate a more pronounced
dependency. The calibration in the ${}^{3}P_{0}$ channel was particularly
challenging at the higher chiral orders and the calibration energies needed to
be shifted to relatively high values at N3LO, as seen in Table 2.

Figure 3: Phase shifts in the channels where OPE is non-perturbative and the
amplitudes are computed using full distorted-wave perturbation theory. The
bands indicate the envelope of the variation due to the two different cutoff
values; 500 MeV (dashed line) and 2500 MeV (solid line). Note that LO and NLO
results coincide for all channels except ${}^{1}S_{0}$, which is why the blue
NLO band appears to be missing in several panels. The black solid lines show
phase shifts from the Nijmegen analysis [62] and the diamond markers indicate
the calibration data at $T_{\mathrm{lab}}$ values from Table 2.

##  III Neutron-Proton Scattering Observables

Here we predict selected $np$ scattering observables up to
$T_{\mathrm{lab}}\approx 100$ MeV using the potentials that were defined and
calibrated in Section II. We compute scattering observables from the partial-
wave amplitudes by first constructing the spin-scattering matrix, $M$, by [64,
65, 56]

$\displaystyle M^{s}_{m^{\prime}_{s}m_{s}}$
$\displaystyle(p_{0},\theta_{\mathrm{cm}},\phi)=\frac{\sqrt{4\pi}}{2ip_{0}}\sum_{j,\ell,\ell^{\prime}}i^{\ell-\ell^{\prime}}(2j+1)\sqrt{2\ell+1}$
$\displaystyle\times\begin{pmatrix}\ell^{\prime}&s&j\\\
m_{s}-m^{\prime}_{s}&m^{\prime}_{s}&-m_{s}\end{pmatrix}\begin{pmatrix}\ell&s&j\\\
0&m_{s}&-m_{s}\end{pmatrix}$ (22) $\displaystyle\times
Y^{\ell^{\prime}}_{m_{s}-m^{\prime}_{s}}(\theta_{\mathrm{cm}},\phi)\left(S^{(\nu)js}_{\ell^{\prime}\ell}(p_{0},p_{0})-\delta_{\ell^{\prime}\ell}\right).$

The angles $\theta_{\mathrm{cm}}\in[0,\pi]$ and $\phi\in[0,2\pi]$ are the
polar and azimuthal scattering angles, respectively where the latter is set to
zero by cylindrical symmetry. The on-shell scattering momentum, $p_{0}$, is
given from the laboratory scattering energy $T_{\mathrm{lab}}$ using Eq. 44 in
Appendix B. We compute $S^{(\nu)js}_{\ell^{\prime}\ell}(p_{0},p_{0})$, i.e.,
the $S-$matrix for a potential up to some chiral order $\nu$, by summing the
perturbatively computed $T$-matrix amplitudes to order $\nu$. Using the
conventions applied in this work, the partial-wave relation between the on-
shell $S$\- and $T$-matrix elements is thus given by

$\displaystyle
S^{(\nu)js}_{\ell^{\prime}\ell}(p_{0},p_{0})=\delta_{\ell^{\prime}\ell}-i\pi
m_{N}p_{0}$
$\displaystyle\times\left[T^{(0)js}_{\ell^{\prime}\ell}(p_{0},p_{0})+\dots+T^{(\nu)js}_{\ell^{\prime}\ell}(p_{0},p_{0})\right].$
(23)

We focus our discussion on the differential $np$ scattering cross section and
two selected polarizations, and calculate these from the spin-scattering
matrix as

$\displaystyle\frac{d\sigma}{d\Omega}$
$\displaystyle=\frac{1}{4}\operatorname{Tr}{MM^{\dagger}}$ (24)
$\displaystyle\frac{d\sigma}{d\Omega}\times P_{b}$
$\displaystyle=\frac{1}{4}\operatorname{Tr}{M\bm{\sigma}_{1n}M^{\dagger}}$
(25) $\displaystyle\frac{d\sigma}{d\Omega}\times A_{yy}$
$\displaystyle=\frac{1}{4}\operatorname{Tr}{M\bm{\sigma}_{1n}\bm{\sigma}_{2n}M^{\dagger}}$
(26)

where $\bm{\sigma}_{in}\equiv\bm{\sigma}_{i}\cdot\hat{\bm{n}}$ for nucleon
$i$, $\bm{\sigma}_{i}$ is the Pauli spin matrices, and $\hat{\bm{n}}$ is
normal to the scattering plane.

Figure 4: Selection of $np$ scattering observables in the energy interval
$T_{\mathrm{lab}}=10$ to 100 MeV. Experimental data from Refs. [66, 67]. The
bands indicate cutoff variation in the same way as in Fig. 3.

Figure 4 shows our prediction for these scattering observables in the energy
range $T_{\mathrm{lab}}=10$ to 100 MeV for the two cutoffs $\Lambda=500$ MeV
and $\Lambda=2500$ MeV. For the lower scattering energies
($T_{\mathrm{lab}}\lesssim 60$ MeV) we observe an order-by-order improvement
for all considered observables. Interestingly, the N3LO predictions do not
always perform better, but in general performs at least as well as N2LO.
Indeed, for $T_{\text{lab}}\approx$ 100 MeV (rightmost panels of Fig. 4), it
appears that the order-by-order improvement in the predictions of the
differential cross section and $P_{b}$ polarization deteriorates and N2LO can
perform better than N3LO. This effect is visible also at the level of phase
shifts shown in Fig. 3. It is not clear at the moment if this is due to
overfitting and (or) an underlying issue with the MWPC that we employ. Our
N3LO predictions are certainly influenced by the adopted values of sub-leading
$\pi N$ LECs [54]. Calculations of other scattering cross observables show
that the order-by-order convergence demonstrated in Fig. 4 is representative
for all elastic $np$ scattering observables in the PC by Long and Yang. Two-
pion exchange is clearly important for achieving a realistic description of
scattering observables with $T_{\mathrm{lab}}\lesssim 100$ MeV.

The total cross section can be straightforwardly computed from the
differential cross section as

$\sigma_{\mathrm{tot}}(p_{0})=2\pi\int_{-1}^{1}d(\cos\theta_{\mathrm{cm}})\
\frac{d\sigma}{d\Omega}(p_{0},\theta_{\mathrm{cm}}),$ (27)

and predictions for scattering energies up to $T_{\mathrm{lab}}=150$ MeV are
shown in Fig. 5. Also for this obvservable, the agreement with experimental
data typically improves order-by-order, at least up to N2LO. The improvement
of N3LO over N2LO is not obvious. At very low energies, the higher-order
predictions for the total cross section are much better than the lower-order
predictions. This result is somewhat peculiar for a low-energy EFT and likely
due to overfitting at the phase shift level. For $T_{\mathrm{lab}}\gtrsim 100$
MeV, roughly corresponding to 220 MeV relative momentum, the agreement with
data even deteriorates at N3LO. This is analogous to what was found for the
angular-differential observables shown in Fig. 4 and consistent with the
observation in Fig. 3 that the phase shifts at N3LO might suffer from
overfitting at the higher energies. Alternatively, the observed decline in
predictive power might indicate the presence of an additional mass scale at
200-300 MeV. Thus, it will be very interesting to study the effects of
accounting for the $\Delta(1232)$-isobar in two-pion exchange in this MWPC.

Figure 5: Total $np$ cross sections computed by integrating the differential
cross sections (27). Panel $(a)$ shows cross sections for a large interval of
scattering energies, $T_{\mathrm{lab}}=5\text{--}150$ MeV. Panels $(b)$ and
$(c)$ expand results at low- and high-energy intervals, respectively. The
bands indicate cutoff variation as in Fig. 3. Experimental data from Refs.
[66, 67].

Next, we analyze how the perturbative breaking of unitarity in $\chi$EFT
affects the predictions of total cross sections. Indeed, the computation of
$S$-matrix elements using Eq. 23, where the order-by-order contributions of
the scattering amplitudes are summed directly to the $S$-matrix, leads to a
perturbative breaking of unitarity. In contrast, amplitudes computed non-
perturbatively, i.e., when the potential terms are summed before solving for
the scattering amplitude (as is done in WPC), are unitary by construction. In
this case, the probability flux in the scattering process is also conserved
exactly and the optical theorem can be safely used to compute the total cross
section as, e.g.,

$\sigma_{\mathrm{tot}}(p_{0})=\frac{2\pi}{p_{0}}\operatorname{Im}\left[a(\theta_{\mathrm{cm}}=0)+b(\theta_{\mathrm{cm}}=0)\right],$
(28)

where $a(\theta_{\mathrm{cm}})$ and $b(\theta_{\mathrm{cm}})$ are Saclay-
amplitudes computed from the $M$-matrix [68].

We use the difference between total cross sections calculated using Eq. 27 and
Eq. 28 to measure the effects of unitarity breaking. In Fig. 6 we show the
relative difference between the cross sections computed using exact
integration and the optical theorem as a function of scattering energy. The
figure demonstrates how unitarity is restored perturbatively as we go to
higher chiral orders. Indeed, the relative difference between the two cross
section calculations is limited to 10% for scattering energies up to 40 MeV at
NLO, 70 MeV at N2LO, and 120 MeV at N3LO, respectively. The bands in the
figure reflect differences coming from using two cutoff values 500 MeV and
2500 MeV. The bands for NLO and N2LO increase smoothly with the scattering
energy. The band at N3LO shows an artifact from the two different calculations
for $\Lambda=2500$ MeV intersecting at some energies leading to very small
relative errors. We also note that the cutoff dependencies for the N2LO and
N3LO calculations do not vanish as the scattering energy approaches zero.

Figure 6: The relative difference between total $np$ cross sections ($\sigma$)
computed by integrating of the differential cross section (27) and the optical
theorem (28). The bands indicate cutoff variation as in Fig. 3. The color
coding for the orders is the same as Fig. 3. The horizontal dashed line marks
a $10$% difference.

We can also discuss this result in terms of the EFT truncation error. For a
given chiral order, we argue that the results from the two different cross
section calculations should not become significantly different until we reach
an energy where the next (omitted) order in the chiral low-energy expansion
becomes relevant. This should correspond to the scattering energy for which
the truncation error is significant. Breaking unitarity implies that the norm
of the partial-wave $S$-matrix in Eq. 23 deviates from unity as
$\left(S^{(\nu)}\right)^{\dagger}S^{(\nu)}=1-\mathcal{C}(Q/\Lambda_{b})^{\nu+1}$,
where we also expect $\mathcal{C}$ to be of natural size. This scaling of
unitarity breaking should be revisited when probability distributions of the
LEC values and the hyperparameters of the EFT truncation error have been
inferred using a Bayesian approach.

##  IV Summary and outlook

This work presents a comprehensive analysis of $np$ scattering observables
(cross sections and polarizations) utilizing an RG-invariant formulation of
$\chi$EFT by Long and Yang. We calibrated the LECs by reproducing Nijmegen
phase shifts at specific scattering energies, and carried out calculation up
to N3LO for two values of the momentum-space cutoffs, $500$ MeV and $2500$
MeV. The PC that we employed is fairly representative of a broad class of
MWPCs in which corrections beyond LO, based on one-pion exchange, are included
perturbatively and the short-range contact potential incorporates counterterms
promoted to renormalize the long-range pion contributions to the scattering
amplitudes. A key result of this paper was a quantitative demonstration that
RG-invariant $\chi$EFT exhibits a steady order-by-order convergence in the
description of scattering observables, starting already at LO. A second key
result was the realistic reproduction of experimental scattering data in an
energy range up to $T_{\mathrm{lab}}=100$ MeV at N2LO. We also found that N3LO
predictions do not always improve over N2LO.

A perturbative approach exposes the deficiencies of any PC, not only the
possible lack of RG-independence. In fact, using a perturbative approach we
found that the accuracy of our N3LO predictions for the total $np$ cross
section declines as one approaches $T_{\mathrm{lab}}\gtrsim 100$ MeV. This
corresponds to a relative scattering momentum of 220 MeV and might suggest the
presence of an additional mass scale at 200–300 MeV. This finding is in
accordance with the known mass splitting between the nucleon and the
$\Delta$(1232) resonance, but is markedly lower than conventional estimates of
the breakdown scale of $\chi$EFT residing in the vicinity of the $\rho$-meson
mass. The latter estimate has also been corroborated in a Bayesian study of
non-perturbative WPC predictions of nucleon-nucleon scattering observables
[69].

Based on our comparison of perturbative and non-perturbative calculations of
phase shifts, we speculated that the magnitudes of the imaginary component of
the non-perturbative phase shift and the $\chi$EFT truncation error are
linked. We also investigated the breaking of unitarity at the level of total
$np$ cross sections. The connection between perturbative unitarity breaking
and the truncation error deserves further attention.

Future work will focus on quantifying posterior probability distributions for
the LECs and the EFT truncation error, making predictions beyond the two-
nucleon system, and the effects of including the $\Delta$(1232) resonance in
the two-pion exchange potential. Fast and accurate emulators [70], adapted to
perturbative computations, will likely be essential for rigorous testing of
RG-invariant $\chi$EFT against nuclear data and to address critical questions
regarding, e.g., the construction of LO, the importance of promoting higher-
order pion exchanges and many-nucleon forces as one increases the mass number,
and the level of fine-tuning in $\chi$EFT.

###### Acknowledgements.

O.T. thanks C.-J. Yang, B. Long, and R. Peng for helpful discussions and for
providing detailed benchmarks. The authors also thank Daniel Phillips for
feedback on a draft version of the manuscript. 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.

## References

  * Ekström  _et al._ [2023] A. Ekström, C. Forssén, G. Hagen, G. R. Jansen, W. Jiang, and T. Papenbrock, What is ab initio in nuclear theory?, Front. Phys. 11, 1129094 (2023), arXiv:2212.11064 [nucl-th] .
  * Hergert [2020] H. Hergert, A Guided Tour of $ab$ $initio$ Nuclear Many-Body Theory, Front. in Phys. 8, 379 (2020), arXiv:2008.05061 [nucl-th] .
  * Epelbaum  _et al._ [2009] E. Epelbaum, H.-W. Hammer, and U.-G. Meissner, Modern Theory of Nuclear Forces, Rev. Mod. Phys. 81, 1773 (2009), arXiv:0811.1338 [nucl-th] .
  * Machleidt and Entem [2011] R. Machleidt and D. R. Entem, Chiral effective field theory and nuclear forces, Phys. Rept. 503, 1 (2011), arXiv:1105.2919 [nucl-th] .
  * Hammer  _et al._ [2020] H. W. Hammer, S. König, and U. van Kolck, Nuclear effective field theory: status and perspectives, Rev. Mod. Phys. 92, 025004 (2020), arXiv:1906.12122 [nucl-th] .
  * Weinberg [1990] S. Weinberg, Nuclear forces from chiral Lagrangians, Phys. Lett. B 251, 288 (1990).
  * Weinberg [1991] S. Weinberg, Effective chiral Lagrangians for nucleon - pion interactions and nuclear forces, Nucl. Phys. B 363, 3 (1991).
  * Entem and Machleidt [2003] D. R. Entem and R. Machleidt, Accurate charge dependent nucleon nucleon potential at fourth order of chiral perturbation theory, Phys. Rev. C 68, 041001 (2003), arXiv:nucl-th/0304018 .
  * Ekström  _et al._ [2013] A. Ekström _et al._ , Optimized Chiral Nucleon-Nucleon Interaction at Next-to-Next-to-Leading Order, Phys. Rev. Lett. 110, 192502 (2013), arXiv:1303.4674 [nucl-th] .
  * Gezerlis  _et al._ [2013] A. Gezerlis, I. Tews, E. Epelbaum, S. Gandolfi, K. Hebeler, A. Nogga, and A. Schwenk, Quantum Monte Carlo Calculations with Chiral Effective Field Theory Interactions, Phys. Rev. Lett. 111, 032501 (2013), arXiv:1303.6243 [nucl-th] .
  * Piarulli  _et al._ [2015] M. Piarulli, L. Girlanda, R. Schiavilla, R. Navarro Pérez, J. E. Amaro, and E. Ruiz Arriola, Minimally nonlocal nucleon-nucleon potentials with chiral two-pion exchange including $\Delta$ resonances, Phys. Rev. C 91, 024003 (2015), arXiv:1412.6446 [nucl-th] .
  * Carlsson  _et al._ [2016] B. D. Carlsson, A. Ekström, C. Forssén, D. F. Strömberg, G. R. Jansen, O. Lilja, M. Lindby, B. A. Mattsson, and K. A. Wendt, Uncertainty analysis and order-by-order optimization of chiral nuclear interactions, Phys. Rev. X 6, 011019 (2016), arXiv:1506.02466 [nucl-th] .
  * Ekström  _et al._ [2015] A. Ekström, G. R. Jansen, K. A. Wendt, G. Hagen, T. Papenbrock, B. D. Carlsson, C. Forssén, M. Hjorth-Jensen, P. Navrátil, and W. Nazarewicz, Accurate nuclear radii and binding energies from a chiral interaction, Phys. Rev. C 91, 051301 (2015), arXiv:1502.04682 [nucl-th] .
  * Jiang  _et al._ [2020] W. G. Jiang, A. Ekström, C. Forssén, G. Hagen, G. R. Jansen, and T. Papenbrock, Accurate bulk properties of nuclei from $A=2$ to $\infty$ from potentials with $\Delta$ isobars, Phys. Rev. C 102, 054301 (2020), arXiv:2006.16774 [nucl-th] .
  * Reinert  _et al._ [2018] P. Reinert, H. Krebs, and E. Epelbaum, Semilocal momentum-space regularized chiral two-nucleon potentials up to fifth order, Eur. Phys. J. A 54, 86 (2018), arXiv:1711.08821 [nucl-th] .
  * Entem  _et al._ [2017] D. R. Entem, R. Machleidt, and Y. Nosyk, High-quality two-nucleon potentials up to fifth order of the chiral expansion, Phys. Rev. C 96, 024004 (2017), arXiv:1703.05454 [nucl-th] .
  * Epelbaum  _et al._ [2020] E. Epelbaum, H. Krebs, and P. Reinert, High-precision nuclear forces from chiral EFT: State-of-the-art, challenges and outlook, Front. in Phys. 8, 98 (2020), arXiv:1911.11875 [nucl-th] .
  * Tews  _et al._ [2020] I. Tews, Z. Davoudi, A. Ekström, J. D. Holt, and J. E. Lynn, New Ideas in Constraining Nuclear Forces, J. Phys. G 47, 103001 (2020), arXiv:2001.03334 [nucl-th] .
  * Tews  _et al._ [2022] I. Tews _et al._ , Nuclear Forces for Precision Nuclear Physics: A Collection of Perspectives, Few Body Syst. 63, 67 (2022), arXiv:2202.01105 [nucl-th] .
  * van Kolck [2020] U. van Kolck, The Problem of Renormalization of Chiral Nuclear Forces, Front. in Phys. 8, 79 (2020), arXiv:2003.06721 [nucl-th] .
  * Nogga  _et al._ [2005] A. Nogga, R. G. E. Timmermans, and U. van Kolck, Renormalization of one-pion exchange and power counting, Phys. Rev. C 72, 054006 (2005), arXiv:nucl-th/0506005 .
  * Frank  _et al._ [1971] W. Frank, D. J. Land, and R. M. Spector, Singular potentials, Rev. Mod. Phys. 43, 36 (1971).
  * Pavon Valderrama and Ruiz Arriola [2006] M. Pavon Valderrama and E. Ruiz Arriola, Renormalization of NN interaction with chiral two pion exchange potential: Non-central phases, Phys. Rev. C 74, 064004 (2006), [Erratum: Phys.Rev.C 75, 059905 (2007)], arXiv:nucl-th/0507075 .
  * Pavon Valderrama and Ruiz Arriola [2005] M. Pavon Valderrama and E. Ruiz Arriola, Renormalization of the deuteron with one pion exchange, Phys. Rev. C 72, 054002 (2005), arXiv:nucl-th/0504067 .
  * Pavon Valderrama [2011] M. Pavon Valderrama, Perturbative Renormalizability of Chiral Two Pion Exchange in Nucleon-Nucleon Scattering: P- and D-waves, Phys. Rev. C 84, 064002 (2011), arXiv:1108.0872 [nucl-th] .
  * Long [2013] B. Long, Improved convergence of chiral effective field theory for 1S0 of NN scattering, Phys. Rev. C 88, 014002 (2013), arXiv:1304.7382 [nucl-th] .
  * Pavón Valderrama  _et al._ [2017] M. Pavón Valderrama, M. Sánchez Sánchez, C. J. Yang, B. Long, J. Carbonell, and U. van Kolck, Power Counting in Peripheral Partial Waves: The Singlet Channels, Phys. Rev. C 95, 054001 (2017), arXiv:1611.10175 [nucl-th] .
  * Valderrama [2011] M. P. Valderrama, Perturbative renormalizability of chiral two pion exchange in nucleon-nucleon scattering, Phys. Rev. C 83, 024003 (2011), arXiv:0912.0699 [nucl-th] .
  * Birse [2006] M. C. Birse, Power counting with one-pion exchange, Phys. Rev. C 74, 014003 (2006), arXiv:nucl-th/0507077 .
  * Long and Yang [2012a] B. Long and C. J. Yang, Short-range nuclear forces in singlet channels, Phys. Rev. C 86, 024001 (2012a), arXiv:1202.4053 [nucl-th] .
  * Sánchez Sánchez  _et al._ [2018] M. Sánchez Sánchez, C. J. Yang, B. Long, and U. van Kolck, Two-nucleon ${}^{1}S_{0}$ amplitude zero in chiral effective field theory, Phys. Rev. C 97, 024001 (2018), arXiv:1704.08524 [nucl-th] .
  * Long and Yang [2012b] B. Long and C.-J. Yang, Renormalizing chiral nuclear forces: Triplet channels, Phys. Rev. C 85, 034002 (2012b).
  * Yang [2016] C. J. Yang, Chiral potential renormalized in harmonic-oscillator space, Phys. Rev. C 94, 064004 (2016), arXiv:1610.01350 [nucl-th] .
  * Mishra  _et al._ [2022] C. Mishra, A. Ekström, G. Hagen, T. Papenbrock, and L. Platter, Two-pion exchange as a leading-order contribution in chiral effective field theory, Phys. Rev. C 106, 024004 (2022), arXiv:2111.15515 [nucl-th] .
  * Peng  _et al._ [2022] R. Peng, S. Lyu, S. König, and B. Long, Constructing chiral effective field theory around unnatural leading-order interactions, Phys. Rev. C 105, 054002 (2022), arXiv:2112.00947 [nucl-th] .
  * Song  _et al._ [2017] Y.-H. Song, R. Lazauskas, and U. van Kolck, Triton binding energy and neutron-deuteron scattering up to next-to-leading order in chiral effective field theory, Phys. Rev. C 96, 024002 (2017), [Erratum: Phys.Rev.C 100, 019901 (2019)], arXiv:1612.09090 [nucl-th] .
  * Yang  _et al._ [2021] C. J. Yang, A. Ekström, C. Forssén, and G. Hagen, Power counting in chiral effective field theory and nuclear binding, Phys. Rev. C 103, 054304 (2021), arXiv:2011.11584 [nucl-th] .
  * Thim  _et al._ [2023] O. Thim, E. May, A. Ekström, and C. Forssén, Bayesian analysis of chiral effective field theory at leading order in a modified Weinberg power counting approach, Phys. Rev. C 108, 054002 (2023), arXiv:2302.12624 [nucl-th] .
  * Yang  _et al._ [2023] C. J. Yang, A. Ekström, C. Forssén, G. Hagen, G. Rupak, and U. van Kolck, The importance of few-nucleon forces in chiral effective field theory, Eur. Phys. J. A 59, 233 (2023), arXiv:2109.13303 [nucl-th] .
  * Long and van Kolck [2008] B. Long and U. van Kolck, Renormalization of Singular Potentials and Power Counting, Annals Phys. 323, 1304 (2008), arXiv:0707.4325 [quant-ph] .
  * Epelbaum and Meissner [2013] E. Epelbaum and U. G. Meissner, On the Renormalization of the One-Pion Exchange Potential and the Consistency of Weinberg‘s Power Counting, Few Body Syst. 54, 2175 (2013), arXiv:nucl-th/0609037 .
  * Epelbaum and Gegelia [2009] E. Epelbaum and J. Gegelia, Regularization, renormalization and ’peratization’ in effective field theory for two nucleons, Eur. Phys. J. A 41, 341 (2009), arXiv:0906.3822 [nucl-th] .
  * Epelbaum  _et al._ [2018] E. Epelbaum, A. M. Gasparyan, J. Gegelia, and U.-G. Meißner, How (not) to renormalize integral equations with singular potentials in effective field theory, Eur. Phys. J. A 54, 186 (2018), arXiv:1810.02646 [nucl-th] .
  * Gasparyan and Epelbaum [2022] A. M. Gasparyan and E. Epelbaum, Is the RG-invariant EFT for few-nucleon systems cutoff independent?, arxiv (2022), arXiv:2210.16225 [nucl-th] .
  * Long and Yang [2011a] B. Long and C.-J. Yang, Renormalizing chiral nuclear forces: A case study of ${}^{3}\phantom{\rule{-1.60004pt}{0.0pt}}{P}_{0}$, Phys. Rev. C 84, 057001 (2011a).
  * Furnstahl  _et al._ [2015] R. J. Furnstahl, N. Klco, D. R. Phillips, and S. Wesolowski, Quantifying truncation errors in effective field theory, Phys. Rev. C 92, 024005 (2015), arXiv:1506.01343 [nucl-th] .
  * Wesolowski  _et al._ [2019] S. Wesolowski, R. J. Furnstahl, J. A. Melendez, and D. R. Phillips, Exploring bayesian parameter estimation for chiral effective field theory using nucleon-nucleon phase shifts, J. Phys. G 46, 045102 (2019).
  * Wesolowski  _et al._ [2021] S. Wesolowski, I. Svensson, A. Ekström, C. Forssén, R. J. Furnstahl, J. A. Melendez, and D. R. Phillips, Rigorous constraints on three-nucleon forces in chiral effective field theory from fast and accurate calculations of few-body observables, Phys. Rev. C 104, 064001 (2021), arXiv:2104.04441 [nucl-th] .
  * Svensson  _et al._ [2022] I. Svensson, A. Ekström, and C. Forssén, Bayesian parameter estimation in chiral effective field theory using the hamiltonian monte carlo method, Phys. Rev. C 105, 014004 (2022).
  * Svensson  _et al._ [2023] I. Svensson, A. Ekström, and C. Forssén, Inference of the low-energy constants in delta-full chiral effective field theory including a correlated truncation error, (2023), arXiv:2304.02004 [nucl-th] .
  * Long and Yang [2011b] B. Long and C. J. Yang, Renormalizing chiral nuclear forces: a case study of 3P0, Phys. Rev. C 84, 057001 (2011b), arXiv:1108.0985 [nucl-th] .
  * Wu and Long [2019] S. Wu and B. Long, Perturbative $nn$ scattering in chiral effective field theory, Phys. Rev. C 99, 024003 (2019).
  * Peng  _et al._ [2020] R. Peng, S. Lyu, and B. Long, Perturbative chiral nucleon–nucleon potential for the ${}^{3}P_{0}$ partial wave, Commun. Theor. Phys. 72, 095301 (2020), arXiv:2011.13186 [nucl-th] .
  * Siemens  _et al._ [2017] D. Siemens, J. Ruiz de Elvira, E. Epelbaum, M. Hoferichter, H. Krebs, B. Kubis, and U. G. Meißner, Reconciling threshold and subthreshold expansions for pion–nucleon scattering, Phys. Lett. B 770, 27 (2017), arXiv:1610.08978 [nucl-th] .
  * Barford and Birse [2003] T. Barford and M. C. Birse, A Renormalization group approach to two-body scattering in the presence of long range forces, Phys. Rev. C 67, 064006 (2003), arXiv:hep-ph/0206146 .
  * Newton [1982] R. G. Newton,  _Scattering theory of waves and particles_ (Springer-Verlag New York, Inc., 175 Fifth Avenue, New York, New York, 10010, U.S.A., 1982).
  * Hussein and Canto [2012] M. S. Hussein and L. F. Canto,  _Scattering Theory Of Molecules, Atoms And Nuclei_ (World Scientific Publishing Company, Singapore, SINGAPORE, 2012).
  * Griesshammer [2022] H. W. Griesshammer, What Can Possibly Go Wrong?, Few Body Syst. 63, 44 (2022), arXiv:2111.00930 [nucl-th] .
  * Vanasse [2013] J. Vanasse, Fully Perturbative Calculation of $nd$ Scattering to Next-to-next-to-leading-order, Phys. Rev. C 88, 044001 (2013), arXiv:1305.0283 [nucl-th] .
  * Erkelenz  _et al._ [1971] K. Erkelenz, R. Alzetta, and K. Holinde, Momentum space calculations and helicity formalism in nuclear physics, Nucl. Phys. A 176, 413 (1971).
  * Yamaguchi [1954] Y. Yamaguchi, Two nucleon problem when the potential is nonlocal but separable. 1., Phys. Rev. 95, 1628 (1954).
  * Stoks  _et al._ [1993] V. G. J. Stoks, R. A. M. Klomp, M. C. M. Rentmeester, and J. J. de Swart, Partial wave analaysis of all nucleon-nucleon scattering data below 350-MeV, Phys. Rev. C 48, 792 (1993).
  * Odell  _et al._ [2023] D. Odell, D. R. Phillips, and U. van Kolck, Effective Field Theory for the Bound States and Scattering of a Heavy Charged Particle and a Neutral Atom, (2023), arXiv:2307.13103 [nucl-th] .
  * Blatt and Biedenharn [1952] J. M. Blatt and L. C. Biedenharn, The Angular Distribution of Scattering and Reaction Cross Sections, Rev. Mod. Phys. 24, 258 (1952).
  * Glöckle [1983] W. Glöckle,  _The Quantum Mechanical Few-body Problem_ (Springer-Verlag, Berlin Heidelberg, 1983).
  * Navarro Pérez  _et al._ [2013] R. Navarro Pérez, J. E. Amaro, and E. Ruiz Arriola, Partial-wave analysis of nucleon-nucleon scattering below the pion-production threshold, Phys. Rev. C 88, 024002 (2013).
  * Pérez  _et al._ [2013] R. N. Pérez, J. E. Amaro, and E. R. Arriola, Coarse-grained potential analysis of neutron-proton and proton-proton scattering below the pion production threshold, Phys. Rev. C 88, 064002 (2013).
  * Bystricky  _et al._ [1978] J. Bystricky, F. Lehar, and P. Winternitz, Formalism of Nucleon-Nucleon Elastic Scattering Experiments, J. Phys. (France) 39, 1 (1978).
  * Melendez  _et al._ [2017] J. A. Melendez, S. Wesolowski, and R. J. Furnstahl, Bayesian truncation errors in chiral effective field theory: nucleon-nucleon observables, Phys. Rev. C 96, 024003 (2017), arXiv:1704.03308 [nucl-th] .
  * Duguet  _et al._ [2023] T. Duguet, A. Ekström, R. J. Furnstahl, S. König, and D. Lee, Eigenvector Continuation and Projection-Based Emulators, (2023), arXiv:2310.19419 [nucl-th] .
  * Epelbaum  _et al._ [1998] E. Epelbaum, W. Gloeckle, and U.-G. Meissner, Nuclear forces from chiral Lagrangians using the method of unitary transformation. 1. Formalism, Nucl. Phys. A 637, 107 (1998), arXiv:nucl-th/9801064 .
  * Epelbaum  _et al._ [2000] E. Epelbaum, W. Gloeckle, and U.-G. Meissner, Nuclear forces from chiral Lagrangians using the method of unitary transformation. 2. The two nucleon system, Nucl. Phys. A 671, 295 (2000), arXiv:nucl-th/9910064 .
  * Haftel and Tabakin [1970] M. I. Haftel and F. Tabakin, Nuclear saturation and the smoothness of nucleon-nucleon potentials, Nucl. Phys. A 158, 1 (1970).
  * Landau [1990] R. H. Landau,  _Quantum mechanics. Vol. 2: A second course in quantum theory_ (1990).
  * Hoppe  _et al._ [2017] J. Hoppe, C. Drischler, R. J. Furnstahl, K. Hebeler, and A. Schwenk, Weinberg eigenvalues for chiral nucleon-nucleon interactions, Phys. Rev. C 96, 054002 (2017), arXiv:1707.06438 [nucl-th] .
  * Stapp  _et al._ [1957] H. P. Stapp, T. J. Ypsilantis, and N. Metropolis, Phase shift analysis of 310-MeV proton proton scattering experiments, Phys. Rev. 105, 302 (1957).

##  Appendix A Nuclear potentials in the Long and Yang power counting

The orders at which potentials appear in the Long and Yang PC in channels
where OPE is treated non-perturbatively are shown in Table 1. Similarly, for
the channels where OPE is treated perturbatively, we follow the PC of Ref.
[52] also shown in Table 1. In this appendix, we list the expressions for the
potentials appearing in Table 1. The potential contributions will be listed
using the following decomposition convention [4]

$\displaystyle V(\bm{p}^{\prime},\bm{p})$
$\displaystyle=V_{C}+\bm{\tau}_{1}\cdot\bm{\tau}_{2}W_{C}$ (29)
$\displaystyle+\left[V_{S}+\bm{\tau}_{1}\cdot\bm{\tau}_{2}W_{S}\right]\bm{\sigma}_{1}\cdot\bm{\sigma}_{2}+$
$\displaystyle+\left[V_{LS}+\bm{\tau}_{1}\cdot\bm{\tau}_{2}W_{LS}\right](-i\bm{S}\cdot\left(\bm{q}\times\bm{k})\right)$
$\displaystyle+\left[V_{T}+\bm{\tau}_{1}\cdot\bm{\tau}_{2}W_{T}\right]\bm{\sigma}_{1}\cdot\bm{q}\bm{\sigma}_{2}\cdot\bm{q}$
$\displaystyle+\left[V_{\sigma L}+\bm{\tau}_{1}\cdot\bm{\tau}_{2}W_{\sigma
L}\right]\bm{\sigma}_{1}\cdot(\bm{q}\times\bm{k})\bm{\sigma}_{2}\cdot(\bm{q}\times\bm{k}),$

where

$\bm{q}=\bm{p}-\bm{p}^{\prime},\quad\bm{k}=\frac{1}{2}\left(\bm{p}+\bm{p}^{\prime}\right),\quad\bm{S}=\frac{1}{2}\left(\bm{\sigma}_{1}+\bm{\sigma}_{2}\right)$
(30)

and $\bm{\sigma}_{i}$ denotes the Pauli spin matrix for the respective
nucleon.

The one-pion exchange potential takes the form

$\displaystyle V^{(0)}_{1\pi}$
$\displaystyle=\left(\bm{\tau}_{1}\cdot\bm{\tau}_{2}\right)\left(\bm{\sigma}_{1}\cdot\bm{q}\bm{\sigma}_{2}\cdot\bm{q}\right)W_{T},$
(31) $\displaystyle W_{T}$
$\displaystyle=-\left(\frac{g_{A}}{4f_{\pi}}\right)^{2}\frac{1}{q^{2}+m^{2}_{\pi}},$
(32)

where $g_{A}=1.29$ is the axial coupling, $f_{\pi}=92.1$ MeV the pion decay
constant, $m_{\pi}=138.039$ MeV is the average pion mass and $q=|\bm{q}|$. For
the two-pion exchange potentials, we employ expressions computed with
dimensional regularization (DR). The leading two-pion exchange potential takes
the form [71, 72, 4]

$\displaystyle V^{(2)}_{2\pi}$
$\displaystyle=\bm{\tau}_{1}\cdot\bm{\tau}_{2}W_{C}+\bm{\sigma}_{1}\cdot\bm{\sigma}_{2}V_{S}+\bm{\sigma}_{1}\cdot\bm{q}\bm{\sigma}_{2}\cdot\bm{q}V_{T},$
(33) $\displaystyle W_{C}$
$\displaystyle=-\frac{L(q)}{384\pi^{2}f^{4}_{\pi}}\Bigg{[}4m^{2}_{\pi}\left(5g^{4}_{A}-4g^{2}_{A}-1\right)+q^{2}\left(23g^{4}_{A}-10g^{2}_{A}-1\right)+\frac{48g^{4}_{A}m^{4}_{\pi}}{w^{2}}\Bigg{]},$
(34) $\displaystyle V_{S}$
$\displaystyle=\frac{3g^{4}_{A}L(q)q^{2}}{64\pi^{2}f^{4}_{\pi}},$ (35)
$\displaystyle V_{T}$
$\displaystyle=-\frac{1}{q^{2}}V_{S}=-\frac{3g^{4}_{A}L(q)}{64\pi^{2}f^{4}_{\pi}},$
(36)

with

$L(q)=\frac{w}{q}\ln\frac{w+q}{2m_{\pi}},\quad w=\sqrt{4m^{2}_{\pi}+q^{2}}.$
(37)

The sub-leading two-pion exchange potential takes the form of Eqs. (4.13) -
(4.20) in [4]. We apply the power counting
$\left(Q/m_{N}\right)=\left(Q/\Lambda_{b}\right)^{2}$ for $(1/m_{N})$
corrections, which means that all terms proportional to $1/m_{N}$ vanish at
order $\left(Q/\Lambda_{b}\right)^{3}$ (N3LO). The non-zero contributions read

$\displaystyle V^{(3)}_{2\pi}$
$\displaystyle=V_{C}+\left(\bm{\tau}_{1}\cdot\bm{\tau}_{2}\right)\left(\bm{\sigma}_{1}\cdot\bm{q}\bm{\sigma}_{2}\cdot\bm{q}\right)W_{T},$
(38) $\displaystyle V_{C}$ $\displaystyle=-\frac{3g^{2}_{A}}{16\pi
f^{4}_{\pi}}\Big{[}2m^{2}_{\pi}(2c_{1}-c_{3})-q^{2}c_{3}\Big{]}\tilde{w}^{2}A(q),$
(39) $\displaystyle W_{T}$
$\displaystyle=-\frac{1}{q^{2}}W_{S}=-\frac{g^{2}_{A}A(q)}{32\pi
f^{4}_{\pi}}c_{4}w^{2},$ (40)

with

$A(q)=\frac{1}{2q}\arctan\frac{q}{2m_{\pi}},\quad\tilde{w}=\sqrt{2m^{2}_{\pi}+q^{2}}.$
(42)

For the $\pi N$ LECs $c_{1},c_{3},c_{4}$, appearing in $V^{(3)}_{2\pi}$, we
employ numerical values determined in a Roy-Steiner analysis at NLO:
$c_{1}=-0.74$ GeV-1, $c_{3}=-3.61$ GeV-1 and $c_{4}=2.44$ GeV-1 [54].

The potential contributions at each order in the channels where OPE is treated
non-perturbatively are listed in Table 3. We denote counterterms in coupled
channels by a $2\times 2$ matrix representing $\ell^{\prime}=j\mp 1$ (rows)
and $\ell=j\mp 1$ (columns). Table 3 expands upon Table I in Ref. [30] to also
explicitly show the perturbative corrections to LECs present at each order.
Table 4 summarizes the number of LECs present at each order, excluding the
three $\pi N$ LECs at N3LO from the total number.

Table 3: Potential contributions at each chiral order in the channels where OPE is treated non-perturbatively. This table complements the information in Table 1. Order | Pion contribution | Contact terms
---|---|---
LO | $V^{(0)}_{1\pi}$ | $V^{(0)}_{\mathrm{ct}}:$
|  |  $C^{(0)}_{{}^{1}S_{0}}$, $\begin{pmatrix}C^{(0)}_{{}^{3}S_{1}}&0\\\ 0&0\end{pmatrix}$, $D^{(0)}_{{}^{3}P_{0}}p^{\prime}p$, $\begin{pmatrix}D^{(0)}_{{}^{3}P_{2}}p^{\prime}p&0\\\ 0&0\end{pmatrix}$
NLO | - |  $V^{(1)}_{\mathrm{ct}}$:
|  |  $D^{(0)}_{{}^{1}S_{0}}(p^{\prime 2}+p^{2})$, $C^{(1)}_{{}^{1}S_{0}}$
N2LO | $V_{2\pi}^{(2)}$ |  $V^{(2)}_{\mathrm{ct}}$:
|  |  $E^{(0)}_{{}^{1}S_{0}}p^{\prime 2}p^{2}$, $D^{(1)}_{{}^{1}S_{0}}(p^{\prime 2}+p^{2})$, $C^{(2)}_{{}^{1}S_{0}}$,
|  |  $\begin{pmatrix}D^{(0)}_{{}^{3}S_{1}}(p^{\prime 2}+p^{2})&D^{(0)}_{SD}p^{2}\\\ D^{(0)}_{SD}p^{\prime 2}&0\end{pmatrix}$, $\begin{pmatrix}C^{(1)}_{{}^{3}S_{1}}&0\\\ 0&0\end{pmatrix}$,
|  |  $E^{(0)}_{{}^{3}P_{0}}p^{\prime}p(p^{\prime 2}+p^{2})$ , $D^{(1)}_{{}^{3}P_{0}}p^{\prime}p$,
|  |  $p^{\prime}p\begin{pmatrix}E^{(0)}_{{}^{3}P_{2}}(p^{\prime 2}+p^{2})&E^{(0)}_{PF}p^{2}\\\ E^{(0)}_{PF}p^{\prime 2}&0\end{pmatrix}$, $\begin{pmatrix}D^{(1)}_{{}^{3}P_{2}}p^{\prime}p&0\\\ 0&0\end{pmatrix}$,
|  |  $D^{(0)}_{{}^{1}P_{1}}p^{\prime}p$, $D^{(0)}_{{}^{3}P_{1}}p^{\prime}p$
N3LO |  $V_{2\pi}^{(3)}$, (includes |  $V^{(3)}_{\mathrm{ct}}$:
|  $\pi N$ LECs: $c_{1},c_{3},c_{4}$) |  $F^{(0)}_{{}^{1}S_{0}}p^{\prime 2}p^{2}(p^{\prime 2}+p^{2})$, $E^{(1)}_{{}^{1}S_{0}}p^{\prime 2}p^{2}$, $D^{(2)}_{{}^{1}S_{0}}(p^{\prime 2}+p^{2})$, $C^{(3)}_{{}^{1}S_{0}}$,
|  |  $\begin{pmatrix}D^{(1)}_{{}^{3}S_{1}}(p^{\prime 2}+p^{2})&D^{(1)}_{SD}p^{2}\\\ D^{(1)}_{SD}p^{\prime 2}&0\end{pmatrix}$, $\begin{pmatrix}C^{(2)}_{{}^{3}S_{1}}&0\\\ 0&0\end{pmatrix}$,
|  |  $E^{(1)}_{{}^{3}P_{0}}p^{\prime}p(p^{\prime 2}+p^{2})$, $D^{(2)}_{{}^{3}P_{0}}p^{\prime}p$,
|  |  $p^{\prime}p\begin{pmatrix}E^{(1)}_{{}^{3}P_{2}}(p^{\prime 2}+p^{2})&E^{(1)}_{PF}p^{2}\\\ E^{(1)}_{PF}p^{\prime 2}&0\end{pmatrix}$, $\begin{pmatrix}D^{(2)}_{{}^{3}P_{2}}p^{\prime}p&0\\\ 0&0\end{pmatrix}$,
|  |  $D^{(1)}_{{}^{1}P_{1}}p^{\prime}p$, $D^{(1)}_{{}^{3}P_{1}}p^{\prime}p$
Table 4: The number of LECs at each order in the Long and Yang PC. Chiral order | New LECs | Pert. correction | Total up to order
---|---|---|---
LO | 4 | – | 4
NLO | 1 | 1 | 6
N2LO | 8 | 5 | 19
N3LO | 1 (+3222Sub-leading $\pi N$ LECs: $c_{1},c_{3},c_{4}$ excluded from the total in the last column.) | 13 | 33

##  Appendix B Numerical implementation of distorted-wave perturbation theory

This appendix gives some more details regarding the implementation of the
equations for higher-order corrections to the scattering amplitude in Eqs. 11,
12 and 13. Since all operator products reduce to the form in Eq. 19, the
implementation can be done in complete analogy with the solution of the
partial-wave Lippmann-Schwinger equation using Gauss-Legendre quadrature [73,
74].

In this appendix we suppress the conserved quantum numbers $s$ and $j$, and
write the resolution of identity in the partial wave basis as

$\mathds{1}=\sum_{\ell}\int_{0}^{\infty}dk\ k^{2}\ket{k,\ell}\bra{k,\ell}.$
(43)

Furthermore, for a stationary proton (mass $m_{p}$) and an incoming neutron
(mass $m_{n}$) with kinetic energy $T_{\mathrm{lab}}$ in the laboratory frame
of reference, the modulus of the c.m.momentum, $p_{0}$, is given by

$p_{0}^{2}=\frac{m_{p}^{2}T_{\mathrm{lab}}(2m_{n}+T_{\mathrm{lab}})}{(m_{n}+m_{p})^{2}+2m_{p}T_{\mathrm{lab}}}.$
(44)

By inserting the resolution of identity in Eq. 19 and discretizing the
integral using Gauss-Legendre quadrature with momentum points and weights,
$\\{k_{i},w_{i}\\}_{i=1}^{N}$, we obtain

$\displaystyle\braket{p^{\prime},\ell^{\prime}}{A_{1}G^{+}_{0}A_{2}}{p,\ell}$
$\displaystyle=\sum_{\ell^{\prime\prime},\ell^{\prime\prime\prime}}\int_{0}^{\infty}dk_{1}\
k_{1}^{2}\int_{0}^{\infty}dk_{2}\
k^{2}_{2}\braket{p^{\prime},\ell^{\prime}}{A_{1}}{k_{1},\ell^{\prime\prime}}\braket{k_{1},\ell^{\prime\prime}}{G^{+}_{0}}{k_{2},\ell^{\prime\prime\prime}}\braket{k_{2},\ell^{\prime\prime\prime}}{A_{2}}{p,\ell}=$
$\displaystyle=\sum_{\ell^{\prime\prime}}\int_{0}^{\infty}dk_{1}\
k_{1}^{2}\braket{p,\ell^{\prime}}{A_{1}}{k_{1},\ell^{\prime\prime}}\frac{m_{N}}{p_{0}^{2}-k_{1}^{2}+i\epsilon}\braket{k_{1},\ell^{\prime\prime}}{A_{2}}{p,\ell}=$
(45)
$\displaystyle=\sum_{\ell^{\prime\prime}}\sum_{i=1}^{N}k^{2}_{i}w_{i}\braket{p,\ell^{\prime}}{A_{1}}{k_{i},\ell^{\prime\prime}}\frac{m_{N}}{p_{0}^{2}-k_{i}^{2}+i\epsilon}\braket{k_{i},\ell^{\prime\prime}}{A_{2}}{p,\ell}.$
(46)

Here, $p_{0}$ denotes the on-shell momentum for a given scattering energy
$T_{\mathrm{lab}}$ given by Eq. 44. Doing some manipulations and converting
the $+i\epsilon$ prescription to a principal value we obtain [65, 74]

$\displaystyle\braket{p^{\prime},\ell^{\prime}}{A_{1}G^{+}_{0}A_{2}}{p,\ell}$
$\displaystyle=\sum_{l^{\prime\prime}}\sum_{i=1}^{N}k_{i}^{2}w_{i}\braket{p^{\prime},\ell^{\prime}}{A_{1}}{k_{i},\ell^{\prime\prime}}\frac{m_{N}}{p_{0}^{2}-k_{i}^{2}}\braket{k_{i},\ell^{\prime\prime}}{A_{2}}{p,\ell}$
$\displaystyle-\braket{p^{\prime},\ell^{\prime}}{A_{1}}{p_{0},\ell^{\prime\prime}}\braket{p_{0},\ell^{\prime\prime}}{A_{2}}{p,\ell}\left[m_{N}p_{0}^{2}\sum_{i=1}^{N}\frac{w_{i}}{p_{0}^{2}-k_{i}^{2}}+\frac{i\pi
m_{N}p_{0}}{2}-m_{N}p_{0}\
\mathrm{arctanh}\left(\frac{p_{0}}{\tilde{\Lambda}}\right)\right].$ (47)

All potentials are regulated using Eq. 16 and at sufficiently high momentum,
$\tilde{\Lambda}$, all potential matrix elements are essentially zero. This
means that the integral in Eq. 46 is well represented by the discretized sum
where the momentum points and weights $\\{k_{i},w_{i}\\}_{i=1}^{N}$ are chosen
using Gauss-Legendre quadrature in the interval $[0,\tilde{\Lambda}]$. The
last term in the bracket in Eq. 47 implements the principal-value integral on
the interval $[\tilde{\Lambda},\infty]$ analytically since the grid is just
doing the integration on $[0,\tilde{\Lambda}]$ [75]. It is possible to have a
grid that extends to numerical infinity, but this generally leads to slower
convergence with $N$. For the calculations in this study, we employ
$\tilde{\Lambda}=\Lambda+1500$ MeV, for both $\Lambda=500$ MeV and
$\Lambda=2500$ MeV, which we find sufficient for numerical convergence.

Equation 47 can be expressed in a simpler form using matrix products, which
speeds up the computations. We define the propagator matrix as

$[G^{+}_{0}]_{ij}=\delta_{ij}F_{i},\quad
F_{i}=\begin{cases}\frac{m_{N}}{p_{0}^{2}-k_{i}^{2}},\quad i=1,...,N\\\
-f(p_{0}),\quad i=N+1,\end{cases}$ (48)

where

$f(p_{0})=m_{N}p_{0}^{2}\sum_{i=1}^{N}\frac{w_{i}}{p_{0}^{2}-k_{i}^{2}}+\frac{i\pi
m_{N}p_{0}}{2}-m_{N}p_{0}\
\mathrm{arctanh}\left(\frac{p_{0}}{\tilde{\Lambda}}\right).$ (49)

Similarly, we make the following definitions of matrices for $A_{\mu}$,
$\mu=1,2$,

$\displaystyle[A_{\mu}^{\ell^{\prime}\ell}]_{i,j}$
$\displaystyle=k_{i}\sqrt{w_{i}}\braket{k_{i},\ell^{\prime}}{A_{\mu}}{k_{j},\ell}k_{j}\sqrt{w_{j}},\quad
i,j=1,\dots,N$ (50) $\displaystyle[A_{\mu}^{\ell^{\prime}\ell}]_{i,j=N+1}$
$\displaystyle=k_{i}\sqrt{w_{i}}\braket{k_{i},\ell^{\prime}}{A_{\mu}}{p_{0},\ell},\quad
i=1,\dots,N$ (51) $\displaystyle[A_{\mu}^{\ell^{\prime}\ell}]_{i=N+1,j}$
$\displaystyle=\braket{p_{0},\ell^{\prime}}{A_{\mu}}{k_{j},\ell}k_{j}\sqrt{w_{j}},\quad
j=1,\dots,N$ (52) $\displaystyle[A_{\mu}^{\ell^{\prime}\ell}]_{i=N+1,j=N+1}$
$\displaystyle=\braket{p_{0},\ell^{\prime}}{A_{\mu}}{p_{0},\ell},$ (53)

effectively including an extra momentum-grid point $k_{N+1}\equiv p_{0}$ with
weight $\sqrt{w_{N+1}}=p_{0}^{-1}$. Using these definitions and defining
$D=A_{1}G^{+}_{0}A_{2}$, Eq. 47 can be written using $(N+1)\times(N+1)$ matrix
products

$[D^{\ell^{\prime}\ell}]_{ij}=\sum_{\ell^{\prime\prime}}\sum_{n,m=1}^{N+1}[A_{1}^{\ell^{\prime}\ell^{\prime\prime}}]_{in}[G^{+}_{0}]_{nm}[A_{2}^{\ell^{\prime\prime}\ell}]_{mj},\quad
i,j=1,\dots,N+1.$ (54)

For coupled channels, we further eliminate the sum over $\ell^{\prime\prime}$
in Eq. 54 by defining $(2N+2)\times(2N+2)$ block-matrices, which for $A_{1}$
reads

$[\bm{A}_{1}]=\begin{pmatrix}[A_{1}^{--}]&[A_{1}^{-+}]\\\
[A_{1}^{+-}]&[A_{1}^{++}]\end{pmatrix}.$ (55)

The $\pm$ notation represents $\ell=j\pm 1$. The propagator is diagonal in
$\ell$ and can be written as

$[\bm{G}^{+}_{0}]=\begin{pmatrix}[\bm{G}^{+}_{0}]&0\\\
0&[\bm{G}^{+}_{0}]\end{pmatrix}.$ (56)

We can finally write Eq. 54 as

$[\bm{D}]=[\bm{A}_{1}][\bm{G}^{+}_{0}][\bm{A}_{2}].$ (57)

Note that the simplification of Eq. 47 to an ordinary matrix product in Eq. 57
is only possible due to the specific structure of having $G^{+}_{0}$ in
between $A_{1}$ and $A_{2}$. This structure gives rise to the last “on-shell”
term in (47) that can be incorporated by adding the grid point
$k_{N+1}=p_{0}$, which then extends the sum in Eq. 47 to $N+1$. Equation 54
can now be used recursively to compute longer products such as
$\braket{p^{\prime},\ell^{\prime}}{A_{1}G^{+}_{0}A_{2}G^{+}_{0}A_{3}}{p,\ell}$.

As an example, the first-order correction to the $T$-matrix in Eq. 11 can be
expressed as the matrix equation

$[\bm{T}^{(1)}]=\left(\mathds{1}+[\bm{T}^{(0)}][\bm{G}^{+}_{0}]\right)[\bm{V}^{(1)}]\left(\mathds{1}+[\bm{G}^{+}_{0}][\bm{T}^{(0)}]\right).$
(58)

##  Appendix C Perturbative phase shifts

In this appendix we discuss how to obtain phase shifts given perturbative
corrections to the $T$-matrix computed from Eqs. 11, 12 and 13. We will follow
the method outlined in Ref. [32] and add some additional details. For
uncoupled scattering channels, the $1\times 1$ $S$-matrix can be parameterized
by

$S=\exp\left(2i\delta\right),$ (59)

where $\delta$ is the phase shift. We expand both the phase shifts and the on-
shell $S$-matrix with the contributions at each chiral order obtaining

$\displaystyle S^{(0)}+S^{(1)}+S^{(2)}+S^{(3)}+\mathcal{O}(Q^{3})=$ (60)
$\displaystyle\exp\left(2i\left[\delta^{(0)}+\delta^{(1)}+\delta^{(2)}+\delta^{(3)}+\mathcal{O}(Q^{3})\right]\right).$
(61)

Performing a Taylor expansion of both sides, and matching chiral orders, gives

$\displaystyle S^{(0)}$ $\displaystyle=\exp\left(2i\delta^{(0)}\right)$ (62)
$\displaystyle S^{(1)}$
$\displaystyle=2i\delta^{(1)}\exp\left(2i\delta^{(0)}\right)$ (63)
$\displaystyle S^{(2)}$
$\displaystyle=\left[2i\delta^{(2)}-2\left(\delta^{(1)}\right)^{2}\right]\exp\left(2i\delta^{(0)}\right)$
(64) $\displaystyle S^{(3)}$
$\displaystyle=\left[2i\delta^{(3)}-4\delta^{(1)}\delta^{(2)}-\frac{4i}{3}\left(\delta^{(1)}\right)^{3}\right]\exp\left(2i\delta^{(0)}\right)$
(65)

From these equations, we straightforwardly obtain explicit expressions for the
LO phase shift $\delta^{(0)}$ (trivial), and all corrections
$\\{\delta^{(\nu)}\\}_{\nu>0}$. We note that all corrections are real valued.
To obtain the total phase shift at, e.g., N2LO, one has to sum
$\delta^{(0)}+\delta^{(1)}+\delta^{(2)}$. The $S$-matrix corrections are
obtained from the $T$-matrix corrections as

$S^{(\nu)}_{\ell^{\prime}\ell}=-i\pi
m_{N}p_{0}T^{(\nu)}_{\ell^{\prime}\ell},\quad\nu>0,$ (66)

for a given on-shell momentum, $p_{0}$.

For coupled channels we use the Stapp-parametrization [76] for the on-shell
$2\times 2$ $S$-matrix

$S=\begin{pmatrix}\cos(2\epsilon)e^{2i\delta_{1}}&i\sin(2\epsilon)e^{i(\delta_{1}+\delta_{2})}\\\
i\sin(2\epsilon)e^{i(\delta_{1}+\delta_{2})}&\cos(2\epsilon)e^{2i\delta_{2}}\end{pmatrix},$
(67)

where the three phase shifts $\delta_{1}$, $\delta_{2}$ and $\epsilon$
parameterize the amplitude for a given channel. We now proceed completely
analogous to the uncoupled case, dividing the $S$-matrix and phase shifts into
chiral orders as

$S=\sum_{\nu=0}^{\infty}=S^{(\nu)},\quad\delta_{1}=\sum_{\nu=0}^{\infty}\delta^{(\nu)}_{1},\quad\delta_{2}=\sum_{\nu=0}^{\infty}\delta^{(\nu)}_{2},\quad\epsilon=\sum_{\nu=0}^{\infty}\epsilon^{(\nu)}.$
(68)

For convenience, we define the functions

$\displaystyle f_{11}(\epsilon,\delta_{1})$
$\displaystyle=\cos(2\epsilon)e^{2i\delta_{1}},$ (69) $\displaystyle
f_{12}(\epsilon,\delta_{1},\delta_{2})$
$\displaystyle=i\sin(2\epsilon)e^{i(\delta_{1}+\delta_{2})},$ (70)
$\displaystyle f_{22}(\epsilon,\delta_{2})$
$\displaystyle=\cos(2\epsilon)e^{2i\delta_{2}},$ (71)

which are the constituents of the matrix in Eq. 67. Inserting the expansions
in Eq. 68 into Eq. 67, Taylor expanding and matching chiral orders, gives the
perturbative corrections to the phase shifts. Expanding the upper left matrix
element of $S$ gives

$\displaystyle S_{11}^{(0)}$ $\displaystyle=f_{11}$ (72) $\displaystyle
S_{11}^{(1)}$
$\displaystyle=\partial_{\epsilon}f_{11}\times\epsilon^{(1)}+\partial_{\delta}f_{11}\times\delta^{(1)}$
(73) $\displaystyle S_{11}^{(2)}$
$\displaystyle=\partial_{\epsilon}f_{11}\times\epsilon^{(2)}+\partial_{\delta}f_{11}\times\delta^{(2)}$
$\displaystyle+g^{(2)}_{11}(\epsilon^{(1)},\delta^{(1)})$ (74) $\displaystyle
S_{11}^{(3)}$
$\displaystyle=\partial_{\epsilon}f_{11}\times\epsilon^{(3)}+\partial_{\delta}f_{11}\times\delta^{(3)}$
$\displaystyle+g^{(3)}_{11}(\epsilon^{(1)},\delta^{(1)},\epsilon^{(2)},\delta^{(2)})$
(75)

where the functions $g^{(\nu)}_{11}$ are introduced to capture all non-linear
terms in the expansion

$\displaystyle g^{(2)}_{11}(\epsilon^{(1)},\delta^{(1)})$
$\displaystyle=\frac{1}{2}\partial^{2}_{\epsilon}f_{11}\times\left(\epsilon^{(1)}\right)^{2}+\frac{1}{2}\partial^{2}_{\delta}f_{11}\times\left(\delta^{(1)}\right)^{2}+\partial_{\epsilon}\partial_{\delta}f_{11}\times\delta^{(1)}\epsilon^{(1)}$
(76) $\displaystyle
g^{(3)}_{11}(\epsilon^{(1)},\delta^{(1)},\epsilon^{(2)},\delta^{(2)})$
$\displaystyle=\partial^{2}_{\epsilon}f_{11}\left(\epsilon^{(1)}\epsilon^{(2)}\right)+\partial_{\epsilon}\partial_{\delta}f_{11}\left(\epsilon^{(1)}\delta^{(2)}+\epsilon^{(2)}\delta^{(1)}\right)$
$\displaystyle+\partial^{2}_{\delta}f_{11}\left(\delta^{(1)}\delta^{(2)}\right)+\frac{1}{6}\partial^{3}_{\epsilon}f_{11}\left(\epsilon^{(1)}\right)^{3}+\frac{1}{2}\partial_{\delta}\partial^{2}_{\epsilon}f_{11}\left(\epsilon^{(1)}\right)^{2}\delta^{(1)}$
$\displaystyle+\frac{1}{2}\partial^{2}_{\delta}\partial_{\epsilon}f_{11}\epsilon^{(1)}\left(\delta^{(1)}\right)^{2}+\frac{1}{6}\partial^{3}_{\delta}f_{11}\left(\delta^{(1)}\right)^{3}.$
(77)

Since $f_{11}$ depends on $\epsilon$ and $\delta_{1}$ the index one is
suppressed. The function $f_{11}$ and all its derivatives are evaluated at
$(\epsilon^{(0)},\delta^{(0)}_{1})$.

For the lower right matrix element described by $f_{22}$ the expressions are
completely analogous to Eqs. 76 and 77, but with $\delta_{2}$ instead of
$\delta_{1}$. For the off-diagonal elements we get

$\displaystyle S_{12}^{(0)}$ $\displaystyle=f_{12}$ (78) $\displaystyle
S_{12}^{(1)}$
$\displaystyle=\partial_{\epsilon}f_{12}\times\epsilon^{(1)}+\partial_{\delta_{1}}f_{12}\times\delta^{(1)}_{1}+\partial_{\delta_{2}}f_{12}\times\delta^{(1)}_{2}$
(79) $\displaystyle S_{12}^{(2)}$
$\displaystyle=\partial_{\epsilon}f_{12}\times\epsilon^{(2)}+\partial_{\delta_{1}}f_{12}\times\delta^{(2)}_{1}+\partial_{\delta_{2}}f_{12}\times\delta^{(2)}_{2}+g^{(2)}_{12}(\epsilon^{(1)},\delta^{(1)}_{1},\delta^{(1)}_{2})$
(80) $\displaystyle S_{12}^{(3)}$
$\displaystyle=\partial_{\epsilon}f_{12}\times\epsilon^{(3)}+\partial_{\delta_{1}}f_{12}\times\delta^{(3)}_{1}+\partial_{\delta_{2}}f_{12}\times\delta^{(3)}_{2}+g^{(3)}_{12}(\epsilon^{(1)},\delta^{(1)}_{1},\delta^{(1)}_{2},\epsilon^{(2)},\delta^{(2)}_{1},\delta^{(2)}_{2}),$
(81)

where the functions $g^{(\nu)}_{12}$ capture the non-linear terms

$\displaystyle g^{(2)}_{12}(\epsilon^{(1)},\delta^{(1)}_{1},\delta^{(1)}_{2})$
$\displaystyle=\frac{1}{2}\partial^{2}_{\epsilon}f_{12}\times\left(\epsilon^{(1)}\right)^{2}+\frac{1}{2}\partial^{2}_{\delta_{1}}f_{12}\times\left(\delta^{(1)}_{1}\right)^{2}+\frac{1}{2}\partial^{2}_{\delta_{2}}f_{12}\times\left(\delta^{(1)}_{2}\right)^{2}$
$\displaystyle+\partial_{\epsilon}\partial_{\delta_{1}}f_{12}\epsilon^{(1)}\delta^{(1)}_{1}+\partial_{\epsilon}\partial_{\delta_{2}}f_{12}\epsilon^{(1)}\delta^{(1)}_{2}+\partial_{\delta_{1}}\partial_{\delta_{2}}f_{12}\delta^{(1)}_{1}\delta^{(1)}_{2}$
(82) $\displaystyle
g^{(3)}_{12}(\epsilon^{(1)},\delta^{(1)}_{1},\delta^{(1)}_{2},\epsilon^{(2)},\delta^{(2)}_{1},\delta^{(2)}_{2})$
$\displaystyle=\partial^{2}_{\epsilon}f_{12}\epsilon^{(1)}\epsilon^{(2)}+\partial^{2}_{\delta_{1}}f_{12}\delta^{(1)}_{1}\delta^{(2)}_{1}+\partial^{2}_{\delta_{2}}f_{12}\delta^{(1)}_{2}\delta^{(2)}_{2}$
$\displaystyle+\partial_{\epsilon}\partial_{\delta_{1}}f_{12}\left(\epsilon^{(1)}\delta^{(2)}_{1}+\epsilon^{(2)}\delta^{(1)}_{1}\right)+\partial_{\epsilon}\partial_{\delta_{2}}f_{12}\left(\epsilon^{(1)}\delta^{(2)}_{2}+\epsilon^{(2)}\delta^{(1)}_{2}\right)$
$\displaystyle+\partial_{\delta_{1}}\partial_{\delta_{2}}f_{12}\left(\delta^{(1)}_{1}\delta^{(2)}_{2}+\delta^{(2)}_{1}\delta^{(1)}_{2}\right)$
$\displaystyle+\frac{1}{2}\partial^{2}_{\epsilon}\partial_{\delta_{1}}f_{12}\left(\epsilon^{(1)}\right)^{2}\delta^{(1)}_{1}+\frac{1}{2}\partial^{2}_{\epsilon}\partial_{\delta_{2}}f_{12}\left(\epsilon^{(1)}\right)^{2}\delta^{(1)}_{2}$
$\displaystyle+\frac{1}{2}\partial_{\epsilon}\partial^{2}_{\delta_{1}}f_{12}\epsilon^{(1)}\left(\delta^{(1)}_{1}\right)^{2}+\frac{1}{2}\partial_{\delta_{2}}\partial^{2}_{\delta_{1}}f_{12}\delta^{(1)}_{2}\left(\delta^{(1)}_{1}\right)^{2}$
$\displaystyle+\frac{1}{2}\partial_{\epsilon}\partial^{2}_{\delta_{2}}f_{12}\epsilon^{(1)}\left(\delta^{(1)}_{2}\right)^{2}+\frac{1}{2}\partial_{\delta_{1}}\partial^{2}_{\delta_{2}}f_{12}\delta^{(1)}_{1}\left(\delta^{(1)}_{2}\right)^{2}+$
$\displaystyle+\partial_{\epsilon}\partial_{\delta_{1}}\partial_{\delta_{2}}f_{12}\epsilon^{(1)}\delta^{(1)}_{1}\delta^{(1)}_{2}$
$\displaystyle+\frac{1}{6}\partial^{3}_{\epsilon}f_{12}\left(\epsilon^{(1)}\right)^{3}+\frac{1}{6}\partial^{3}_{\delta_{1}}f_{12}\left(\delta^{(1)}_{1}\right)^{3}+\frac{1}{6}\partial^{3}_{\delta_{2}}f_{12}\left(\delta^{(1)}_{2}\right)^{3}.$
(83)

The function $f_{12}$ and all its derivatives are evaluated at
$(\epsilon^{(0)},\delta^{(0)}_{1},\delta^{(0)}_{2})$. Note that all the
functions $g^{(\nu)}_{**}$ vanish if the NLO corrections
$(\delta^{(1)}_{1},\delta^{(1)}_{2},\epsilon^{(1)})$ are zero. This is the
case for all coupled channels where OPE is treated non-perturbatively as seen
in Table 3. Furthermore, in all channels where OPE is treated perturbatively
the LO phase shifts are all zero, which makes many of the terms in the
expressions for $g^{(\nu)}_{**}$ vanish due to vanishing derivatives. Thus, in
both the perturbative and non-perturbative cases, Eqs. 76, 77, 82 and 83 can
be simplified substantially. The phase shift corrections
$(\epsilon^{(\nu)},\delta^{(\nu)}_{1},\delta^{(\nu)}_{2})$ for $\nu=1,2,3.$
are finally obtained by solving a system of linear equations

$\displaystyle\mathrm{\textbf{NLO}}:$
$\displaystyle\quad\begin{pmatrix}S_{11}^{(1)}\\\ S_{12}^{(1)}\\\
S_{22}^{(1)}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\partial_{\epsilon}f_{11}&\partial_{\delta_{1}}f_{11}&0\\\
\partial_{\epsilon}f_{12}&\partial_{\delta_{1}}f_{12}&\partial_{\delta_{2}}f_{12}\\\
\partial_{\epsilon}f_{22}&0&\partial_{\delta_{2}}f_{22}\end{pmatrix}\begin{pmatrix}\epsilon^{(1)}\\\
\delta^{(1)}_{1}\\\ \delta^{(1)}_{2}\end{pmatrix}$ (84)
$\displaystyle\mathrm{\textbf{N${}^{2}$LO}}:$
$\displaystyle\quad\begin{pmatrix}S_{11}^{(2)}-g^{(2)}_{11}\\\
S_{12}^{(2)}-g^{(2)}_{12}\\\ S_{22}^{(2)}-g^{(2)}_{22}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\partial_{\epsilon}f_{11}&\partial_{\delta_{1}}f_{11}&0\\\
\partial_{\epsilon}f_{12}&\partial_{\delta_{1}}f_{12}&\partial_{\delta_{2}}f_{12}\\\
\partial_{\epsilon}f_{22}&0&\partial_{\delta_{2}}f_{22}\end{pmatrix}\begin{pmatrix}\epsilon^{(2)}\\\
\delta^{(2)}_{1}\\\ \delta^{(2)}_{2}\end{pmatrix}$ (85)
$\displaystyle\mathrm{\textbf{N${}^{3}$LO}}:$
$\displaystyle\quad\begin{pmatrix}S_{11}^{(3)}-g^{(3)}_{11}\\\
S_{12}^{(3)}-g^{(3)}_{12}\\\ S_{22}^{(2)}-g^{(3)}_{22}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\partial_{\epsilon}f_{11}&\partial_{\delta_{1}}f_{11}&0\\\
\partial_{\epsilon}f_{12}&\partial_{\delta_{1}}f_{12}&\partial_{\delta_{2}}f_{12}\\\
\partial_{\epsilon}f_{22}&0&\partial_{\delta_{2}}f_{22}\end{pmatrix}\begin{pmatrix}\epsilon^{(3)}\\\
\delta^{(3)}_{1}\\\ \delta^{(3)}_{2}\end{pmatrix}.$ (86)