$\psi^{+}_{{\mathbf{p}}}({\mathbf{r}})$ obtained by solving the Schrödinger
equation with LO isospin-symmetric chiral potential

$V_{{}^{1}S_{0}}({\mathbf{r}})=-\alpha_{\pi}\frac{e^{-M_{\pi}r}}{r}+C\,\delta^{(3)}({\mathbf{r}})\,,\qquad\alpha_{\pi}=\frac{g_{A}^{2}M_{\pi}^{2}}{4\pi
F_{\pi}^{2}}\,,\qquad M_{\pi}=\frac{M_{\pi^{0}}+2M_{\pi^{\pm}}}{3}\,,$ (7.5)

with $C$ tuned to reproduce either the $np$ or $nn$ scattering length (the
difference being negligible at the chosen kinematic point). To solve for the
wave functions we have used the methods described in Ref. [23], and to compute
the amplitude, isolating its singular component, we have used the approach
described in detail in Ref. [43].

Our amplitude satisfies Watson’s theorem [106] and at the kinematic point
$|{\mathbf{p}}|=25\,\text{MeV}$, $|{\mathbf{p}}^{\prime}|=30\,\text{MeV}$ we
find161616The amplitude $\mathcal{A}_{\nu}$ is related to the $S$-matrix
element for the process $n(\mathbf{p})\ n(-\mathbf{p})\to
p(\mathbf{p}^{\prime})\ p(-\mathbf{p}^{\prime})\ e(\mathbf{p}_{e})\
e(-\mathbf{p}_{e})$ by
$S_{\nu}=i(2\pi)^{4}\,\delta^{(4)}(p_{f}-p_{i})\,(4G_{F}^{2}V_{ud}^{2}m_{\beta\beta}\
\bar{u}_{L}(\mathbf{p}_{e})u_{L}^{c}(-\mathbf{p}_{e}))\,{\mathcal{A}}_{\nu}$ .

$\displaystyle\mathcal{A}_{\nu}(|{\mathbf{p}}|,|{\mathbf{p}}^{\prime}|)\times
e^{-i(\delta_{{}^{1}S_{0}}(|{\mathbf{p}}|)+\delta_{{}^{1}S_{0}}(|{\mathbf{p}}^{\prime}|))}$
$\displaystyle=-\left(2.271-0.075\ \tilde{\cal C}_{1}(4\,M_{\pi})\right)\times
10^{-2}\,\text{MeV}^{-2}$ $\displaystyle=-1.95(5)_{\tilde{\cal C}_{1}}\times
10^{-2}\,\text{MeV}^{-2}\,,$ (7.6)

where in the second line we have used our prediction $\tilde{\cal
C}_{1}(4M_{\pi})=4.2(6)$. At the chosen kinematic point and chiral
renormalization scale we find that the contact-term contribution interferes
destructively with the long-range neutrino exchange and reduces the amplitude
by about $15\%$. Since, as discussed in Refs. [40, 43], the effect of a
contact term of natural size that affects $\Delta I=0$ transitions such as
$nn\to pp$ at the $(10\text{--}20)\%$ level, becomes amplified to the level of
$(50\text{--}70)\%$ in $\Delta I=2$ nuclear transitions due to a node in the
matrix element density, it is possible that a more pronounced effect occurs in
nuclei of interest for $0\nu\beta\beta$ searches, but we stress that our
observation is based on the $\overline{\text{MS}}$ scheme at $\mu=4M_{\pi}$
and thus need not carry over to other scales and schemes.

Figure 14: Left: values of $\tilde{\cal C}_{1}(\Lambda_{p},n)$ as function of
$\Lambda_{p}$ for 3 choices of the regulator $n=2$ (solid), $n=3$ (dashed),
and $n=4$ (dotted). Right: values of $\tilde{\cal C}_{1}(R_{S})$ as function
of the coordinate-space regulator $R_{S}$.

We advocate using the synthetic data (7.6) to determine the LEC $\tilde{\cal
C}_{1}$ in any regularization/renormalization scheme applicable in many-body
nuclear calculations [26, 27, 107, 28, 108]. As an example, we show here how
to perform such a determination for often-applied momentum- and coordinate-
space regulators. In momentum space, the strong potential is regulated through
an exponential cutoff of the form

$V_{\mathrm{strong}}(\mathbf{p}^{\prime},\mathbf{p})\rightarrow\mathrm{exp}\left[-\left(\frac{\mathbf{p}^{\prime\,2}}{\Lambda_{p}^{2}}\right)^{n}\right]\,V_{\mathrm{strong}}(\mathbf{p}^{\prime},\mathbf{p})\,\mathrm{exp}\left[-\left(\frac{\mathbf{p}^{2}}{\Lambda_{p}^{2}}\right)^{n}\right]\,,$
(7.7)

in terms of a momentum cutoff $\Lambda_{p}$. The LS equation is solved
numerically for different values of $\Lambda_{p}$ and $n$ and the strong
counter term, $C(\Lambda_{p},n)$ in the ${}^{1}S_{0}$ channel, is determined
from a fit to the $N\\!N$ scattering lengths. We then insert the long- and
short-range neutrino potential between initial- and final-state scattering
states and determine $\tilde{C}_{1}(\Lambda_{p},n)$ by fitting to Eq. (7.6).
An analogous procedure is followed in coordinate space. In this case, contact
interactions are regulated by a local Gaussian regulator

$\delta^{(3)}({\mathbf{r}})\rightarrow\frac{1}{\left(\sqrt{\pi}R_{S}\right)^{3}}\exp\left(-\frac{r^{2}}{R_{S}^{2}}\right).$
(7.8)

The $pp$ and $nn$ wave functions obtained for a fixed value of $R_{S}$ are
then used to compute the amplitude (7.3).

Figure 15: Top left: Long-distance (dashed), short-distance (dotted), and
total (solid) contributions to
$\mathcal{A}_{\nu}(|{\mathbf{p}}|=25\,\text{MeV},|{\mathbf{p}}^{\prime}|)$ for
$\Lambda_{p}=2\,\text{fm}^{-1}$ (blue) and $\Lambda_{p}=20\,\text{fm}^{-1}$
(red) (the total contribution is shown in black, as it is identical for the
two cutoffs). Top right: ratio of short-to-long-distance contributions for
$\Lambda_{p}=2\,\text{fm}^{-1}$ (blue), $\Lambda_{p}=5\,\text{fm}^{-1}$
(green, dotted), and $\Lambda_{p}=20\,\text{fm}^{-1}$ (red, dashed). Bottom:
same as top, but with the coordinate space regulator $R_{S}$, with
$R^{-1}_{S}=1.25\,\text{fm}^{-1}$ (blue), $R^{-1}_{S}=2\,\text{fm}^{-1}$
(green), and $R^{-1}_{S}=20\,\text{fm}^{-1}$ (red).

In Fig. 14 we depict the values of the counter term $\tilde{\cal C}_{1}$ in
the two schemes. In the left panel, we show $\tilde{\cal
C}_{1}(\Lambda_{p},n)$ for a wide range of $\Lambda_{p}$ and three choices of
$n$, showing a clear logarithmic dependence of $\Lambda_{p}$ and a mild
dependence on $n$. In the right panel, we show the same plot in the $R_{S}$
scheme, which again exhibits a logarithmic dependence on $R_{S}$, with power
corrections becoming important at small $R^{-1}_{S}$.

The relative importance of the short-distance contribution depends on the
regulator and on the kinematics. In Fig. 15 we show the amplitude
$\mathcal{A}_{\nu}$ for $|\mathbf{p}|=25$ MeV, while varying
$|\mathbf{p}^{\prime}|$ between $26$ and $46$ MeV. The region includes the
fitting point $|\mathbf{p}^{\prime}|=30$ MeV, and at these small momenta the
LO chiral wave functions provide a good description of the ${}^{1}S_{0}$ phase
shift. The plots in the top panel use the momentum-space regulator
$\Lambda_{p}$, while the bottom plots use the coordinate-space $R_{S}$ scheme.
In the top-left panel, we have chosen two specific regulators
$\Lambda_{p}=2\,\text{fm}^{-1}$ (blue) and $\Lambda_{p}=20\,\text{fm}^{-1}$
(red) and kept $n=2$ as the dependence on $n$ is mild. We show individually
the long-distance (dashed), short-distance (dotted), and sum (solid)
contributions. The sum does not depend on $\Lambda_{p}$, but the counter term
goes from constructive ($+15\%$) to destructive ($-20\%$) between the two
choices of the regulators. This is shown more clearly in the right panel of
Fig. 15, where the ratio of short-to-long-distance contributions is given for
three choices of regulators $\Lambda_{p}=\\{2,\,5,\,20\\}\,\,\text{fm}^{-1}$
in blue, green, and red respectively. This confirms explicitly that the
separation of long- and short-distance contributions is not physical.

The bottom panel shows the same plots in the $R_{S}$ scheme, using
$R^{-1}_{S}=1.25\,\text{fm}^{-1}$ (green), $R^{-1}_{S}=2\,\text{fm}^{-1}$
(blue), and $R^{-1}_{S}=20\,\text{fm}^{-1}$ (red). Also in this case, the
interference goes from destructive at large cutoff to constructive at small
cutoff, with the short-range amplitude being approximately zero at
$2\,\text{fm}^{-1}$. Of course, the total amplitude agrees between the
$\Lambda_{p}$ and $R_{S}$ regularization schemes as it should. These plots
show that, as already emphasized above, we cannot say whether the new short-
distance contribution will add constructively or destructively in ab-initio
calculations of nuclear transitions. This question can only be answered within
a specific regularization scheme and choice for the strong $N\\!N$ potential.

##  8 Conclusions

In this work we have presented the details of the new method developed in Ref.
[59] to estimate the LO contact-term contribution to the amplitude $nn\to
ppe^{-}e^{-}$ through the exchange of light Majorana neutrinos. This is
currently the missing ingredient in order to construct complete LO $nn\to pp$
transition operators in chiral EFT, to be used in ab-initio nuclear-structure
calculations of matrix elements in nuclei of experimental interest for
$0\nu\beta\beta$ searches.

Our approach to estimate the contact term is based on a representation of the
$nn\to ppe^{-}e^{-}$ amplitude as the momentum integral of the neutrino
propagator ($1/k^{2}$) times the generalized forward Compton scattering
amplitude $n(p_{1})n(p_{2})W^{+}(k)\to\
p(p_{1}^{\prime})p(p_{2}^{\prime})W^{-}(k)$, in close analogy to the
Cottingham formula for the electromagnetic mass splittings of pions and
nucleons. To extract the contact term from this integral representation, we
have constructed model-independent descriptions of the integrand in the low-
momentum region, using chiral EFT, and in the high-momentum region, using the
OPE. In the low- and intermediate-momentum region we have kept only the
elastic contribution, i.e., the effect of two-nucleon intermediate states,
with the most important momentum dependence generated by single-nucleon form
factors as well as the $N\\!N$ amplitude itself. While we do not have a strict
dispersive derivation as in the case of hadron mass splittings, this approach
has the potential to match the $(20\text{--}30)\%$ accuracy of the elastic
approximation observed there, an expectation that we verified by studying CIB
in low-energy $N\\!N$ scattering within the same framework and reproducing,
within uncertainties, the CIB contribution to the ${}^{1}S_{0}$ $N\\!N$
scattering lengths. This phenomenological success gives us confidence that our
method is sound and the uncertainty estimate realistic.

The phenomenological validation is particularly important given that our
method does introduce model-dependent input in the intermediate-momentum
region, albeit anchored to known constraints from QCD at low and high momenta.
The extraction of the LO contact term, in a given scheme, then proceeds by
matching the full amplitude obtained in this way to the LO chiral EFT
amplitude. We considered several sources of uncertainty, chief among them the
missing contributions from inelastic intermediate states, such as $N\\!N\pi$.
The intermediate steps of our analysis have been performed in dimensional
regularization with the $\overline{\text{MS}}$ scheme, while the final result
can be expressed in terms of the scheme-independent renormalized amplitude
$\mathcal{A}_{\nu}(|\mathbf{p}|,|\mathbf{p}^{\prime}|)$ at a set of kinematic
points near threshold, where LO chiral EFT is expected to give an excellent
approximation. Using our synthetic data, as discussed in Section 7, one can
then determine the contact term in any regularization and renormalization
scheme, in particular the ones employed in nuclear-structure calculations for
isotopes of experimental interest for $0\nu\beta\beta$ searches. This
application is timely in view of the remarkable progress in ab-initio
calculations of $0\nu\beta\beta$ decay rates of light and intermediate-mass
nuclei [43, 29, 30, 31, 109], ranging from 6He to 48Ca and 76Ge, starting from
microscopic nuclear forces obtained from chiral EFT. So far these decay rates
include the long-distance neutrino-exchange contributions while omitting the
contact term, which can now be remedied using our synthetic data for $nn\to
ppe^{-}e^{-}$, allowing, for the first time, for complete LO calculations of
nuclear $0\nu\beta\beta$ decay rates. For heavier nuclei such as 136Xe, which
are still beyond the capabilities of ab-initio techniques, the impact of the
short-range term could be studied indirectly, e.g., by matching ab-initio
results and nuclear models for nuclei accessible to both approaches (this
strategy was used in Ref. [110] for the axial-vector current).

As a benchmark point, we find for the amplitude ${\cal A}_{\nu}$ in the
$\overline{\text{MS}}$ scheme at $\mu_{\chi}=4M_{\pi}$, with initial- and
final-state nucleon momenta $|\mathbf{p}|=25\,\text{MeV}$ and
$|\mathbf{p}^{\prime}|=30\,\text{MeV}$, that the contact-term contribution
adds destructively to the neutrino exchange at the 15% level. However, as
discussed in detail in Section 7, this statement depends on the
renormalization scheme and scale, e.g., the two cutoff schemes studied there
can lead to both constructive and destructive effects depending on the choice
of scale. This illustrates that the separation into short- and long-distance
contributions is unphysical, see Fig. 15 for the decomposition in the
$\Lambda_{p}$ and $R_{S}$ schemes as a function of the final-state momentum
and the cutoff scales. Further, as discussed in Refs. [40, 43], while a
contact term of natural size affects $\Delta I=0$ transitions such as $nn\to
ppe^{-}e^{-}$ at the (10–20)% level, its effect is amplified to the level of
(50–70)% in $\Delta I=2$ nuclear transitions due to a node in the matrix
element density. The size of the effect in realistic $0\nu\beta\beta$ nuclear
transitions can now be addressed, greatly reducing a crucial uncertainty in
the interpretation of future experimental searches [111, 112, 113, 114, 115,
116, 117].

Improving the accuracy of the results presented here would require at least a
thorough study of the inelastic $N\\!N\pi$ channel, but at the same time a
more systematic connection to the Cottingham formula would likely become
necessary as well: while the identification of the leading elastic effects is
rather intuitive, the extension to subleading corrections is not. Additional
experimental input could, in principle, separate ${\cal C}_{1}$ from ${\cal
C}_{2}$, e.g., via CIB in nuclei, but such an extraction also requires the
development of a suitable theoretical framework. While our results allow for
first phenomenological estimates of the impact of the contact term on
$0\nu\beta\beta$ decay rates, they thus also define a benchmark for future
lattice-QCD calculations [52, 53, 46, 54, 55, 56, 57]. In addition to
comparing the final result for ${\cal A}_{\nu}$, there could also be aspects
of the matching strategy and spectral representation, as described in Section
2, that might prove synergistic between the two approaches. Finally, having
concentrated on the LO contact term for light Majorana exchange in this paper,
we remark that contact-term contributions arise at LO for other operators
mediating $0\nu\beta\beta$ decay as well, both at dimension $7$ and dimension
$9$ [37], or through exchange of massive sterile neutrinos [44], and
generalizations of the strategies presented here could help constrain the
contact terms in these mechanisms.

Acknowledgments—We thank Jon Engel, Evgeny Epelbaum, Michael Graesser, Bira
van Kolck, and André Walker-Loud for discussions at various stages of this
work. The work of VC and EM is supported by the US Department of Energy
through the Los Alamos National Laboratory. Los Alamos National Laboratory is
operated by Triad National Security, LLC, for the National Nuclear Security
Administration of U.S. Department of Energy (Contract No. 89233218CNA000001).
WD is supported by U.S. Department of Energy Office of Science, under contract
DE-SC0009919. MH is supported by an Eccellenza Grant (Project No.
PCEFP2_181117) of the Swiss National Science Foundation. JdV is supported by
the RHIC Physics Fellow Program of the RIKEN BNL Research Center. The DOE DBD
Topical Nuclear Theory Collaboration inspired this work and partially
supported some of us during its genesis. MH thanks the T-2 group at LANL for
their hospitality and support during a visit when this project was initiated.

##  Appendix A Half-off-shell $\boldsymbol{T}$ matrix

In this appendix we describe the calculation of the HOS $T$-matrix elements
defined in Eq. (4.8) and the associated form factors in Eq. (4.17). These
quantities affect the $nn\to pp$ amplitude through the ratio
$r(|{\mathbf{k}}|)$, introduced in Eq. (4.19), parameterizing the higher-order
corrections to the forward Compton amplitude. We present the calculations in
order of increasing difficulty, starting from ${/\pi}$EFT, in which analytic
results can be obtained, moving to chiral EFT, and finally using $N\\!N$
interactions from potential models that successfully fit the $N\\!N$
scattering data.

###  A.1 The half-off-shell form factor in pionless EFT

In this appendix we discuss the HOS form factors in the framework of
${/\pi}$EFT. Since analytic results are readily available, the ${/\pi}$EFT
analysis will provide several insights on the general problem at hand.

####  A.1.1 Half-off-shell $\boldsymbol{T}$ matrix in pionless EFT at NLO

Working in dimensional regularization with power divergence subtraction [25],
the resummation of bubble diagrams with $N\\!N$ vertices involving $C_{0}$ and
$C_{2}$ leads to ($E={\mathbf{p}}^{2}/m_{N}$,
${\mathbf{q}}^{2}\neq{\mathbf{p}}^{2}$)

$\langle{\mathbf{q}}|\hat{T}_{S}(E)|{\mathbf{p}}\rangle=\frac{\tilde{C}(E)+\frac{C_{2}}{2}({\mathbf{q}}^{2}-{\mathbf{p}}^{2})}{1-\tilde{C}(E)I_{0}}\,,\qquad\tilde{C}(E)=C_{0}+C_{2}\,{\mathbf{p}}^{2}\,,$
(A.1)

where we used the definitions (3.3) with $\hat{V}_{\pi}\to 0$ and

$I_{0}\left(E\right)=\Big{(}\frac{\mu}{2}\Big{)}^{4-d}\int\frac{d^{d-1}q}{(2\pi)^{d-1}}\,\frac{m_{N}}{{\mathbf{p}}^{2}-{\mathbf{q}}^{2}+i\epsilon}=-\frac{m_{N}}{4\pi}\,(\mu+i|{\mathbf{p}}|)\,.$
(A.2)

This leads to

$f_{S}({\mathbf{q}},{\mathbf{p}})=1+\frac{C_{2}}{2C_{0}}\frac{{\mathbf{q}}^{2}-{\mathbf{p}}^{2}}{1+(C_{2}/C_{0}){\mathbf{p}}^{2}}\simeq
1+\frac{C_{2}}{2C_{0}}\big{(}{\mathbf{q}}^{2}-{\mathbf{p}}^{2}\big{)}\,.$
(A.3)

As mentioned before, these HOS form factors are different depending on whether
they involve the initial- or final-state on-shell momenta, cf. Eq. (2). By
themselves these form factors are clearly not physical. First, one notes that
they depend on the renormalization scale $\mu$, as $C_{2}/C_{0}$ is
$\mu$-dependent. The form factors vary significantly with $\mu$, while having
the expected qualitative features at small $|{\mathbf{k}}|$. Second, as shown
below, $f_{S}$ can be changed by performing field redefinitions, as expected
for any off-shell quantity. However, in the actual application insertions of
$f_{S}$ combine into a physical contribution to the $nnW^{+}(k)\to ppW^{-}(k)$
amplitude, or equivalently to the integrand in the sum of ${\cal
M}^{<}_{B+\bar{B}}$ and ${\cal M}^{<}_{C}$ in Eq. (4.2). In the end, the net
effect is a modification of the function ${\cal I}_{C}^{<}(|{\mathbf{k}}|)$
and hence a change in $r(|{\mathbf{k}}|)$.

####  A.1.2 Behavior of $\boldsymbol{f_{S}}$ and $\boldsymbol{{\cal
I}_{C}^{<}(|{\mathbf{k}}|)}$ under field redefinitions

We consider the non-linear field redefinition given in Ref. [79], namely

$N\to N+\eta N(N^{\dagger}N)\,,$ (A.4)

where $\eta$ is an arbitrary parameter of mass dimension $-3$. This
transformation induces a derivative $N\\!N$ interaction from the kinetic term
as well as a two-nucleon current operator, schematically

$\displaystyle{\cal L}$ $\displaystyle\to{\cal L}\ +\
\eta\,(N^{\dagger}N)\,N^{\dagger}\left(i\partial_{0}+\frac{{\mathbf{\partial}}^{2}}{2m_{N}}\right)N\
+\
\eta\,N^{\dagger}\left(i\partial_{0}+\frac{{\mathbf{\partial}}^{2}}{2m_{N}}\right)\left[N\,(N^{\dagger}N)\right]\,,$
$\displaystyle N^{\dagger}J_{\mu}\tau^{+}N$ $\displaystyle\to
N^{\dagger}J_{\mu}\tau^{+}N+2\eta\
N^{\dagger}J_{\mu}\tau^{+}N\,(N^{\dagger}N)\,.$ (A.5)

The resulting $\eta$-dependent four-nucleon vertex is [79]

$\langle N_{1^{\prime}}N_{2^{\prime}}|i{\cal
L}|N_{1}N_{2}\rangle=i\eta\left[\left(\omega_{1}-\frac{{\mathbf{p}}_{1}^{2}}{2m_{N}}\right)+\left(\omega_{2}-\frac{{\mathbf{p}}_{2}^{2}}{2m_{N}}\right)+\left(\omega_{1^{\prime}}-\frac{{\mathbf{p}}_{1^{\prime}}^{2}}{2m_{N}}\right)+\left(\omega_{2^{\prime}}-\frac{{\mathbf{p}}_{2^{\prime}}^{2}}{2m_{N}}\right)\right]\,.$
(A.6)

This modifies the HOS form factors (A.3) as follows

$f_{S}({\mathbf{q}},{\mathbf{p}})\simeq
1+\frac{C_{2}+2(\eta/m_{N})}{2C_{0}}\big{(}{\mathbf{q}}^{2}-{\mathbf{p}}^{2}\big{)}\,.$
(A.7)

These $\eta$-dependent HOS form factors combine with new contributions to the
current matrix element in such a way to leave no $\eta$-dependent term in the
amplitude $nnW^{+}(k)\to ppW^{-}(k)$ itself. It is simple to verify that
cancellations happen through groups of diagrams represented in Fig. 16.
Therefore, the integrands in the amplitudes ${\cal M}^{<}_{B+\bar{B}}$ and
${\cal M}^{<}_{C}$ in Eq. (4.2) (and hence ${\cal I}_{C}^{<}(|{\mathbf{k}}|)$
via Eqs. (4.7) and (4.18)) are individually invariant under the field
redefinition considered here, i.e., they do not depend on $\eta$. We have used
this simple field redefinition of Ref. [79] to illustrate a more generally
valid point.

Figure 16:  Group of diagrams contributing to $nnW^{+}(k)\to ppW^{-}(k)$ in
which the effects of field redefinitions cancels out. The red circles denote
vertices proportional to $\eta$, as given in Eq. (A.5). The left-most legs
denote on-shell protons. The right-most legs denote off-shell neutrons. The
wavy lines denote weak current insertions.

####  A.1.3 Impact of $\boldsymbol{f_{S}}$ on $\boldsymbol{{\cal
I}_{C}^{<}(|{\mathbf{k}}|)}$

Inserting the form factor (A.3) into Eq. (4.18) one obtains (using
$|{\mathbf{p}}^{\prime}|=|{\mathbf{p}}|$)

${\cal I}_{C}^{<}(|{\mathbf{k}}|)=\frac{1}{8|{\bf k}|}\,\theta(|{\bf
k}|-2|{\bf p}|)+\frac{i}{8\pi|{\bf k}|}\,\log\left|\frac{1+2\frac{|{\bf
p}|}{{|{\bf k}|}}}{1-2\frac{|{\bf p}|}{{|{\bf
k}|}}}\right|+\frac{1}{4\pi}\,\frac{C_{2}}{C_{0}}\Big{(}\mu+i|{\mathbf{p}}|\Big{)}+\mathcal{O}(C_{2}^{2})\,.$
(A.8)

For the real part of ${\cal I}_{C}^{<}(|{\mathbf{k}}|)$, relevant for the
matching, this implies

${\rm Re}\,{\cal I}_{C}^{<}(|{\mathbf{k}}|)=\frac{1}{8|{\bf
k}|}\left[1+\frac{|{\mathbf{k}}|}{\pi}\left(\frac{2\mu
C_{2}}{C_{0}}\right)\right]\approx\frac{1}{8|{\bf
k}|}\left[1-\frac{|{\mathbf{k}}|r_{0}}{\pi}\right]\,.$ (A.9)

In the last step we have used the ${/\pi}$EFT expressions of $C_{0,2}(\mu)$
and taken $\mu\sim M_{\pi}\gg 1/a$. The size of the effective range
$r_{0}=2.73~{}{\rm fm}=1/(72\,\text{MeV})$ implies sizable deviations from the
leading $1/(8|{\mathbf{k}}|)$ behavior already at $|{\mathbf{k}}|\leq
M_{\pi}$.

While the analysis above is promising, there are two unsatisfactory aspects:
(i) Even though it does not enter the matching, the imaginary part of ${\cal
I}_{C}^{<}(|{\mathbf{k}}|)$ depends on the scale $\mu$ through the ratio
$C_{2}/C_{0}$. (ii) The real part of ${\cal I}_{C}^{<}(|{\mathbf{k}}|)$ is
only approximately independent of $\mu$, as we have neglected $1/(\mu a)\ll 1$
in the expression

$\frac{\mu C_{2}}{C_{0}}=-\frac{r_{0}}{2}\frac{1}{1-\frac{1}{\mu a}}\,.$
(A.10)

Both issues can be addressed by studying how $f_{S}$ impacts all components of
the generalized Compton amplitude, i.e., the integrand of ${\cal
M}^{<}_{B+\bar{B}}$ and ${\cal M}^{<}_{C}$ in Eq. (4.2), and not just the
integrand of ${\cal M}^{<}_{C}$ as done above. The final result is very
simple: the contributions combine to give a real-valued and scale-independent
shift to ${\cal I}_{C}^{<}(|{\mathbf{k}}|)$ equal to $\delta{\cal
I}_{C}^{<}(|{\mathbf{k}}|)=-r_{0}/(8\pi)$, which agrees with the result in the
second line of Eq. (A.9).

The argument proceeds as follows. Omitting a common factor proportional to the
weak current insertions, namely $-(1+3g_{A}^{2})/{\mathbf{k}}^{2}$, the
integrand corresponding to ${\cal M}^{<}_{B+\bar{B}}+{\cal M}^{<}_{C}$ in Eq.
(4.2) has the following structure:

$\displaystyle
m_{N}\Bigg{(}T(E^{\prime})\,\frac{f_{S}({\mathbf{p}}^{\prime},{\mathbf{k}}+{\mathbf{p}})}{{\mathbf{p}}^{\prime
2}-({\mathbf{k}}+{\mathbf{p}})^{2}+i\epsilon}+\frac{f_{S}({\mathbf{k}}+{\mathbf{p}}^{\prime},{\mathbf{p}})}{{\mathbf{p}}^{2}-({\mathbf{k}}+{\mathbf{p}}^{\prime})^{2}+i\epsilon}\,T(E)$
$\displaystyle+T(E^{\prime})\,m_{N}\,\int\frac{d^{d-1}q}{(2\pi)^{d-1}}\,\frac{f_{S}({\mathbf{p}}^{\prime},{\mathbf{k}}+{\mathbf{q}})}{{\mathbf{p}}^{\prime
2}-({\mathbf{k}}+{\mathbf{q}})^{2}+i\epsilon}\
\frac{f_{S}({\mathbf{q}},{\mathbf{p}})}{{\mathbf{p}}^{2}-{\mathbf{q}}^{2}+i\epsilon}\
T(E)\Bigg{)}\,,$ (A.11)

where $T(E)$ denotes the LO on-shell $N\\!N$ $T$ matrix in the ${}^{1}S_{0}$
channel. Substituting the form factors (A.3), one obtains

$-m_{N}\Bigg{(}\frac{C_{2}}{2C_{0}}\Big{(}T(E^{\prime})+T(E)\Big{)}+\frac{C_{2}}{2C_{0}^{2}}\Big{[}T(E^{\prime})\
C_{0}I_{0}(E)T(E)\ +\ T(E^{\prime})C_{0}I_{0}(E^{\prime})\
T(E)\Big{]}\Bigg{)}\,,$ (A.12)

where $I_{0}(E)$ is given in Eq. (A.2) and $T(E)=C_{0}/(1-C_{0}I_{0}(E))$.
Using the relation

$T(E)-C_{0}=C_{0}I_{0}(E)T(E)\,,$ (A.13)

one sees that the terms proportional to $C_{2}/C_{0}$ in Eq. (A.12) cancel and
one is left with the scale-independent quantity $T(E^{\prime})\
(-C_{2}/C_{0}^{2})\ T(E)$. Upon matching the overall factors of $m_{N}$, this
leads to

$\delta{\cal
I}_{C}^{<}(|{\mathbf{k}}|)=-\frac{C_{2}}{m_{N}C_{0}^{2}}=-\frac{r_{0}}{8\pi}\,.$
(A.14)

To conclude, we discuss several lessons that can be abstracted from the
${/\pi}$EFT analysis: first, the large fractional correction in ${\cal
I}_{C}^{<}(|{\mathbf{k}}|)$ (and hence $a_{<}(|{\mathbf{k}}|)$) at relatively
low $|{\mathbf{k}}|\sim M_{\pi}$ (in ${/\pi}$EFT) is due to the large
effective range in the ${}^{1}S_{0}$ channel. Second, in ${/\pi}$EFT, the
inclusion of HOS form factor effects from ${\cal M}^{<}_{B+\bar{B}}$ was
needed to cancel a subleading $\mathcal{O}(1/(\mu a))$ scale dependence in the
integrand corresponding to ${\cal M}^{<}_{C}$, proportional to ${\cal
I}_{C}^{<}(|{\mathbf{k}}|)$. In practice one gets the leading corrections by
just including the HOS form factors in the ${\cal M}^{<}_{C}$ integrand. This
is consistent with the previous result that ${\cal M}^{<}_{B+\bar{B}}$ drops
out of the LO matching formula and the expectation that the HOS form factor
makes the integrand in ${\cal M}^{<}_{B+\bar{B}}$ even more convergent.
Accordingly, this suggests that the inclusion of the HOS form factors in the
integrand of ${\cal M}^{<}_{C}$ captures the bulk of the physical suppression
of the $nnW^{+}(k)\to ppW^{-}(k)$ amplitude at small $|{\mathbf{k}}|$, while
the inclusion of HOS form factors in ${\cal M}^{<}_{B+\bar{B}}$ would lead to
effects that are below the current level of accuracy.

###  A.2 Half-off-shell effects in chiral EFT at NLO

In this appendix we extend the analysis to chiral EFT, thus keeping pion-
exchange effects. To take into account off-shell effects in $\hat{T}_{S}$, we
start from Eq. (3.3),

$\hat{T}_{S}^{\chi}=(1+\hat{V}_{\pi}\hat{G}^{(\pi)}_{+}(E))\hat{V}_{S}(1-\hat{G}^{(\pi)}_{+}(E)\hat{V}_{S})^{-1}(1+\hat{G}^{(\pi)}_{+}(E)\hat{V}_{\pi})\,,$
(A.15)

where
$\langle{\mathbf{q}}|\hat{V}_{S}|{\mathbf{p}}\rangle=C+C_{2}\frac{{\mathbf{q}}^{2}+{\mathbf{p}}^{2}}{2}$
is the NLO short-range potential, so that we can write

$\hat{V}_{S}=(1+m_{N}EC_{2}/C)\,\hat{V}(C)-\frac{m_{N}C_{2}}{2C}\left\\{\hat{V}(C),\,\left(\hat{G}^{(0)}_{+}(E)\right)^{-1}\right\\}\,.$
(A.16)

Here $\hat{V}(C)$ is the LO part of the short-distance potential and the last
term gives rise to off-shell effects proportional to $C_{2}$. Expanding
$\hat{T}_{S}^{\chi}$ to first order in $C_{2}$ we have

$\displaystyle\langle{\mathbf{q}}|\hat{T}_{S}^{\chi}|{\mathbf{k}}\rangle$
$\displaystyle=\chi^{+}_{{\mathbf{q}}}(0)\chi^{+}_{{\mathbf{k}}}(0)\left[K_{E}^{(0)}+K_{E}^{(1)}\right]-\frac{C_{2}}{C}K_{E}^{(0)}\left[\chi^{+}_{{\mathbf{q}}}(0)\frac{{\mathbf{p}}^{2}-{\mathbf{k}}^{2}}{2}+\chi^{+}_{{\mathbf{k}}}(0)\frac{{\mathbf{p}}^{2}-{\mathbf{q}}^{2}}{2}\right]\,,$
(A.17) $\displaystyle K_{E}^{(0)}$
$\displaystyle=\frac{C}{1-CG^{(\pi)}_{E}(0,0)}\,,\qquad
K_{E}^{(1)}=\left(K_{E}^{(0)}\right)^{2}\frac{C_{2}}{C^{2}}\left[{\mathbf{p}}^{2}-m_{N}\Big{(}\frac{\mu}{2}\Big{)}^{4-d}\int\frac{d^{d-1}q}{(2\pi)^{d-1}}V_{\pi}(q)\right]\,,$

with ${\mathbf{p}}^{2}=m_{N}E$. The scale-dependent term in $K_{E}^{(1)}$ can
be absorbed into an NLO correction to $C\to C+C^{(1)}$, after which
$K_{E}^{(0,1)}$ become scale independent and the NLO on-shell matrix element
of $\hat{T}_{S}^{\chi}$ takes the form
$T_{S}^{\chi}({\mathbf{p}},{\mathbf{p}})=\left(\chi^{+}_{{\mathbf{p}}}(0)\right)^{2}\left[K_{E}^{(0)}+K_{E}^{(1)}\right]$.

Inserting the above off-shell matrix elements
$\langle{\mathbf{q}}|\hat{T}_{S}^{\chi}|{\mathbf{k}}\rangle$ into Eq. (4.2) we
can write the low-momentum component of the $\mathcal{M}_{C}$ amplitude as

$\displaystyle\mathcal{M}_{C}^{<}$
$\displaystyle=\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\
\langle{\mathbf{p}}^{\prime}|\left(\hat{T}^{\chi}_{S}(E^{\prime})\hat{G}_{+}^{(0)}(E^{\prime})\right)\,\hat{O}_{\chi}^{LL}({\mathbf{k}})\,\left(\hat{G}_{+}^{(0)}(E)\,\hat{T}^{\chi}_{S}(E)\right)|{\mathbf{p}}\rangle$
$\displaystyle=\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}f({\mathbf{k}})\chi^{+}_{{\mathbf{p}}}(0)\chi^{+}_{{\mathbf{p}}^{\prime}}(0)$
$\displaystyle\times\Bigg{\\{}\left[K_{E^{\prime}}^{(0)}+K_{E^{\prime}}^{(1)}\right]\left[K_{E}^{(0)}+K_{E}^{(1)}\right]\int\frac{d^{3}{\mathbf{q}}}{(2\pi)^{3}}\frac{\chi^{+}_{{\mathbf{q}}}(0)}{E^{\prime}-{\mathbf{q}}^{2}/m_{N}+i\epsilon}\frac{\chi^{+}_{{\mathbf{q}}+{\mathbf{k}}}(0)}{E-({\mathbf{q}}+{\mathbf{k}})^{2}/m_{N}+i\epsilon}$
$\displaystyle\qquad-\frac{m_{N}C_{2}}{2C}K_{E}^{(0)}K_{E^{\prime}}^{(0)}\left(G^{(\pi)}_{E}(0,0)+G^{(\pi)}_{E^{\prime}}(0,0)\right)\Bigg{\\}}\,,$
(A.18)

where
$\langle{\mathbf{q}}|\hat{O}_{\chi}^{LL}({\mathbf{k}})|{\mathbf{q}}^{\prime}\rangle=(2\pi)^{3}\delta^{3}({\mathbf{q}}-{\mathbf{q}}^{\prime}+{\mathbf{k}})f({\mathbf{k}})$,
with
$f({\mathbf{k}})=-\frac{g^{2}_{V}({\mathbf{k}})+3h_{GT}({\mathbf{k}})}{{\mathbf{k}}^{2}}$.
The first term in curly brackets contributes to the most singular part of
$\mathcal{M}_{C}$ (the parts that diverge when integrating over
${\mathbf{k}}$). In particular, the terms with
$\chi_{{\mathbf{q}},{\mathbf{q}}+{\mathbf{k}}}^{+}(0)\to 1$ give rise to the
part of $\mathcal{I}_{C}$ in Eq. (4.18), with $f_{S}\to 1$. In other words,
these pieces are responsible for the contributions involving the on-shell
$T_{S}$ matrix. Instead, the contributions from the terms proportional to
$1-\chi_{{\mathbf{q}},{\mathbf{q}}+{\mathbf{k}}}^{+}(0)$ correspond to
diagrams involving pion exchanges within the bubble, which are convergent.

The last line in Eq. (A.18) arises from the off-shell piece of $\hat{T}_{S}$,
the second term in Eq. (A.17), and also contributes to the most divergent
parts of $\mathcal{M}_{C}$. As in pionless EFT, these terms are $\mu$
dependent by themselves and have to be combined with pieces from the
$\mathcal{M}_{B}$ diagrams in order to obtain a $\mu$-independent answer. This
part of the amplitude can be written as

$\displaystyle\mathcal{M}^{<}_{B+\bar{B}}$
$\displaystyle=\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\Bigg{[}\langle{\mathbf{p}}^{\prime}|\left(\hat{T}^{\chi}_{S}(E^{\prime})\,\hat{G}_{+}^{(0)}(E^{\prime})\right)\,\hat{O}_{\chi}^{LL}({\mathbf{k}})\,\left(I+\hat{G}_{+}^{(0)}(E)\,\hat{T}^{\chi}_{\pi}(E)\right)|{\mathbf{p}}\rangle$
$\displaystyle\qquad+\langle{\mathbf{p}}^{\prime}|\left(I+\hat{T}^{\chi}_{\pi}(E^{\prime})\hat{G}_{+}^{(0)}(E)\right)\,\hat{O}_{\chi}^{LL}({\mathbf{k}})\,\left(\hat{G}_{+}^{(0)}(E)\hat{T}^{\chi}_{S}(E)\right)|{\mathbf{p}}\rangle\Bigg{]}$
$\displaystyle=\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}f({\mathbf{k}})\Bigg{\\{}\frac{d^{3}{\mathbf{q}}}{(2\pi)^{3}}\Bigg{[}\frac{\langle\chi_{{\mathbf{p}}^{\prime}}|{\mathbf{q}}\rangle\chi^{+}_{{\mathbf{q}}+{\mathbf{k}}}(0)}{E-({\mathbf{q}}+{\mathbf{k}})^{2}/m_{N}+i\epsilon}\chi^{+}_{{\mathbf{p}}}(0)\left(K_{E}^{(0)}+K_{E}^{(1)}\right)$
$\displaystyle\qquad+\left(K_{E^{\prime}}^{(0)}+K_{E^{\prime}}^{(1)}\right)\chi^{+}_{{\mathbf{p}}^{\prime}}(0)\frac{\chi^{+}_{{\mathbf{q}}}(0)\langle{\mathbf{q}}+{\mathbf{k}}|\chi_{{\mathbf{p}}}\rangle}{E^{\prime}-{\mathbf{q}}^{2}/m_{N}+i\epsilon}\Bigg{]}$
$\displaystyle-\frac{m_{N}C_{2}}{2C}\chi^{+}_{{\mathbf{p}}^{\prime}}(0)\chi^{+}_{{\mathbf{p}}}(0)\left(K_{E}^{(0)}+K_{E^{\prime}}^{(0)}\right)\Bigg{\\}}\,,$
(A.19)

where the terms within square brackets are finite, while the last line comes
from the off-shell part of $\hat{T}_{S}^{\chi}$ , which is again $\mu$
dependent. Combining the terms that have explicit $C_{2}$ dependence with
those in $\mathcal{M}_{C}$, we obtain

$\displaystyle\mathcal{M}^{<}_{B+\bar{B}}+\mathcal{M}^{<}_{C}\Big{|}_{C_{2}}$
$\displaystyle=-\frac{m_{N}C_{2}}{2C}\chi^{+}_{{\mathbf{p}}^{\prime}}(0)\chi^{+}_{{\mathbf{p}}}(0)\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}f({\mathbf{k}})$
$\displaystyle\qquad\times\left[K_{E}^{(0)}+K_{E^{\prime}}^{(0)}+K_{E}^{(0)}K_{E^{\prime}}^{(0)}\left(G^{(\pi)}_{E}(0,0)+G^{(\pi)}_{E^{\prime}}(0,0)\right)\right]$
$\displaystyle=-\frac{m_{N}C_{2}}{C^{2}}\chi^{+}_{{\mathbf{p}}^{\prime}}(0)\chi^{+}_{{\mathbf{p}}}(0)K_{E}^{(0)}K_{E^{\prime}}^{(0)}\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}f({\mathbf{k}})\,,$
(A.20)

which is now $\mu$ independent and has the factors of
$\chi^{+}_{{\mathbf{p}}}(0)K_{E}$ we expect for the $\mathcal{M}_{C}$-type
diagrams. Using

$\mathcal{M}_{C}^{<}=\chi^{+}_{{\mathbf{p}}^{\prime}}(0)K_{E^{\prime}}^{(0)}\left[m_{N}^{2}\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}f({\mathbf{k}})\mathcal{I}_{C}\right]\chi^{+}_{{\mathbf{p}}}(0)K_{E}^{(0)}\,,$
(A.21)

the terms in Eq. (A.20) provide a shift to $\mathcal{I}^{<}_{C}$ equal to
$\delta\mathcal{I}^{<}_{C}=-C_{2}/(m_{N}C^{2})$, which is the same as in
${/\pi}$EFT, up to the fact that $C_{2}/C^{2}$ cannot simply be expressed in
terms of the effective range parameter $r_{0}$.

###  A.3 Half-off-shell effects in $\boldsymbol{N\\!N}$ potential models

Finally, in this appendix we describe the calculation of the HOS form factor
and the ratio $r(|{\mathbf{k}}|)$ in potential models. As a benchmark
resonance model for the ${}^{1}S_{0}$ channel interaction we will use the
three-Yukawa potential, known as the Reid potential [88], which has also been
studied in the context of analyzing the convergence of the $N\\!N$ EFT [87].
This is a simplified version of the general resonance saturation model for
$N\\!N$ interactions [118]. In coordinate space the potential takes the form

$V=V_{\pi}+V_{S}\,,\qquad
V_{\pi}(r)=-\alpha_{\pi}\frac{e^{-M_{\pi}r}}{r}\,,\qquad
V_{S}(r)=-\alpha_{\sigma}\frac{e^{-M_{\sigma}r}}{r}+\alpha_{\rho}\frac{e^{-M_{\rho}r}}{r}\,,$
(A.22)

with the short-range term modeled through the attractive $\sigma$-meson and
repulsive $\rho$-meson contributions. Once the mass parameters are fixed
($M_{\sigma}=4M_{\pi}$, $M_{\rho}=7M_{\pi}$ in Ref. [88] and
$M_{\sigma}=500\,\text{MeV}$, $M_{\rho}=770\,\text{MeV}$ in Ref. [87]), the
couplings $\alpha_{\sigma,\rho}$ are tuned to reproduce the ${}^{1}S_{0}$
phase shift or the scattering length and effective range [87], producing an
overall good representation of the data.

In this class of models it is relatively simple to compute the HOS $T$-matrix
elements defined in Eq. (4.8) and the associated form factors in Eq. (4.17):
one needs to solve the LS equation or alternatively the Schrödinger equation
using $V_{S}$ in Eq. (A.22). We have performed the calculation in both
methods, LS equation in momentum space and Schrödinger equation in coordinate
space. Useful cross-checks in the calculation are provided by the property
that the HOS $T$ matrix has the same phase-shift dependence as the on-shell
$T$ matrix $T({\mathbf{p}},{\mathbf{p}})$ [86]

$T({\mathbf{q}},{\mathbf{p}})=|T({\mathbf{q}},{\mathbf{p}})|\
e^{\delta(|{\mathbf{p}}|)}\,.$ (A.23)

We have also performed the calculation of $T({\mathbf{q}},{\mathbf{p}})$ and
$f_{S}({\mathbf{q}},{\mathbf{p}})$ with the AV18 potential [89], after
subtracting the one-pion exchange. This is a representative of the class of
high-quality potentials that reproduce the $N\\!N$ scattering data for
$|{\mathbf{p}}|$ out to several hundred MeV.

Figure 17:  Left panel: HOS form factor $f_{S}({\mathbf{q}},{\mathbf{p}})$ at
$|{\mathbf{p}}|=1\,\text{MeV}$ for Kaplan–Steele (blue), Reid (green, dotted),
and AV18 (red, dashed) potential. Right panel: $r(|{\mathbf{k}}|)$ computed at
$|{\mathbf{p}}|=|{\mathbf{p}}^{\prime}|=1\,\text{MeV}$ for Kaplan–Steele
(blue), Reid (green, dotted), and AV18 (red, dashed) potential.

Representative results from this analysis are summarized in Fig. 17. In the
left panel we show $f_{S}({\mathbf{q}},{\mathbf{p}})$ versus $|{\mathbf{q}}|$
at $|{\mathbf{p}}|=1\,\text{MeV}$ for the three potentials considered here,
Kaplan–Steele, Reid, and AV18. In the right panel we report the ratio
$r(|{\mathbf{k}}|)$, see Eq. (4.19), relevant for the calculation of ${\cal
A}_{C}^{<,\rm sing}$. The difference in $r(|{\mathbf{k}}|)$ provides some
indication of the level of model dependence in our matching calculation. While
locally the differences can be large, the fall-off of the ${\cal A}_{C}^{<,\rm
sing}$ integrand multiplying $r(|{\mathbf{k}}|)$ dilutes the impact of any
feature in $r(|{\mathbf{k}}|)$ above $500\,\text{MeV}$, so that ultimately the
impact on the extracted LEC is relatively minor. We expect this conclusion to
hold in general, beyond the limited set of potential models explored here.

##  Appendix B Left–right correlator

In the main text we defined the “left–left” correlator, which naturally
appears when dealing with the weak currents. In the case of electromagnetic
currents and for matching purposes, we also need the “left–right” correlator.
We define

$\displaystyle\hat{\Pi}_{\mu\nu}^{LL}(k,0)$ $\displaystyle=\int
d^{4}r\,e^{ik\cdot r}\ T\Big{(}\bar{u}_{L}\gamma_{\mu}d_{L}(r/2)\
\bar{u}_{L}\gamma_{\mu}d_{L}(-r/2)\Big{)}\,,$
$\displaystyle\hat{\Pi}_{\mu\nu}^{LR}(k,0)$ $\displaystyle=\int
d^{4}r\,e^{ik\cdot r}\ T\Big{(}\bar{u}_{L}\gamma_{\mu}d_{L}(r/2)\
\bar{u}_{R}\gamma_{\mu}d_{R}(-r/2)\Big{)}\,,$ (B.1)

and

$\displaystyle\hat{O}^{LL}_{\alpha\beta}(0)$
$\displaystyle=\hat{O}_{\alpha\beta}(0)=\bar{u}_{L}(0)\gamma_{\alpha}T^{a}d_{L}(0)\
\bar{u}_{L}(0)\gamma_{\beta}T^{a}d_{L}(0)\,,$
$\displaystyle\hat{O}^{LR}_{\alpha\beta}(0)$
$\displaystyle=\bar{u}_{L}(0)\gamma_{\alpha}T^{a}d_{L}(0)\
\bar{u}_{R}(0)\gamma_{\beta}T^{a}d_{R}(0)\,.$ (B.2)

With these definitions, the time-ordered products at large $k^{2}$ can be
written as, see Eq. (4.24),

$\displaystyle\hat{\Pi}_{\mu\nu}^{LL}(k,0)$
$\displaystyle=-\frac{4ig_{s}^{2}}{(k^{2}+i\epsilon)^{3}}\times\Big{[}(k_{\mu}k_{\nu}-k^{2}g_{\mu\nu})g^{\alpha\beta}\,\hat{O}^{LL}_{\alpha\beta}(0)$
$\displaystyle+g_{\mu\nu}k^{\alpha}k^{\beta}\hat{O}^{LL}_{\alpha\beta}(0)+k^{2}\hat{O}^{LL}_{\nu\mu}(0)-k_{\nu}k^{\alpha}\,\hat{O}^{LL}_{\alpha\mu}(0)-k_{\mu}k^{\alpha}\,\hat{O}^{LL}_{\nu\alpha}(0)\Big{]}\,,$
$\displaystyle\hat{\Pi}_{\mu\nu}^{LR}(k,0)$
$\displaystyle=\frac{4ig_{s}^{2}}{(k^{2}+i\epsilon)^{3}}\times\Big{[}(k_{\mu}k_{\nu}-k^{2}g_{\mu\nu})g^{\alpha\beta}\,\hat{O}^{LR}_{\alpha\beta}(0)$
$\displaystyle+g_{\mu\nu}k^{\alpha}k^{\beta}\hat{O}^{LR}_{\alpha\beta}(0)+k^{2}\hat{O}^{LR}_{\nu\mu}(0)-k_{\nu}k^{\alpha}\,\hat{O}^{LR}_{\alpha\mu}(0)-k_{\mu}k^{\alpha}\,\hat{O}^{LR}_{\nu\alpha}(0)\Big{]}\,.$
(B.3)

Similarly, we obtain:

$\displaystyle g^{\mu\nu}\,\hat{\Pi}_{\mu\nu}^{LL}(k,0)$
$\displaystyle=\frac{-4ig_{s}^{2}}{(k^{2}+i\epsilon)^{3}}\
\Big{(}-2k^{2}(\hat{O}^{LL})_{\alpha}^{\alpha}+2k^{\alpha}k^{\beta}\hat{O}^{LL}_{\alpha\beta}\Big{)}\to\frac{6ig_{s}^{2}}{(k^{2}+i\epsilon)^{2}}(\hat{O}^{LL})_{\alpha}^{\alpha}\,,$
$\displaystyle g^{\mu\nu}\,\hat{\Pi}_{\mu\nu}^{LR}(k,0)$
$\displaystyle=\frac{4ig_{s}^{2}}{(k^{2}+i\epsilon)^{3}}\
\Big{(}-2k^{2}(\hat{O}^{LR})_{\alpha}^{\alpha}+2k^{\alpha}k^{\beta}\hat{O}^{LR}_{\alpha\beta}\Big{)}\to\frac{-6ig_{s}^{2}}{(k^{2}+i\epsilon)^{2}}(\hat{O}^{LR})_{\alpha}^{\alpha}\,,$
(B.4)

where the respective last step holds under symmetric integration. Finally,
using Fierz identities we can express $(\hat{O}^{LL})_{\alpha}^{\alpha}$ and
$(\hat{O}^{LR})_{\alpha}^{\alpha}$ in terms of the four-quark scalar operator
basis used in Ref. [42], obtaining

$\displaystyle g^{\mu\nu}\,\hat{\Pi}_{\mu\nu}^{LL}(k,0)$
$\displaystyle=\frac{ig_{s}^{2}}{(k^{2}+i\epsilon)^{2}}\ 2O_{1}\,,$
$\displaystyle g^{\mu\nu}\,\hat{\Pi}_{\mu\nu}^{LR}(k,0)$
$\displaystyle=\frac{ig_{s}^{2}}{(k^{2}+i\epsilon)^{2}}\
2\bar{O}_{LR}\,,\qquad\bar{O}_{LR}=\frac{1}{2}\left(O_{4}-3O_{5}\right)\,.$
(B.5)

These results lead to:

$\hat{O}^{LL}_{>}({\mathbf{k}})=\frac{3g_{s}^{2}}{4}\frac{1}{|{\mathbf{k}}|^{5}}\
O_{1}\,,\qquad\hat{O}^{LR}_{>}({\mathbf{k}})=\frac{3g_{s}^{2}}{4}\frac{1}{|{\mathbf{k}}|^{5}}\
\frac{1}{2}\left(O_{4}-3O_{5}\right)\,.$ (B.6)

As discussed in Ref. [42], the four-quark local operators $O_{1,4.5}$ can be
treated at low energy in chiral EFT. The discussion of $O_{1}$ is given in
Section 4.2 and we focus here on $O_{4,5}$. At LO only non-derivative four-
nucleon and pion–pion vertices appear,

$O_{4,5}\ \to\ g_{4,5}^{NN}\,\bar{p}n\bar{p}n\ +\
\frac{1}{2}g_{4,5}^{\pi\pi}\,F_{\pi}^{2}\,\pi^{-}\pi^{-}\,+\ldots$ (B.7)

The $\pi\pi$ couplings induce a pion-range transition operator, while the
$N\\!N$ couplings give a short-range transition operator. The leading chiral
expression for the $nn\to pp$ matrix element of $\bar{O}_{LR}$ is given by

$\langle f_{-}|\,\bar{O}_{LR}\,|i_{+}\rangle=\hat{\cal
A}_{A}+\chi^{+}_{{\mathbf{p}}^{\prime}}({\mathbf{0}})\,K_{E^{\prime}}\,\hat{\cal
A}_{B}+\ \hat{\cal
A}_{\bar{B}}\,K_{E}\,\chi^{+}_{{\mathbf{p}}}({\mathbf{0}})+\chi^{+}_{{\mathbf{p}}^{\prime}}({\mathbf{0}})\,K_{E^{\prime}}\,\left(\hat{\cal
A}_{C}+\frac{2g_{LR}^{NN}}{C^{2}}\right)\,K_{E}\,\chi^{+}_{{\mathbf{p}}}({\mathbf{0}})\,,$
(B.8)

where $\hat{\cal A}_{A,B,C}$ have the same formal expression as ${\cal
A}_{A,B,C}$ [43] with the replacement [42]

$V^{{}^{1}S_{0}}_{\nu\,{\rm L}}({\mathbf{q}})\ \to\
V^{{}^{1}S_{0}}_{\bar{O}_{LR}}({\mathbf{q}})\
=\frac{1}{2}\,g_{A}^{2}\,g_{LR}^{\pi\pi}\
\frac{{\mathbf{q}}^{2}}{({\mathbf{q}}^{2}+M_{\pi}^{2})^{2}}\,,$ (B.9)

where [42]

$g_{LR}^{NN}=\frac{1}{2}\left(g_{4}^{NN}-3g_{5}^{NN}\right)\sim\mathcal{O}((4\pi)^{2})\,,\qquad
g_{LR}^{\pi\pi}=\frac{1}{2}\left(g_{4}^{\pi\pi}-3g_{5}^{\pi\pi}\right)\sim\mathcal{O}((4\pi
F_{\pi})^{2})\,.$ (B.10)

The scaling of the pion couplings follows from naive dimensional analysis and
is confirmed by a lattice QCD calculation [90]. On the other hand, the scaling
of the $N\\!N$ couplings follows from the RGE in chiral EFT [42], which we
report here for completeness in terms of the rescaled couplings
$\tilde{g}_{4,5}^{NN}$ defined by
$g_{4,5}^{NN}=(m_{N}C/(4\pi))^{2}\,\tilde{g}_{4,5}^{NN}$:

$\frac{d\,\tilde{g}_{4,5}^{NN}}{d\log\nu_{\chi}}=\frac{g_{A}^{2}}{4}\,g_{4,5}^{\pi\pi}\,.$
(B.11)

The scaling for the nucleon couplings follows from the above RGE and the
scaling $m_{N}C/(4\pi)\sim 1/Q$. As a consequence of the above discussion, the
coupling $g_{LR}^{NN}$ depends on both the short-distance QCD renormalization
scale $\mu$ and on the chiral renormalization scale $\nu_{\chi}$.

Using the above relations we arrive at

${\cal
M}_{LR}^{>}=\mathcal{A}^{>}_{A}+\chi^{+}_{{\mathbf{p}}^{\prime}}({\mathbf{0}})\,K_{E^{\prime}}\,\mathcal{A}^{>}_{B}+{\mathcal{A}}^{>}_{\bar{B}}\,K_{E}\,\chi^{+}_{{\mathbf{p}}}({\mathbf{0}})+\chi^{+}_{{\mathbf{p}}^{\prime}}({\mathbf{0}})\,K_{E^{\prime}}\,\mathcal{A}^{>}_{C}\,K_{E}\,\chi^{+}_{{\mathbf{p}}}({\mathbf{0}})\,,$
(B.12)

with

$\displaystyle{\cal A}_{A,B,\bar{B}}^{>}$
$\displaystyle=\frac{3\alpha_{s}}{2\pi}\,\int_{\Lambda}^{\infty}\
d|{\mathbf{k}}|\ \frac{1}{|{\mathbf{k}}|^{3}}\ \hat{\cal A}_{A,B,\bar{B}}\,,$
$\displaystyle{\cal A}_{C}^{>}$
$\displaystyle=\frac{3\alpha_{s}}{2\pi}\,\int_{\Lambda}^{\infty}\
d|{\mathbf{k}}|\ \frac{1}{|{\mathbf{k}}|^{3}}\
\left(\frac{2g_{LR}^{NN}}{C^{2}}+\hat{\cal A}_{C}\right)\,.$ (B.13)

For matching purposes only ${\cal A}_{C}^{>}$ matters. In the numerical
estimates we will retain only the term proportional to $g_{LR}^{NN}$ and the
singular part of $\hat{\cal A}_{C}$, proportional to $g_{LR}^{\pi\pi}$.
Further, for convenience in the matching analysis we express the couplings
$g_{LR}^{NN}$ and $g_{LR}^{\pi\pi}$ in terms of dimensionless quantities as
follows

$\displaystyle g_{LR}^{NN}$
$\displaystyle=\left(\frac{m_{N}}{4\pi}C\right)^{2}\,(4\pi
F_{\pi})^{2}\,\bar{g}_{LR}^{NN}\,,$ $\displaystyle\bar{g}_{LR}^{NN}$
$\displaystyle\sim\mathcal{O}(1)\,,$ $\displaystyle g_{LR}^{\pi\pi}$
$\displaystyle=(4\pi F_{\pi})^{2}\,\bar{g}_{LR}^{\pi\pi}\,,$
$\displaystyle\bar{g}_{LR}^{\pi\pi}$ $\displaystyle\sim\mathcal{O}(1)\,.$
(B.14)

In numerical estimates we will use $\bar{g}_{LR}^{\pi\pi}=8.23$ (in the
$\overline{\rm MS}$ scheme at $\mu=2\,\text{GeV}$, based on the lattice QCD
results of Ref. [90]) and assume the range $\bar{g}_{LR}^{NN}\in[-10,+10]$.

##  Appendix C Estimating inelastic effects

As a representative of inelastic effects, we estimate the contribution from
the diagram in Fig. 8.171717This corresponds to diagram $(G)$ in Fig. 3 of
Ref. [38]. The contribution to the generalized Compton amplitude in Eq. (2.9)
is already a two-loop diagram, which with appropriate momentum routing can be
written as

$\displaystyle\delta{\cal T}(k,p_{\rm ext})$ $\displaystyle=iC^{2}\
\frac{(1+3g_{A}^{2})m_{N}^{2}}{16F_{\pi}^{2}}$
$\displaystyle\times\int\frac{d^{3}{\mathbf{q}}}{(2\pi)^{3}}\int\frac{d^{3}{\mathbf{q}}^{\prime}}{(2\pi)^{3}}\,\frac{1}{{\mathbf{p}}^{2}-{\mathbf{q}}^{2}+i\epsilon}\,\frac{1}{{\mathbf{p}}^{\prime
2}-{\mathbf{q}}^{\prime
2}+i\epsilon}\frac{1}{-(k^{0})^{2}+({\mathbf{q}}-{\mathbf{q}}^{\prime}+{\mathbf{k}})^{2}+M_{\pi}^{2}}\,.$
(C.1)

Upon changing variables to ${\mathbf{q}}$ and
${\mathbf{l}}={\mathbf{q}}-{\mathbf{q}}^{\prime}$ one has

$\delta{\cal T}(k,p_{\rm ext})=iC^{2}\
\frac{(1+3g_{A}^{2})m_{N}^{2}}{16F_{\pi}^{2}}\
\int\frac{d^{3}{\mathbf{l}}}{(2\pi)^{3}}\
\frac{1}{-(k^{0})^{2}+({\mathbf{l}}+{\mathbf{k}})^{2}+M_{\pi}^{2}-i\epsilon}\
{\cal I}_{C}({\mathbf{l}}^{2},{\mathbf{p}}^{2},{\mathbf{p}}^{\prime 2})\,.$
(C.2)

Recalling that ${\cal I}_{C}\to 1/(8|{\mathbf{l}}|)$ for
${\mathbf{p}}={\mathbf{p}}^{\prime}\to 0$, one obtains

$\delta{\cal T}(k,p_{\rm ext})=iC^{2}\
\frac{(1+3g_{A}^{2})m_{N}^{2}}{128F_{\pi}^{2}}\
\int\frac{d^{3}{\mathbf{l}}}{(2\pi)^{3}}\
\frac{1}{-(k^{0})^{2}+({\mathbf{l}}+{\mathbf{k}})^{2}+M_{\pi}^{2}-i\epsilon}\
\frac{1}{|{\mathbf{l}}|}\,.$ (C.3)

The integral is logarithmically divergent and we regulate the divergence with
a cutoff scale $\nu\sim\Lambda_{\chi}$. We will mimic the effect of the
associated (and unknown) counter term by varying the scale $\nu$ in order to
obtain a rough estimate of the corresponding inelastic effect.

The inelastic contribution to ${\cal A}_{C}^{<}$ is given by

$\delta{\cal A}_{C}^{<}=2\,\int\frac{d^{4}k}{(2\pi)^{4}}\frac{\delta{\cal
T}(k,p_{\rm ext})}{k^{2}+i\epsilon}\,.$ (C.4)

We reduce this integral by first integrating over $k^{0}$ via the residues at
the pion and neutrino poles, and second performing the angular integral in
$d^{3}{\mathbf{l}}$, obtaining

$\delta{\cal A}_{C}^{<}=\frac{m_{N}^{2}}{(4\pi)^{2}}\frac{1}{(4\pi
F_{\pi})^{2}}\frac{1}{4}\ \int\,d|{\mathbf{k}}|\,(1+3g_{A}^{2})\
\int_{0}^{\nu}d|{\mathbf{l}}|\
\log\left[\frac{|{\mathbf{k}}|+\sqrt{(|{\mathbf{k}}|+|{\mathbf{l}}|)^{2}+M_{\pi}^{2}}}{|{\mathbf{k}}|+\sqrt{(|{\mathbf{k}}|-|{\mathbf{l}}|)^{2}+M_{\pi}^{2}}}\right]\,.$
(C.5)

The integral in $d|{\mathbf{l}}|$ can be done analytically, leading to

$\displaystyle\delta{\cal A}_{C}^{<}$
$\displaystyle=\frac{m_{N}^{2}}{(4\pi)^{2}}\ \int\,d|{\mathbf{k}}|\ \delta
a_{<}(|{\mathbf{k}}|)\,,$ $\displaystyle\delta a_{<}(|{\mathbf{k}}|)$
$\displaystyle=\frac{1+3g_{A}^{2}}{2(4\pi F_{\pi})^{2}}\
\left[|{\mathbf{k}}|\left(1+\log\frac{\nu\,M_{\pi}}{(|{\mathbf{k}}|+\sqrt{{\mathbf{k}}^{2}+M_{\pi}^{2}})^{2}}\right)+M_{\pi}F\bigg{(}\frac{|{\mathbf{k}}|}{M_{\pi}}\bigg{)}\right]\,,$
$\displaystyle F(x)$
$\displaystyle=\sqrt{1-x^{2}}\,\left[\arctan\left(\frac{x^{2}}{\sqrt{1-x^{4}}}\right)-\arctan\left(\frac{x}{\sqrt{1-x^{2}}}\right)\right]\,.$
(C.6)

In the $M_{\pi}\to 0$ limit the integrand becomes

$\delta a_{<}(|{\mathbf{k}}|)=\frac{1+3g_{A}^{2}}{2(4\pi
F_{\pi})^{2}}\,|{\mathbf{k}}|\left(1+\log\frac{\nu}{4|{\mathbf{k}}|}\right)\,.$
(C.7)

The logarithmic divergence signals that there must be, at this order, a
counter term for the generalized Compton amplitude (the divergence arose at
that level, before integrating over $d^{3}{\mathbf{k}}$). Naive dimensional
analysis suggests that this counter term should be of the same order of the
coefficient of $\log\nu$. To obtain a numerical estimate we vary the scale
within $\nu\sim(0.5\text{--}1.0)\,\text{GeV}$. The new contribution to the
integrand is clearly suppressed by two chiral orders compared to the leading
terms, as long as $|{\mathbf{k}}|\ll\Lambda_{\chi}$, as one is comparing the
LO $1/|{\mathbf{k}}|$ to $|{\mathbf{k}}|/(4\pi F_{\pi})^{2}$. However, this
leads to a badly divergent integral in $|{\mathbf{k}}|$. In order to estimate
the contribution to the integral, and therefore to the LEC,

$\delta\tilde{C}_{1}=\frac{1}{2}\ \delta\bar{\cal A}_{C}^{<}\,,$ (C.8)

we explore two options: first, one may use the chiral EFT form of the
integrand and cut it off at $\Lambda_{\chi}\sim 1\,\text{GeV}$, or, second,
one may replace $1+3g_{A}^{2}\to
g_{V}^{2}({\mathbf{k}}^{2})+3g_{A}^{2}({\mathbf{k}}^{2})$ in Eq. (C.6) and
integrate up to $\Lambda$. Numerically, the second option leads to
$|\delta\tilde{\cal C}_{1}|=0.15$, while the first option leads to
$|\delta\tilde{\cal C}_{1}|=0.35$ (the largest absolute value is obtained for
$\nu=0.5\,\text{GeV}$). Based on these observations, we will take into account
the effects of inelastic channels by adding a $\pm 0.5$ uncertainty to our
estimate of $\tilde{\cal C}_{1}(\mu_{\chi})$.

##  Appendix D Evaluation of $\boldsymbol{Z^{<}}$

Using ${\cal T}_{<}^{\pi\pi}(k,0)$ from Eqs. (6.15), (6.16), and (6.17) we
obtain

$Z^{<}=\frac{3i}{2F^{2}}\,\int\frac{d^{4}k}{(2\pi)^{4}}\
\frac{1}{k^{2}+i\epsilon}\left(\frac{M_{V}^{2}}{M_{V}^{2}-k^{2}}\right)^{2}\,.$
(D.1)

We evaluate this integral in two ways. First, performing the usual Wick
rotation and introducing the variable $k_{E}^{2}=(k^{0})^{2}+{\mathbf{k}}^{2}$
to reduce the integral to181818One could also perform the integrals in analogy
to the derivation of the Cottingham formula and would obtain the same result
in this case, because after the change of variables the integrand does not
depend on $k^{0}$ and
$\int_{-k_{E}}^{+k_{E}}dk^{0}\sqrt{k_{E}^{2}-(k^{0})^{2}}=(\pi k_{E}^{2})/2$.

$Z^{<}=\frac{3}{2}\frac{1}{(4\pi
F)^{2}}\,\int_{0}^{\infty}dk_{E}^{2}\left(\frac{M_{V}^{2}}{M_{V}^{2}+k_{E}^{2}}\right)^{2}=\frac{3}{2}\frac{M_{V}^{2}}{(4\pi
F)^{2}}\,\int_{0}^{\infty}dx\left(\frac{1}{1+x}\right)^{2}=\frac{3}{2}\frac{M_{V}^{2}}{(4\pi
F)^{2}}\,.$ (D.2)

Since the integral is convergent, we have taken the integration limit in
$k_{E}^{2}$ to $+\infty$. Note that the overall normalization of the integral
above agrees with the literature on the pion mass difference, see, e.g., Ref.
[63].

Next, we perform the integral (D.1) by first carrying out the $k^{0}$
integration via the residues’ theorem, which allows us to identify the
integrand in the variable $|{\mathbf{k}}|$ in Eq. (6.21). The integrand has
single poles at $k^{0}=\pm(|{\mathbf{k}}|-i\epsilon)$ and double poles at
$k^{0}=\pm(\omega_{V}-i\epsilon)$, where
$\omega_{V}=\sqrt{{\mathbf{k}}^{2}+M_{V}^{2}}$. The integration gives:

$\displaystyle Z^{<}$
$\displaystyle=\frac{3}{4F^{2}}\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\
\frac{1}{|{\mathbf{k}}|}\left(1-\frac{|{\mathbf{k}}|}{\omega_{V}}\left(1+\frac{M_{V}^{2}}{2\omega_{V}^{2}}\right)\right)$
$\displaystyle=\frac{3}{2}\frac{1}{(4\pi F)^{2}}\
4\int_{0}^{\Lambda}d|{\mathbf{k}}||{\mathbf{k}}|\left(1-\frac{|{\mathbf{k}}|}{\omega_{V}}\left(1+\frac{M_{V}^{2}}{2\omega_{V}^{2}}\right)\right)\,,$
(D.3)

which is the result used in Eq. (6.21). As a check, we perform the integration
up to $\Lambda\to\infty$, obtaining

$Z^{<}(\Lambda\to\infty)=\frac{3}{2}\frac{M_{V}^{2}}{(4\pi F)^{2}}\
4\int_{0}^{\infty}dxx\left(1-\frac{x}{2}\frac{3+2x^{2}}{(1+x^{2})^{3/2}}\right)=\frac{3}{2}\frac{M_{V}^{2}}{(4\pi
F)^{2}}\,.$ (D.4)

In the last step we used the fact that the integral in the $x$ variable
evaluates to $1/4$. This result is consistent with Eq. (D.2), obtained with
the more commonly used integration variable $k_{E}^{2}$. The advantage of the
standard approach is that the result directly applies to an arbitrary pion
vector form factor, while collecting the residues in the $k^{0}$ integration
assumes the form dictated by vector meson dominance. For our application this
approximation is completely sufficient, but could be improved by using a
dispersive representation of $F_{\pi}^{V}(k^{2})$ [119, 120] and treating the
Cauchy kernel along the same lines.

##  Appendix E Details on CIB in $\boldsymbol{N\\!N}$ scattering near
threshold

###  E.1 Low-energy $\boldsymbol{N\\!N}$ scattering in the presence of Coulomb
interactions

We begin by recalling some results on the Coulomb-modified effective range
expansion for $pp$ scattering [97, 98, 99]. We specialize to the $S$-wave and
indicate the pure Coulomb phase shift by $\sigma$ and the full phase shift by
$\sigma+\nu$. Further, we denote the total scattering amplitude by $T$, the
purely Coulomb component by $T_{C}$, and the Coulomb-modified strong amplitude
by $T_{SC}$

$\displaystyle T$
$\displaystyle=T_{C}+T_{SC}=\frac{4\pi}{m_{p}}\frac{e^{2i(\sigma(k)+\nu(k))}-1}{2ik}\,,$
$\displaystyle T_{C}$
$\displaystyle=\frac{4\pi}{m_{p}}\frac{e^{2i\sigma(k)}-1}{2ik}\,,$
$\displaystyle T_{SC}$
$\displaystyle=\frac{4\pi}{m_{p}}e^{i2\sigma(k)}\frac{\left(e^{2i\nu(k)}-1\right)}{2ik}=\frac{4\pi}{m_{p}}e^{i2\sigma(k)}\frac{1}{k\left({\rm
cot}(\nu(k))-i\right)}\,.$ (E.1)

The phase shift $\nu(k)$ obeys a modified effective range expansion [97, 98,
99]

$\frac{1}{k\left({\rm
cot}(\nu(k))-i\right)}=\frac{C_{0}^{2}(\eta)}{-\frac{1}{a_{C}}-r_{0}\frac{k^{2}}{2}-\alpha
m_{p}h(\eta)+\ldots}\,,$ (E.2)

where

$\eta=\frac{\alpha m_{p}}{2k}\,,\qquad
C_{0}^{2}(\eta)=\frac{2\pi\eta}{e^{2\pi\eta}-1}\,,\qquad h(\eta)={\rm
Re}\left(\frac{\Gamma^{\prime}(i\eta)}{\Gamma(i\eta)}\right)-\log\eta\,.$
(E.3)

The standard effective range expansion is recovered formally by setting
$\alpha\to 0$, which also implies $C_{0}^{2}(\eta)\to 1$. For $\alpha\neq 0$,
the function $h(\eta)$ goes smoothly to zero for $k\to 0$, so that $a_{C}$ can
be interpreted as a scattering length. However, the Sommerfeld factor
$C_{0}^{2}(\eta)$—the square of the Coulomb wave function at the origin—goes
to zero for $k\to 0$ due to Coulomb repulsion and this prevents the usual
threshold analysis of the amplitude. It will prove useful to work with the
momentum-dependent quantity

$\tilde{a}_{pp}(k)=\frac{C_{0}^{2}(\eta)}{\frac{1}{a_{C}}+\alpha
m_{p}h(\eta)}\,,$ (E.4)

which is a proxy for the very low-energy Coulomb-subtracted $pp$ scattering
amplitude.

For the CIB analysis we will consider the quantities

$\tilde{a}_{\rm CIB}(k)=\frac{a_{nn}+\tilde{a}_{pp}(k)}{2}-a_{np}$ (E.5)

as well as

${a}_{\rm CIB}=\frac{a_{nn}+{a}^{C}_{pp}}{2}-a_{np}\,,$ (E.6)

and we will use the following values for the scattering lengths determined
from $N\\!N$ data [103, 104, 105, 101, 102]:

$a_{np}=-23.74(2)\,\textrm{fm}\,,\qquad a_{nn}=-18.9(4)\,\textrm{fm}\,,\qquad
a_{pp}^{C}=-7.817(4)\,\textrm{fm}\,.$ (E.7)

In Fig. 18 we show $\tilde{a}_{\rm CIB}(k)$ using the input from Eq. (E.7).

Figure 18:  $\tilde{a}_{\rm CIB}(k)$ defined in Eq. (E.5) versus $k$, using
experimental input from Eq. (E.7).

###  E.2 Low-energy $\boldsymbol{N\\!N}$ amplitudes in chiral EFT with isospin
breaking

Using the techniques of Ref. [23], we extract the phase shift $\nu(k)$ from
the asymptotic behavior of the solution to the Schrödinger equation with
inclusion of Coulomb potential and the LO chiral potential. The ${}^{1}S_{0}$
chiral potential includes the LO channel-dependent contact interaction
(denoted by $C_{nn,pp,np}$ in the following), as well as the channel-dependent
Yukawa potential with the electromagnetically induced pion mass splitting. By
including these effects in the Schrödinger equation we effectively iterate the
isospin-breaking corrections to all orders. This is mandatory for the Coulomb
potential, while for the other terms a perturbative treatment should be
numerically close to our analysis, see Ref. [43] for additional details.

The phase shift $\nu(k)$ is related to the complex coefficients $y_{pp}(k)$
and $z_{pp}(k)$ controlling the large-$r$ behavior of the regular and
irregular solutions of the Schrödinger equation with Coulomb and Yukawa
potentials

$\displaystyle\psi_{\text{reg}}(r)$ $\displaystyle\to y_{pp}(k)\,\frac{\
e^{i(kr-\eta\log(2kr)+\sigma(k))}}{kr}+{\rm c.c.}\,,$
$\displaystyle\psi_{\text{irr}}(r)$ $\displaystyle\to z_{pp}(k)\,\frac{\
e^{i(kr-\eta\log(2kr)+\sigma(k))}}{kr}+{\rm c.c.}$ (E.8)

By matching the low-energy form of the amplitude obtained in EFT to the
(modified) effective range form, one obtains

$\displaystyle a_{np}$ $\displaystyle=a_{np}^{\pi}+G(C_{np},y_{np},z_{np})\,,$
$\displaystyle a_{nn}$ $\displaystyle=a_{nn}^{\pi}+G(C_{nn},y_{nn},z_{nn})\,,$
$\displaystyle\tilde{a}_{pp}(k)$
$\displaystyle=a_{pp}^{\pi}+G(C_{pp},y_{pp},z_{pp})\,,$ $\displaystyle
G(C_{ij},y_{ij},z_{ij})$ $\displaystyle=\frac{m_{ij}}{16\pi y_{ij}^{*2}}\
\frac{1}{-\frac{1}{C_{ij}(\mu)}+\Delta_{ij}\big{(}\frac{\mu}{\lambda}\big{)}+\frac{z^{*}(\lambda)}{y^{*}}}\,.$
(E.9)

In the above relations $a_{ij}^{\pi}$ is the Yukawa-induced scattering length
(obtained in terms of the phase of $y_{ij}(k)$ at small $k$). $m_{ij}$ is the
appropriate channel-dependent mass for the nucleon system ($m_{nn}=m_{n}$,
$m_{pp}=m_{p}$, $m_{np}=2m_{n}m_{p}/(m_{n}+m_{p})$). The quantities $z_{ij}$
depend on the regulator $\lambda$ introduced when imposing boundary conditions
for the irregular solution of the Schrödinger equation [23]. The results
quoted below are obtained with $1/\lambda=0.001\,\text{fm}$ and we have
checked stability in the range $1/\lambda=0.001\to 0.05\,\text{fm}$. The
dependence on $\lambda$ is canceled by the matching factors
$\Delta_{ij}(\mu/\lambda)$, to connect to the $\overline{\rm MS}$ scheme, in
which the contact couplings $C_{ij}(\mu)$ are defined. The explicit form of
these scheme-changing factors is

$\displaystyle\Delta_{nn}$
$\displaystyle=-\alpha_{\pi}\frac{m_{nn}^{2}}{8\pi}\left(\log\frac{\mu^{2}}{\lambda^{2}}+2\gamma_{E}-1\right)\,,\qquad\alpha_{\pi}=\frac{g_{A}^{2}M_{\pi^{0}}^{2}}{16\pi
F_{\pi}^{2}}\,,$ $\displaystyle\Delta_{np}$
$\displaystyle=-\alpha_{\pi}\left(1+2\frac{M_{\pi^{\pm}}^{2}-M_{\pi^{0}}^{2}}{M_{\pi^{0}}^{2}}\right)\frac{m_{np}^{2}}{8\pi}\left(\log\frac{\mu^{2}}{\lambda^{2}}+2\gamma_{E}-1\right)\,,$
$\displaystyle\Delta_{pp}$
$\displaystyle=-(\alpha_{\pi}-\alpha)\,\frac{m_{pp}^{2}}{8\pi}\left(\log\frac{\mu^{2}}{\lambda^{2}}+2\gamma_{E}-1\right)\,.$
(E.10)

In the presence of isospin breaking, the short-range ${}^{1}S_{0}$ $N\\!N$
couplings can be decomposed as follows

$\displaystyle C_{np}$ $\displaystyle=C+\frac{e^{2}}{3}\left({\cal
C}_{1}+{\cal C}_{2}\right)\,,$ $\displaystyle C_{nn}$
$\displaystyle=C-\frac{e^{2}}{6}\left({\cal C}_{1}+{\cal
C}_{2}\right)+\frac{1}{2}{\cal C}_{\rm CSB}\,,$ $\displaystyle C_{pp}$
$\displaystyle=C-\frac{e^{2}}{6}\left({\cal C}_{1}+{\cal
C}_{2}\right)-\frac{1}{2}{\cal C}_{\rm CSB}\,,$ (E.11)

where the CIB combination ${\cal C}_{1}+{\cal C}_{2}$ arises form $I=2$
interactions, while the CSB term ${\cal C}_{\rm CSB}$ stems from $I=1$
interactions. The latter term can originate from strong isospin-breaking
$\propto(m_{u}-m_{d})$ or electromagnetic interactions. For the CIB coupling,
renormalization arguments (cancellation of divergences) enforces the LO
scaling $\sim e^{2}/Q^{2}$ [100, 40, 43]. For the CSB coupling, the Coulomb
potential also induces a UV divergence that requires it to scale as $\sim
e^{2}/Q^{2}$.

###  E.3 Fitting the couplings

Equation (E.9) can be used to extract the three contact couplings
$C_{np,nn,pp}$ (or equivalently $C$, ${\cal C}_{1}+{\cal C}_{2}$, ${\cal
C}_{\rm CSB}$). The output is very stable when using momenta near threshold,
in the range $k\in[1\,\text{MeV},10\,\text{MeV}]$. For definiteness, we quote
results for $k=1\,\text{MeV}$. Using the reference scale
$\mu=M_{\pi}\equiv(M_{\pi^{0}}+2M_{\pi^{\pm}})/3$, this analysis leads to

$\displaystyle C_{ij}(\mu)$
$\displaystyle=-\frac{1}{(\Lambda_{ij}(\mu))^{2}}\,,$
$\displaystyle\Lambda_{nn,pp,np}(M_{\pi})$
$\displaystyle=(95.4,88.7,98.9)\,\text{MeV}\,,$ $\displaystyle C(\mu)$
$\displaystyle=-\frac{1}{(\Lambda_{C}(\mu))^{2}}\,,$
$\displaystyle\Lambda_{C}(M_{\pi})$ $\displaystyle=94.4\,\text{MeV}\,,$
$\displaystyle\frac{e^{2}}{2}({\cal C}_{1}+{\cal C}_{2})(\mu)$
$\displaystyle=\frac{1}{(\Lambda_{\rm CIB}(\mu))^{2}}\,,$
$\displaystyle\Lambda_{\rm CIB}(M_{\pi})$
$\displaystyle=247\,\text{MeV}\,,\qquad(\tilde{\cal C}_{1}+\tilde{\cal
C}_{2})(M_{\pi})=5.1\,,$ $\displaystyle{\cal C}_{\rm CSB}(\mu)$
$\displaystyle=\frac{1}{(\Lambda_{\rm CSB}(\mu))^{2}}\,,$
$\displaystyle\Lambda_{\rm CSB}(M_{\pi})$ $\displaystyle=254\,\text{MeV}\,.$
(E.12)

The result indicates that both CIB and CSB effects provide corrections to the
isosymmetric coupling $C$ at the (10–15)% level.

###  E.4 Predicting the CIB combination of scattering lengths

Alternatively, we can use Eq. (E.9) to validate our calculation of $({\cal
C}_{1}+{\cal C}_{2})(\mu)$ as follows:

  1. 1.

We use one linear combination of Eq. (E.9) to extract the isospin-conserving
coupling $C$.

Figure 19:  $\tilde{a}^{\rm th(1)}_{\rm CIB}$ (lower curve) and
$\tilde{a}^{\rm th(2)}_{\rm CIB}$ (upper curve), both defined in Eq. (E.13),
versus the dimensionless coupling $\tilde{\cal C}_{1}+\tilde{\cal C}_{2}$. The
curves are evaluated at $k_{0}=1\,\text{MeV}$ for $\mu=M_{\pi}$ (left panel)
and $\mu=2M_{\pi}$ (right panel). The horizontal line represents the
experimental value $a_{\rm CIB}(k=1\,\text{MeV})=14.25\,\text{fm}$. For
reference, our theoretical ranges for the CIB LEC are $(\tilde{\cal
C}_{1}+\tilde{\cal C}_{2})(M_{\pi})=2.9(1.2)$ and $(\tilde{\cal
C}_{1}+\tilde{\cal C}_{2})(2M_{\pi})=5.4(1.2)$.

  2. 2.

We construct the CIB combinations $\tilde{a}_{\rm CIB}(k)$ and $a_{\rm CIB}$.
The theoretical expressions constructed from the right-hand side of Eq. (E.9)
depend in general on $C$, $({\cal C}_{1}+{\cal C}_{2})$, and ${\cal C}_{\rm
CSB}$. We will fix $C$ to the fit value (see previous step) and consider the
dependence of the theoretical expressions $\tilde{a}^{\rm th}_{\rm CIB}(k)$
and $a_{\rm CIB}^{\rm th}$ on the $I=2$ coupling $({\cal C}_{1}+{\cal
C}_{2})(\mu)$. Using our range for $({\cal C}_{1}+{\cal C}_{2})(\mu)$ we can
predict $\tilde{a}_{\rm CIB}(k)$ and $a_{\rm CIB}$, and compare to the
experimental value.

  3. 3.

Working to first order in isospin-breaking quantities, the observable
$\tilde{a}^{\rm th}_{\rm CIB}(k)$ should be insensitive to the $I=1$ coupling
${\cal C}_{\rm CSB}$. However, in Eq. (E.9) we are effectively iterating the
insertions of electromagnetic isospin-breaking sources to all orders
(including ${\cal C}_{\rm CSB}$), since we are iterating the full couplings
$C_{ij}$. As indicated by ${\cal C}_{\rm CSB}/C\sim 15\%$, second-order
effects due to two insertions of ${\cal C}_{\rm CSB}$ might not be negligible.
To quantify the impact of this on our CIB analysis, in Fig. 19 we show for
$k_{0}=1\,\text{MeV}$, $\mu=M_{\pi}$ (left panel), and $\mu=2M_{\pi}$ (right
panel) the two functions

$\displaystyle\tilde{a}^{\rm th(1)}_{\rm CIB}(\tilde{\cal C}_{1}+\tilde{\cal
C}_{2})$ $\displaystyle=\tilde{a}^{\rm th}_{\rm CIB}(k=k_{0})\Big{|}_{C=C_{\rm
fit},\,{\cal C}_{\rm CSB}=0}\,,$ $\displaystyle\tilde{a}^{\rm th(2)}_{\rm
CIB}(\tilde{\cal C}_{1}+\tilde{\cal C}_{2})$ $\displaystyle=\tilde{a}^{\rm
th}_{\rm CIB}(k=k_{0})\Big{|}_{C=C_{\rm fit},\,{\cal C}_{\rm CSB}={\cal
C}_{\rm CSB,\,fit}}\,,$ (E.13)

where in both diagrams the horizontal line represents the experimental value
of $\tilde{a}_{\rm CIB}(k=1\,\text{MeV})=14.25\,\text{fm}$.

  4. 4.

A few comments are in order: (i) The plots clearly show the need for a LO CIB
coupling. Not including $(\tilde{\cal C}_{1}+\tilde{\cal C}_{2})$ would lead
to a scale-dependent prediction for $\tilde{a}_{\rm CIB}(k)$. Moreover, at
$\mu=M_{\pi}$ the theory would predict $\tilde{a}_{\rm
CIB}(k=1\,\text{MeV})\sim 30\,\text{fm}$, more than a factor of two larger
than the experimental value. (ii) The two choices of treating ${\cal C}_{\rm
CSB}$ lead to a shift $\Delta\tilde{a}_{\rm CIB}\sim 2.5\,\text{fm}$ at
$\mu=M_{\pi}$ (growing to $\sim 4\,\text{fm}$ at $\mu=2M_{\pi}$), subdominant
but not negligible. (iii) One could use the intercept of the curves in Fig. 19
with the horizontal experimental constraint to perform an alternative
extraction of $(\tilde{\cal C}_{1}+\tilde{\cal C}_{2})$. Due to the missing
CSB contribution, using $\tilde{a}^{\rm th(1)}_{\rm CIB}$ implies a value of
$(\tilde{\cal C}_{1}+\tilde{\cal C}_{2})(M_{\pi})\sim 3.7$, to be contrasted
with the $5.1$ from the previous analysis. (iv) Finally, we note that using
$\tilde{a}^{\rm th(2)}_{\rm CIB}$ produces results that are more stable under
variation of $\mu$ (see below).

  5. 5.

Our theoretical ranges for the CIB LEC are $(\tilde{\cal C}_{1}+\tilde{\cal
C}_{2})(M_{\pi})=2.9(1.2)$ and $(\tilde{\cal C}_{1}+\tilde{\cal
C}_{2})(2M_{\pi})=5.4(1.2)$. These values can be used in conjunction with
$a^{\rm th(1,2)}_{\rm CIB}(\tilde{\cal C}_{1}+\tilde{\cal C}_{2})$ to predict
$a_{\rm CIB}$. Using $\tilde{a}^{\rm th(1)}_{\rm CIB}$ we find

$\displaystyle\tilde{a}_{\rm CIB}(k=1\,\text{MeV})$
$\displaystyle=17^{+5}_{-4}\,\text{fm}$ $\displaystyle\mu$
$\displaystyle=M_{\pi}\,,$ $\displaystyle\tilde{a}_{\rm CIB}(k=1\,\text{MeV})$
$\displaystyle=15.3^{+5}_{-4}\,\text{fm}$ $\displaystyle\mu$
$\displaystyle=2M_{\pi}\,,$ (E.14)

while the second variant $\tilde{a}^{\rm th(2)}_{\rm CIB}$ gives

$\displaystyle\tilde{a}_{\rm CIB}(k=1\,\text{MeV})$
$\displaystyle=19.6^{+4.5}_{-4}\,\text{fm}$ $\displaystyle\mu$
$\displaystyle=M_{\pi}\,,$ $\displaystyle\tilde{a}_{\rm CIB}(k=1\,\text{MeV})$
$\displaystyle=19.0^{+4.5}_{-4}\,\text{fm}$ $\displaystyle\mu$
$\displaystyle=2M_{\pi}\,.$ (E.15)

These values compare well with $\tilde{a}^{\rm exp}_{\rm
CIB}(k=1\,\text{MeV})=14.25\,\text{fm}$. Using $k=10\,\text{MeV}$, the
experimental value decreases to $13.0\,\text{fm}$ and the prediction is also
reduced by slightly over $1\,\text{fm}$.

  6. 6.

We can repeat the same analysis using the momentum-independent combination of
scattering lengths

$a_{\rm CIB}=\frac{a_{nn}+a_{pp}^{C}}{2}-a_{np}\approx 10.4\,\text{fm}\,,$
(E.16)

by constructing theoretical quantities in complete analogy to Eq. (E.13). The
results are illustrated in Fig. 20. Using $a^{\rm th(1)}_{\rm CIB}$ we find

$\displaystyle a_{\rm CIB}$ $\displaystyle=13.5^{+5}_{-4}\,\text{fm}$
$\displaystyle\mu$ $\displaystyle=M_{\pi}\,,$ $\displaystyle a_{\rm CIB}$
$\displaystyle=12^{+5}_{-4}\,\text{fm}$ $\displaystyle\mu$
$\displaystyle=2M_{\pi}\,,$ (E.17)

and for $a^{\rm th(2)}_{\rm CIB}$

$\displaystyle a_{\rm CIB}$ $\displaystyle=16.0^{+4.5}_{-4}\,\text{fm}$
$\displaystyle\mu$ $\displaystyle=M_{\pi}\,,$ $\displaystyle a_{\rm CIB}$
$\displaystyle=15.4^{+4.5}_{-4}\,\text{fm}$ $\displaystyle\mu$
$\displaystyle=2M_{\pi}\,.$ (E.18)

Figure 20:  $a^{\rm th(1)}_{\rm CIB}$ (lower curve) and $a^{\rm th(2)}_{\rm
CIB}$ (upper curve) versus the dimensionless coupling $\tilde{\cal
C}_{1}+\tilde{\cal C}_{2}$. The curves are evaluated at $\mu=M_{\pi}$ (left
panel) and $\mu=2M_{\pi}$ (right panel). The horizontal line represents the
experimental value $a_{\rm CIB}=10.35\,\text{fm}$. For reference, our
theoretical ranges for the CIB LEC are $(\tilde{\cal C}_{1}+\tilde{\cal
C}_{2})(M_{\pi})=2.9(1.2)$ and $(\tilde{\cal C}_{1}+\tilde{\cal
C}_{2})(2M_{\pi})=5.4(1.2)$.

In summary, we have derived a prediction of the CIB scattering length that
overshoots the experimental value by about $35\%$ and has a comparable
uncertainty, thus comparing well with data. This is a significant
phenomenological success of our theoretical approach and supports the validity
of our uncertainty estimates.

## References

  * [1] W. Furry,  Phys. Rev. 56, 1184 (1939).
  * [2] J. Schechter and J. Valle,  Phys. Rev. D 25, 2951 (1982).
  * [3] S. Davidson, E. Nardi, and Y. Nir,  Phys. Rept. 466, 105 (2008), 0802.2962.
  * [4] W. Rodejohann,  Int. J. Mod. Phys. E 20, 1833 (2011), 1106.1334.
  * [5] KamLAND-Zen, A. Gando et al.,  Phys. Rev. Lett. 117, 082503 (2016), 1605.02889,  [Addendum: Phys. Rev. Lett. 117, 109903 (2016)].
  * [6] NEMO-3, R. Arnold et al.,  Phys. Rev. D 94, 072003 (2016), 1606.08494.
  * [7] EXO, J. Albert et al.,  Phys. Rev. Lett. 120, 072701 (2018), 1707.08707.
  * [8] Majorana, C. Aalseth et al.,  Phys. Rev. Lett. 120, 132502 (2018), 1710.11608.
  * [9] CUORE, D. Adams et al.,  Phys. Rev. Lett. 124, 122501 (2020), 1912.10966.
  * [10] GERDA, M. Agostini et al.,  Phys. Rev. Lett. 125, 252502 (2020), 2009.06079.
  * [11] I. Avignone, Frank T., S. R. Elliott, and J. Engel,  Rev. Mod. Phys. 80, 481 (2008), 0708.1033.
  * [12] J. Menéndez, A. Poves, E. Caurier, and F. Nowacki,  Nucl. Phys. A 818, 139 (2009), 0801.3760.
  * [13] J. Vergados, H. Ejiri, and F. Šimkovic,  Rept. Prog. Phys. 75, 106301 (2012), 1205.0649.
  * [14] F. Šimkovic, V. Rodin, A. Faessler, and P. Vogel,  Phys. Rev. C 87, 045501 (2013), 1302.1509.
  * [15] N. López Vaquero, T. R. Rodríguez, and J. L. Egido,  Phys. Rev. Lett. 111, 142501 (2013), 1401.0650.
  * [16] J. Barea, J. Kotila, and F. Iachello,  Phys. Rev. C 91, 034304 (2015), 1506.08530.
  * [17] R. Sen’kov and M. Horoi,  Phys. Rev. C 93, 044334 (2016), 1512.06157.
  * [18] J. Engel and J. Menéndez,  Rept. Prog. Phys. 80, 046301 (2017), 1610.06548.
  * [19] M. J. Dolinski, A. W. Poon, and W. Rodejohann,  Ann. Rev. Nucl. Part. Sci. 69, 219 (2019), 1902.04097.
  * [20] S. Weinberg,  Physica A96, 327 (1979).
  * [21] S. Weinberg,  Phys. Lett. B251, 288 (1990).
  * [22] S. Weinberg,  Nucl. Phys. B363, 3 (1991).
  * [23] D. B. Kaplan, M. J. Savage, and M. B. Wise,  Nucl. Phys. B 478, 629 (1996), nucl-th/9605002.
  * [24] D. B. Kaplan, M. J. Savage, and M. B. Wise,  Phys. Lett. B424, 390 (1998), nucl-th/9801034.
  * [25] D. B. Kaplan, M. J. Savage, and M. B. Wise,  Nucl. Phys. B534, 329 (1998), nucl-th/9802075.
  * [26] E. Epelbaum, H.-W. Hammer, and U.-G. Meißner,  Rev. Mod. Phys. 81, 1773 (2009), 0811.1338.
  * [27] R. Machleidt and D. Entem,  Phys. Rept. 503, 1 (2011), 1105.2919.
  * [28] H.-W. Hammer, S. König, and U. van Kolck,  Rev. Mod. Phys. 92, 025004 (2020), 1906.12122.
  * [29] J. Yao et al.,  Phys. Rev. Lett. 124, 232501 (2020), 1908.05424.
  * [30] A. Belley, C. G. Payne, S. R. Stroberg, T. Miyagi, and J. D. Holt,  Phys. Rev. Lett. 126, 042502 (2021), 2008.06588.
  * [31] S. Novario et al.,  (2020), 2008.09696.
  * [32] I. Towner,  Phys. Rept. 155, 263 (1987).
  * [33] S. Pastore et al.,  Phys. Rev. C 97, 022501 (2018), 1709.03592.
  * [34] P. Gysbers et al.,  Nature Phys. 15, 428 (2019), 1903.00047.
  * [35] G. Prézeau, M. Ramsey-Musolf, and P. Vogel,  Phys. Rev. D 68, 034016 (2003), hep-ph/0303205.
  * [36] J. Menéndez, D. Gazit, and A. Schwenk,  Phys. Rev. Lett. 107, 062501 (2011), 1103.3622.
  * [37] V. Cirigliano, W. Dekens, J. de Vries, M. L. Graesser, and E. Mereghetti,  JHEP 12, 082 (2017), 1708.09390.
  * [38] V. Cirigliano, W. Dekens, E. Mereghetti, and A. Walker-Loud,  Phys. Rev. C 97, 065501 (2018), 1710.01729,  [Erratum: Phys. Rev. C 100, 019903 (2019)].
  * [39] S. Pastore et al.,  Phys. Rev. C 97, 014606 (2018), 1710.05026.
  * [40] V. Cirigliano et al.,  Phys. Rev. Lett. 120, 202001 (2018), 1802.10097.
  * [41] L.-J. Wang, J. Engel, and J. M. Yao,  Phys. Rev. C 98, 031301 (2018), 1805.10276.
  * [42] V. Cirigliano, W. Dekens, J. de Vries, M. L. Graesser, and E. Mereghetti,  JHEP 12, 097 (2018), 1806.02780.
  * [43] V. Cirigliano et al.,  Phys. Rev. C 100, 055504 (2019), 1907.11254.
  * [44] W. Dekens, J. de Vries, K. Fuyuto, E. Mereghetti, and G. Zhou,  JHEP 06, 097 (2020), 2002.07182.
  * [45] M. Pavón Valderrama and D. R. Phillips,  Phys. Rev. Lett. 114, 082502 (2015), 1407.0437.
  * [46] V. Cirigliano, W. Detmold, A. Nicholson, and P. Shanahan,  (2020), 2003.08493.
  * [47] G. A. Miller and J. E. Spencer,  Annals Phys. 100, 562 (1976).
  * [48] F. Šimkovic, A. Faessler, H. Müther, V. Rodin, and M. Stauf,  Phys. Rev. C 79, 055501 (2009), 0902.0331.
  * [49] J. Engel, J. Carlson, and R. B. Wiringa,  Phys. Rev. C 83, 034317 (2011), 1101.0554.
  * [50] O. Benhar, R. Biondi, and E. Speranza,  Phys. Rev. C 90, 065504 (2014), 1401.2030.
  * [51] F. Šimkovic, A. Faessler, V. Rodin, P. Vogel, and J. Engel,  Phys. Rev. C 77, 045503 (2008), 0710.2055.
  * [52] X. Feng, L.-C. Jin, X.-Y. Tuo, and S.-C. Xia,  Phys. Rev. Lett. 122, 022001 (2019), 1809.10511.
  * [53] X.-Y. Tuo, X. Feng, and L.-C. Jin,  Phys. Rev. D 100, 094511 (2019), 1909.13525.
  * [54] NPLQCD, W. Detmold and D. Murphy,  (2020), 2004.07404.
  * [55] X. Feng, L.-C. Jin, Z.-Y. Wang, and Z. Zhang,  Phys. Rev. D 103, 034508 (2021), 2005.01956.
  * [56] Z. Davoudi et al.,  Phys. Rept. 900, 1 (2021), 2008.11160.
  * [57] Z. Davoudi and S. V. Kadam,  Phys. Rev. Lett. 126, 152003 (2021), 2012.02083.
  * [58] T. R. Richardson, M. R. Schindler, S. Pastore, and R. P. Springer,  (2021), 2102.02184.
  * [59] V. Cirigliano, W. Dekens, J. de Vries, M. Hoferichter, and E. Mereghetti,  Phys. Rev. Lett. 126, 172002 (2021), 2012.11602.
  * [60] W. Cottingham,  Annals Phys. 25, 424 (1963).
  * [61] H. Harari,  Phys. Rev. Lett. 17, 1303 (1966).
  * [62] G. Ecker, J. Gasser, A. Pich, and E. de Rafael,  Nucl. Phys. B 321, 311 (1989).
  * [63] W. A. Bardeen, J. Bijnens, and J. Gérard,  Phys. Rev. Lett. 62, 1343 (1989).
  * [64] J. F. Donoghue, B. R. Holstein, and D. Wyler,  Phys. Rev. D 47, 2089 (1993).
  * [65] R. Baur and R. Urech,  Phys. Rev. D 53, 6552 (1996), hep-ph/9508393.
  * [66] J. F. Donoghue and A. F. Pérez,  Phys. Rev. D 55, 7075 (1997), hep-ph/9611331.
  * [67] J. Gasser and H. Leutwyler,  Nucl. Phys. B 94, 269 (1975).
  * [68] J. Gasser and H. Leutwyler,  Phys. Rept. 87, 77 (1982).
  * [69] A. Walker-Loud, C. E. Carlson, and G. A. Miller,  Phys. Rev. Lett. 108, 232301 (2012), 1203.0254.
  * [70] A. Thomas, X. Wang, and R. Young,  Phys. Rev. C 91, 015209 (2015), 1406.4579.
  * [71] F. Erben, P. Shanahan, A. Thomas, and R. Young,  Phys. Rev. C 90, 065205 (2014), 1408.6628.
  * [72] J. Gasser, M. Hoferichter, H. Leutwyler, and A. Rusetsky,  Eur. Phys. J. C 75, 375 (2015), 1506.06747,  [Erratum: Eur. Phys. J. C 80, 353 (2020)].
  * [73] J. Gasser, H. Leutwyler, and A. Rusetsky,  Phys. Lett. B 814, 136087 (2021), 2003.13612.
  * [74] J. Gasser, H. Leutwyler, and A. Rusetsky,  Eur. Phys. J. C 80, 1121 (2020), 2008.05806.
  * [75] S. Borsanyi et al.,  Science 347, 1452 (2015), 1406.4088.
  * [76] D. A. Brantley et al.,  (2016), 1612.07733.
  * [77] CSSM, QCDSF, UKQCD, R. Horsley et al.,  J. Phys. G 46, 115004 (2019), 1904.02304.
  * [78] S. Weinberg,  Phys. Rev. Lett. 43, 1566 (1979).
  * [79] R. Furnstahl, H. Hammer, and N. Tirfessa,  Nucl. Phys. A 689, 846 (2001), nucl-th/0010078.
  * [80] M. L. Goldberger and K. M. Watson,  Collision Theory (Wiley, New York, 1964).
  * [81] R. J. Hill, P. Kammel, W. J. Marciano, and A. Sirlin,  Rept. Prog. Phys. 81, 096301 (2018), 1708.08462.
  * [82] J. Arrington and I. Sick,  Phys. Rev. C 76, 035201 (2007), nucl-th/0612079.
  * [83] I. Lorenz, U.-G. Meißner, H.-W. Hammer, and Y. B. Dong,  Phys. Rev. D 91, 014023 (2015), 1411.1704.
  * [84] M. Hoferichter, B. Kubis, J. Ruiz de Elvira, H.-W. Hammer, and U.-G. Meißner,  Eur. Phys. J. A 52, 331 (2016), 1609.06722.
  * [85] Z. Ye, J. Arrington, R. J. Hill, and G. Lee,  Phys. Lett. B 777, 8 (2018), 1707.09063.
  * [86] M. Srivastava and D. Sprung,  Off-Shell Behavior of the Nucleon-Nucleon Interaction (Springer Sicence + Business Media, New York, 1975).
  * [87] D. B. Kaplan and J. V. Steele,  Phys. Rev. C 60, 064002 (1999), nucl-th/9905027.
  * [88] J. Reid, Roderick V.,  Annals Phys. 50, 411 (1968).
  * [89] R. B. Wiringa, V. G. J. Stoks, and R. Schiavilla,  Phys. Rev. C51, 38 (1995), nucl-th/9408016.
  * [90] A. Nicholson et al.,  Phys. Rev. Lett. 121, 172501 (2018), 1805.02634.
  * [91] M. J. Savage,  Phys. Rev. C 59, 2293 (1999), nucl-th/9811087.
  * [92] V. Cirigliano, W. Dekens, M. Graesser, and E. Mereghetti,  Phys. Lett. B 769, 460 (2017), 1701.01443.
  * [93] J. Bijnens and G. Ecker,  Ann. Rev. Nucl. Part. Sci. 64, 149 (2014), 1405.6488.
  * [94] G. Colangelo, M. Hoferichter, M. Procura, and P. Stoffer,  JHEP 09, 074 (2015), 1506.01386.
  * [95] G. Colangelo, M. Hoferichter, M. Procura, and P. Stoffer,  Phys. Rev. Lett. 118, 232001 (2017), 1701.06554.
  * [96] G. Colangelo, M. Hoferichter, M. Procura, and P. Stoffer,  JHEP 04, 161 (2017), 1702.07347.
  * [97] H. Bethe,  Phys. Rev. 76, 38 (1949).
  * [98] J. D. Jackson and J. M. Blatt,  Rev. Mod. Phys. 22, 77 (1950).
  * [99] X. Kong and F. Ravndal,  Nucl. Phys. A 665, 137 (2000), hep-ph/9903523.
  * [100] E. Epelbaum and U.-G. Meißner,  Phys. Lett. B 461, 287 (1999), nucl-th/9902042,  [Erratum: Phys. Lett. B 467, 308 (1999)].
  * [101] J. Bergervoet, P. van Campen, W. van der Sanden, and J. J. de Swart,  Phys. Rev. C 38, 15 (1988).
  * [102] P. Reinert, H. Krebs, and E. Epelbaum,  Eur. Phys. J. A 54, 86 (2018), 1711.08821.
  * [103] S. Klarsfeld, J. Martorell, and D. Sprung,  J. Phys. G 10, 165 (1984).
  * [104] R. Machleidt,  Phys. Rev. C 63, 024001 (2001), nucl-th/0006014.
  * [105] Q. Chen et al.,  Phys. Rev. C 77, 054002 (2008).
  * [106] K. M. Watson,  Phys. Rev. 95, 228 (1954).
  * [107] H.-W. Hammer, A. Nogga, and A. Schwenk,  Rev. Mod. Phys. 85, 197 (2013), 1210.4273.
  * [108] K. Hebeler,  Phys. Rept. 890, 1 (2021), 2002.09548.
  * [109] J. M. Yao et al.,  Phys. Rev. C 103, 014315 (2021), 2010.08609.
  * [110] M. Hoferichter, J. Menéndez, and A. Schwenk,  Phys. Rev. D 102, 074018 (2020), 2007.08529.
  * [111] NEMO-3, SuperNEMO, T. Le Noblet,  J. Phys. Conf. Ser. 1342, 012029 (2020).
  * [112] LEGEND, A. J. Zsigmond,  J. Phys. Conf. Ser. 1468, 012111 (2020).
  * [113] CUPID, B. Schmidt et al.,  J. Phys. Conf. Ser. 1468, 012129 (2020), 1911.10426.
  * [114] nEXO, A. Pocar,  J. Phys. Conf. Ser. 1468, 012131 (2020).
  * [115] CANDLES, K. Tetsuno et al.,  J. Phys. Conf. Ser. 1468, 012132 (2020).
  * [116] SNO+, V. Lozza,  J. Phys. Conf. Ser. 1468, 012135 (2020).
  * [117] KamLAND-Zen, Y. Gando,  J. Phys. Conf. Ser. 1468, 012142 (2020).
  * [118] E. Epelbaum, U.-G. Meißner, W. Glöckle, and C. Elster,  Phys. Rev. C 65, 044001 (2002), nucl-th/0106007.
  * [119] G. Colangelo, M. Hoferichter, and P. Stoffer,  JHEP 02, 006 (2019), 1810.00007.
  * [120] G. Colangelo, M. Hoferichter, and P. Stoffer,  Phys. Lett. B 814, 136073 (2021), 2010.07943.