LA-UR-21-20994

Determining the leading-order contact term in
neutrinoless double $\boldsymbol{\beta}$ decay

Vincenzo Cirigliano,a Wouter Dekens,b Jordy de Vries,c,d,e,f
Martin Hoferichter,g Emanuele Mereghettia

a Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545,
USA

b Department of Physics, University of California at San Diego, La Jolla, CA
92093, USA

c Institute for Theoretical Physics Amsterdam and Delta Institute for
Theoretical Physics, University of Amsterdam, Science Park 904, 1098 XH
Amsterdam, The Netherlands

d Nikhef, Theory Group, Science Park 105, 1098 XG, Amsterdam, The Netherlands

e Amherst Center for Fundamental Interactions, Department of Physics,
University of Massachusetts, Amherst, MA 01003, USA

f RIKEN BNL Research Center, Brookhaven National Laboratory, Upton, New York
11973-5000, USA

g Albert Einstein Center for Fundamental Physics, Institute for Theoretical
Physics, University of Bern, Sidlerstrasse 5, 3012 Bern, Switzerland

We present a method to determine the leading-order (LO) contact term
contributing to the $nn\to ppe^{-}e^{-}$ amplitude through the exchange of
light Majorana neutrinos. Our approach is based on the representation of the
amplitude as the momentum integral of a known kernel (proportional to the
neutrino propagator) 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 analogy to the Cottingham
formula for the electromagnetic contribution to hadron masses. We construct
model-independent representations of the integrand in the low- and high-
momentum regions, through chiral EFT and the operator product expansion,
respectively. We then construct a model for the full amplitude by
interpolating between these two regions, using appropriate nucleon factors for
the weak currents and information on nucleon–nucleon ($N\\!N$) scattering in
the ${}^{1}S_{0}$ channel away from threshold. By matching the amplitude
obtained in this way to the LO chiral EFT amplitude we obtain the relevant LO
contact term and discuss various sources of uncertainty. We validate the
approach by computing the analog $I=2$ $N\\!N$ contact term and by
reproducing, within uncertainties, the charge-independence-breaking
contribution to the ${}^{1}S_{0}$ $N\\!N$ scattering lengths. While our
analysis is performed in the $\overline{\text{MS}}$ scheme, we express our
final result in terms of the scheme-independent renormalized amplitude
$\mathcal{A}_{\nu}(|\mathbf{p}|,|\mathbf{p}^{\prime}|)$ at a set of kinematic
points near threshold. We illustrate for two cutoff schemes how, using our
synthetic data for $\mathcal{A}_{\nu}$, one can determine the contact-term
contribution in any regularization scheme, in particular the ones employed in
nuclear-structure calculations for isotopes of experimental interest.

###### Contents

  1. 1 Introduction
  2. 2 Integral representation and matching strategy
    1. 2.1 Generalities
    2. 2.2 Matching strategy
    3. 2.3 Spectral representation
  3. 3 Chiral EFT result to leading order
  4. 4 Full theory parameterization
    1. 4.1 Low- and intermediate-momentum region: ${\cal A}^{<}$
    2. 4.2 High-momentum region: ${\cal A}^{>}$
  5. 5 Matching and extraction of $\boldsymbol{{\cal C}_{1}}$
    1. 5.1 Results
    2. 5.2 Discussion of uncertainties
  6. 6 Vector–vector amplitude and $\boldsymbol{{\cal C}_{1}+{\cal C}_{2}}$
    1. 6.1 Effective Lagrangians
    2. 6.2 Pion two-point function and the low-energy constant $Z$
    3. 6.3 $nn\to pp$ vector-like amplitude in chiral EFT and full theory
      1. 6.3.1 ${\cal A}_{VV}$ in chiral EFT
      2. 6.3.2 ${\cal A}_{VV}$ in the full theory
    4. 6.4 Matching
    5. 6.5 Charge-independence-breaking contribution to $N\\!N$ scattering
  7. 7 Synthetic data for $\boldsymbol{nn\to pp}$ near threshold
  8. 8 Conclusions
  9. A Half-off-shell $\boldsymbol{T}$ matrix
    1. A.1 The half-off-shell form factor in pionless EFT
      1. A.1.1 Half-off-shell $T$ matrix in pionless EFT at NLO
      2. A.1.2 Behavior of $f_{S}$ and ${\cal I}_{C}^{<}(|{\mathbf{k}}|)$ under field redefinitions
      3. A.1.3 Impact of $f_{S}$ on ${\cal I}_{C}^{<}(|{\mathbf{k}}|)$
    2. A.2 Half-off-shell effects in chiral EFT at NLO
    3. A.3 Half-off-shell effects in $N\\!N$ potential models
  10. B Left–right correlator
  11. C Estimating inelastic effects
  12. D Evaluation of $\boldsymbol{Z^{<}}$
  13. E Details on CIB in $\boldsymbol{N\\!N}$ scattering near threshold
    1. E.1 Low-energy $N\\!N$ scattering in the presence of Coulomb interactions
    2. E.2 Low-energy $N\\!N$ amplitudes in chiral EFT with isospin breaking
    3. E.3 Fitting the couplings
    4. E.4 Predicting the CIB combination of scattering lengths

##  1 Introduction

Neutrinoless double $\beta$ decay ($0\nu\beta\beta$) is the process in which
two neutrons in a nucleus convert into two protons by emitting two electrons
and no neutrinos [1]. This process is by far the most sensitive laboratory
probe of lepton number violation (LNV) and its observation would prove that
neutrinos are Majorana fermions [2], constrain neutrino mass parameters, and
provide experimental validation for leptogenesis scenarios [3, 4]. If
$0\nu\beta\beta$ decay is caused by the exchange of light Majorana neutrinos,
as we assume throughout this paper, the amplitude is proportional to the
effective neutrino mass $m_{\beta\beta}=\sum_{i}U_{ei}^{2}m_{i}$, where the
sum runs over light neutrino masses $m_{i}$ and $U_{ei}$ are elements of the
neutrino mixing matrix. $0\nu\beta\beta$ decay is a complicated process
encompassing aspects from particle, nuclear, and atomic physics, with the
interpretation of current experimental limits [5, 6, 7, 8, 9, 10] and of
potential future discoveries limited by substantial uncertainties in the
calculation of hadronic and nuclear matrix elements [11, 12, 13, 14, 15, 16,
17, 18, 19].

It has been realized in recent years that chiral effective field theory (EFT)
[20, 21, 22, 23, 24, 25] can play a central role in addressing these
uncertainties. Nuclear structure, ab-initio calculations based on chiral-EFT
interactions [26, 27, 28] have recently become available for some
phenomenologically relevant nuclei [29, 30, 31]. In addition, the issue of
$g_{A}$ quenching in single $\beta$ decays has been demonstrated to arise from
the combination of two-nucleon weak currents and strong correlations in the
nucleus [32, 33, 34], and the few-nucleon amplitudes used as input in nuclear
structure calculations have been scrutinized in chiral EFT for various sources
of LNV [35, 36, 37, 38, 39, 40, 41, 42, 43, 44]. In the context of light-
Majorana-neutrino exchange, using naive dimensional counting, the leading
contribution in the chiral-EFT expansion arises from a neutrino-exchange
diagram, in which the LNV arises from insertion of the $\Delta L=2$ effective
neutrino mass $m_{\beta\beta}$. When considering the ${}^{1}S_{0}$ channel, in
analogy to the nucleon–nucleon ($N\\!N$) potential itself [23, 24, 25] and
external currents [45], this conclusion no longer holds when demanding
manifest renormalizability of the amplitude. In fact, it has been shown that
renormalization requires the promotion of an $nn\to ppe^{-}e^{-}$ contact
operator to LO [40, 43], which encodes the exchange of neutrinos with
energy/momentum greater than the nuclear scale and thus cannot be resolved in
chiral EFT. As discussed in greater detail in Refs. [43, 46], the new coupling
encodes a non-factorizable two-nucleon effect, beyond the factorizable one-
nucleon corrections captured by the radii of weak form factors, which also
give a short-range neutrino potential. Moreover, the new short-range coupling
is not captured by the so-called short-range correlations [47, 48, 49, 50], as
it is needed even when one works with fully correlated wave functions, i.e.,
exact solutions of the Schrödinger equation with the appropriate strong
potential. The situation is analogous to single $\beta$ decay, where two-
nucleon weak currents and short-range correlations are both present, and the
combination of both leads to the apparent quenching of $g_{A}$ [33, 34].

To leading order in chiral EFT, a contact interaction is needed only in the
${}^{1}S_{0}$ channel and not in higher partial waves [43]. However, it is
worth emphasizing that the effect of the contact term in the ${}^{1}S_{0}$
channel is amplified in nuclear matrix elements by the cancellation between
the contribution of $N\\!N$ pairs in the ${}^{1}S_{0}$ channel and in states
with higher total angular momentum. This is seen quite dramatically in the
matrix element densities for light nuclear transitions studied in Refs. [40,
43]. These ab-initio results in light nuclei are in qualitative agreement with
the behavior observed in heavy nuclei, first discussed in Ref. [51]. A
complete discussion of the $nn\to pp$ transition operator in chiral EFT can be
found in Ref. [43] (leading order) and Refs. [38, 41] (higher orders).

The value of the short-range coupling in the ${}^{1}S_{0}$ channel then has to
be either extracted from other processes related by chiral symmetry or
calculated from first principles in lattice QCD [52, 53, 46, 54, 55, 56, 57]
(see Ref. [58] for a large-$N_{c}$ analysis). Currently, however, the size of
this contact operator is unknown, leading to substantial uncertainties in the
interpretation of $0\nu\beta\beta$ decays besides the nuclear-structure ones,
especially given that its impact is enhanced in $\Delta I=2$ nuclear
transitions due to a node in the matrix element density [40, 43]. In this work
we present in some detail the method used to obtain a first estimate of the
complete $nn\to ppe^{-}e^{-}$ amplitude including this contact-term
contribution [59].

The hadronic component of the light-Majorana-neutrino-exchange amplitude has
the structure

${\cal
A}_{\nu}\propto\int\frac{d^{4}k}{(2\pi)^{4}}\frac{g_{\alpha\beta}}{k^{2}+i\epsilon}\int
d^{4}x\,e^{ik\cdot x}\langle
pp|T\\{j_{\text{w}}^{\alpha}(x)j_{\text{w}}^{\beta}(0)\\}|nn\rangle\,,$ (1.1)

and is controlled by the two-nucleon matrix element of the time-ordered
product $T\\{j_{\text{w}}^{\alpha}(x)j_{\text{w}}^{\beta}(0)\\}$ of two weak
currents. Such matrix elements with the weak current replaced by the
electromagnetic current $j_{\text{em}}^{\alpha}(x)$ appear in the
electromagnetic contributions to hadron masses or scattering processes, in
which case a relation exists between the forward Compton scattering amplitude
and its contraction with a massless propagator, as given in Eq. (1.1). The
relation, known as the Cottingham formula [60, 61], has been used to estimate
the electromagnetic contributions to the masses of pions [62, 63, 64, 65, 66]
and nucleons [67, 68, 69, 70, 71, 72, 73, 74]. Since the matrix elements in
these cases have precisely the same structure as required for the light-
Majorana-neutrino-exchange contribution to the $0\nu\beta\beta$ decay $nn\to
ppe^{-}e^{-}$, our method aims to constrain the corresponding amplitude by
generalizing the Cottingham approach to the two-nucleon system, and then
determine the contact-term contribution by matching to chiral EFT.

The application of the Cottingham approach to the pion and nucleon mass
difference has a long history and in both cases the by far dominant
contribution arises from elastic intermediate states. The pion-pole
contribution gives more than $80\%$ of the pion mass difference [62, 63, 64,
65, 66] and, similarly, the nucleon pole provides the bulk of the
electromagnetic part of the proton–neutron mass difference
$m_{p-n}^{\text{el}}=0.75(2)\,\text{MeV}$. Despite the tension between the
estimate of the inelastic contributions in lattice QCD,
$m_{p-n}^{\text{inel}}=0.28(11)\,\text{MeV}$ [75, 76, 77], and from nucleon
structure functions, $m_{p-n}^{\text{inel}}=-0.17(16)\,\text{MeV}$ [73, 74,
72, 67], also in this case the elastic estimate is accurate at the $30\%$
level. The main complication in the generalization to $0\nu\beta\beta$ decay
arises from the two-particle nature of initial and final states. First, due to
the proliferation of kinematic variables and scalar functions in a Lorentz
decomposition of the general amplitude, it becomes extremely cumbersome to try
and set up a strict derivation of the elastic contribution via dispersion
relations. Second, the $N\\!N$ scattering amplitude itself gives rise to an
additional source of momentum dependence that adds to the physics included in
terms of pion and nucleon form factors in the standard Cottingham approach.
Accordingly, we do not attempt a comprehensive analysis of all scalar
functions describing the full two-particle problem, but instead include the
most important intermediate states in terms of the respective form factors—in
close analogy to the elastic results for the pion and nucleon Cottingham
formula—as well as the momentum dependence of the $N\\!N$ scattering
amplitude. To validate this approach, we also consider the two-nucleon matrix
element with two electromagnetic currents, which can be accessed
experimentally in terms of charge independence breaking (CIB) in the $N\\!N$
scattering lengths. Comparison with data then allows us to confirm the
expectation of an accuracy around $30\%$ if only elastic contributions are
kept, as suggested by the Cottingham result for the proton–neutron mass
difference. A determination at this level already has a valuable impact in
bounding the size of the contact-term contribution to $0\nu\beta\beta$ decay.

The derivation is organized as follows: in Section 2 we present the general
integral representation of the amplitude and our matching strategy. In Section
3 we recast the LO chiral EFT amplitude in a form suitable for matching
purposes. In Section 4 we present our construction of our full $nn\to pp$
amplitude, followed by matching to the EFT result and extraction of the
contact term in Section 5. In Section 6 we present the analysis for the two-
nucleon $I=2$ electromagnetic amplitude and the validation of the method
through comparison with the experimental data on CIB in the $N\\!N$ scattering
lengths. This section is fairly technical and can be omitted by readers
primarily interested in $0\nu\beta\beta$ decay. In Section 7 we return to the
LNV amplitude $nn\to pp$ and present synthetic data at kinematic points near
threshold, illustrating how these can be used to extract the contact term in
any regularization and renormalization scheme. We present our concluding
remarks in Section 8. Details on the half-off-shell behavior of the $N\\!N$
scattering amplitude (Appendix A), the operator product expansion (OPE)
(Appendix B), the size of typical inelastic contributions (Appendix C), the
electromagnetic pion mass splitting (Appendix D), and the CIB $N\\!N$
scattering lengths (Appendix E) are provided in the appendices.

##  2 Integral representation and matching strategy

###  2.1 Generalities

Including the effect of LNV from the dimension-five Weinberg operator [78],
the low-energy effective Lagrangian at scale $\mu\gtrsim\Lambda_{\chi}\sim
1\,\text{GeV}$ is given by

${\cal L}_{\rm eff}={\cal L}_{\rm QCD}\ -\left\\{\ 2\sqrt{2}G_{F}V_{ud}\
\bar{u}_{L}\gamma^{\mu}d_{L}\,\bar{e}_{L}\gamma_{\mu}\nu_{eL}\ +\
\frac{1}{2}m_{\beta\beta}^{*}\ {\nu}_{eL}^{T}C\nu_{eL}\ -\ C_{L}\,O_{L}+{\rm
h.c.}\right\\}\,.$ (2.1)

The second term in Eq. (2.1) represents the Fermi charged-current weak
interaction. The last two terms encode LNV through the neutrino Majorana mass,
given by $m_{\beta\beta}=\sum_{i}U_{ei}^{2}m_{i}$ in terms of mass eigenstates
and elements of the neutrino mixing matrix,111The effective mass probed in
$0\nu\beta\beta$ decay is often defined as
$m_{\beta\beta}=\big{|}\sum_{i}U_{ei}^{2}m_{i}\big{|}$, but to simplify the
notation at the Lagrangian level we formally keep its phase. and a dimension-
nine $\Delta L=2$ operator generated at the electroweak threshold:

$O_{L}=\bar{e}_{L}e_{L}^{c}\ \bar{u}_{L}\gamma_{\mu}d_{L}\
\bar{u}_{L}\gamma^{\mu}d_{L}\equiv\bar{e}_{L}e_{L}^{c}\ O_{1}\,,$ (2.2)

with $e_{L}^{c}=C\bar{e}_{L}^{T}$. Since
$C_{L}=(8V_{ud}^{2}G_{F}^{2}m_{\beta\beta})/M_{W}^{2}\times(1+\mathcal{O}(\alpha_{s}/\pi))$,
the effect of the latter term on the $0\nu\beta\beta$ amplitude is suppressed
by $(k_{F}/M_{W})^{2}$ (where $k_{F}\sim\mathcal{O}(100)\,\text{MeV}$ is the
typical Fermi momentum of nucleons in a nucleus) compared to light-neutrino
exchange and can be safely neglected at this stage. However, the isotensor
four-quark local operator $O_{1}$ itself will play an important role in the
following analysis.

The interactions of Eq. (2.1) induce $\Delta L=2$ transitions (such as
$\pi^{-}\pi^{-}\to e^{-}e^{-}$, $nn\to ppe^{-}e^{-}$, ${}^{76}{\rm
Ge}\to\,^{76}{\rm Se}\,e^{-}e^{-}$, ${}^{136}{\rm Xe}\to\,^{136}{\rm
Ba}\,e^{-}e^{-}$, …) through the non-local effective action obtained by
contracting the neutrino fields in the two weak vertices,

$S_{\rm eff}^{\Delta L=2}=\frac{8G_{F}^{2}V_{ud}^{2}m_{\beta\beta}}{2!}\int
d^{4}xd^{4}y\
S(x-y)\times\bar{e}_{L}(x)\gamma^{\mu}\gamma^{\nu}e_{L}^{c}(y)\times
T\Big{(}\bar{u}_{L}\gamma_{\mu}d_{L}(x)\
\bar{u}_{L}\gamma_{\nu}d_{L}(y)\Big{)}\,\,,$ (2.3)

where

$S(r)=\int\frac{d^{4}k}{(2\pi)^{4}}\frac{e^{-ik\cdot r}}{k^{2}+i\epsilon}$
(2.4)

is the scalar massless propagator, a remnant of the neutrino propagator.
Computing matrix elements of $S_{\rm eff}^{\Delta L=2}$ in hadronic and
nuclear states is a notoriously difficult task. The multi-scale nature of the
problem can be seen more explicitly by going to the Fourier representation
222In deriving Eq. (2.5) from Eq. (2.3) we approximate
$\bar{e}_{L}(x)\gamma^{\mu}\gamma^{\nu}e_{L}^{c}(y)\simeq\bar{e}_{L}(x)\gamma^{\mu}\gamma^{\nu}e_{L}^{c}(x)=g^{\mu\nu}\bar{e}_{L}(x)e_{L}^{c}(x)$,
which amounts to neglecting the difference in electron momenta, a safe
assumption given that $|p_{e1}-p_{e2}|/k_{F}\ll 1$.

$\langle e_{1}e_{2}h_{f}|S_{\rm eff}^{\Delta
L=2}|h_{i}\rangle=\frac{8G_{F}^{2}V_{ud}^{2}m_{\beta\beta}}{2!}\int d^{4}x\
\langle
e_{1}e_{2}|\bar{e}_{L}(x)e_{L}^{c}(x)|0\rangle\int\frac{d^{4}k}{(2\pi)^{4}}\frac{g^{\mu\nu}\overline{T}_{\mu\nu}(k,p_{\rm
ext},x)}{k^{2}+i\epsilon}\,\,,\ \ $ (2.5)

where

$\displaystyle\overline{T}_{\mu\nu}(k,p_{\rm ext},x)$ $\displaystyle=\langle
h_{f}(p_{f})|\ \hat{\Pi}_{\mu\nu}^{LL}(k,x)\ |h_{i}(p_{i})\rangle\,,$
$\displaystyle\hat{\Pi}_{\mu\nu}^{LL}(k,x)$ $\displaystyle=\int
d^{4}r\,e^{ik\cdot r}\ T\Big{\\{}J_{\mu}^{L}(x+r/2)\
J_{\nu}^{L}(x-r/2)\Big{\\}}\,,\qquad
J^{L}_{\mu}=\bar{u}_{L}\gamma_{\mu}d_{L}\,,$ (2.6)

and $p_{\rm ext}$ denotes generically the hadronic external momenta $p_{f}$
and $p_{i}$. Using translational invariance one has

$\displaystyle\overline{T}_{\mu\nu}(k,p_{\rm ext},x)$ $\displaystyle=\langle
h_{f}(p_{f})|\ e^{ix\cdot P}\ \hat{\Pi}_{\mu\nu}^{LL}(k,0)\ e^{-ix\cdot P}\
|h_{i}(p_{i})\rangle$ $\displaystyle=e^{ix\cdot(p_{f}-p_{i})}\langle
h_{f}(p_{f})|\ \hat{\Pi}_{\mu\nu}^{LL}(k,0)\ |h_{i}(p_{i})\rangle$
$\displaystyle=e^{ix\cdot(p_{f}-p_{i})}\ \overline{T}_{\mu\nu}(k,p_{\rm
ext},0)\,.$ (2.7)

Therefore, defining333In what follows we will suppress the space-time label in
the correlator: $\hat{\Pi}^{LL}_{\mu\nu}(k,0)\to\hat{\Pi}^{LL}_{\mu\nu}(k)$.

$T_{\mu\nu}(k,p_{\rm ext})\equiv\overline{T}_{\mu\nu}(k,p_{\rm ext},0)=\langle
h_{f}(p_{f})|\ \hat{\Pi}_{\mu\nu}^{LL}(k,0)\ |h_{i}(p_{i})\rangle\,,$ (2.8)

one arrives at

$\displaystyle\langle e_{1}e_{2}h_{f}|S_{\rm eff}^{\Delta L=2}|h_{i}\rangle$
$\displaystyle=(2\pi)^{4}\,\delta^{(4)}(p_{e1}+p_{e2}+p_{f}-p_{i})\Big{(}4G_{F}^{2}V_{ud}^{2}m_{\beta\beta}\
\bar{u}_{L}(p_{1})u_{L}^{c}(p_{2})\Big{)}\times{\cal A}_{\nu}\,,$
$\displaystyle{\cal A}_{\nu}$
$\displaystyle=2\,\int\frac{d^{4}k}{(2\pi)^{4}}\frac{{\cal T}(k,p_{\rm
ext})}{k^{2}+i\epsilon}\,,\qquad{\cal T}(k,p_{\rm ext})\equiv
g^{\mu\nu}T_{\mu\nu}(k,p_{\rm ext})\,.$ (2.9)

The hadronic amplitude ${\cal A}_{\nu}$ in Eq. (2.9) receives contributions
from neutrino virtualities $k^{2}$ ranging from the weak scale all the way
down to the infrared (IR) scale of nuclear bound states.

To estimate the LO contact term arising in chiral EFT, we will employ the
representation (2.9) to obtain the amplitude in the “full theory,” and then
match to the appropriate EFT expression. Since the contact term arises in the
${}^{1}S_{0}$ channel, we will take as external states $nn$ and $pp$ in the
${}^{1}S_{0}$ state and $T_{\mu\nu}(k,p_{\rm ext})$ will be thought of as the
generalized forward Compton amplitude

$n(p_{1})\ n(p_{2})\ W^{+}(k)\ \to\ p(p_{1}^{\prime})\ p(p_{2}^{\prime})\
W^{-}(k)\,.$ (2.10)

Since the low-energy constants (LECs) do not depend on the IR details, we will
perform the matching calculation at the simplest kinematic point, in which the
two electrons are emitted with zero three-momentum in the center-of-mass frame
of the incoming neutron pair [40, 43]. Explicitly we have

$\displaystyle p_{1\mu}$ $\displaystyle=(E,{\mathbf{p}})\,,$ $\displaystyle
p_{2\mu}$ $\displaystyle=(E,-{\mathbf{p}})\,,$ $\displaystyle E$
$\displaystyle=\sqrt{{\mathbf{p}}^{2}+m_{n}^{2}}\,,$ $\displaystyle
p^{\prime}_{1\mu}$ $\displaystyle=(E^{\prime},{\mathbf{p}}^{\prime})\,,$
$\displaystyle p^{\prime}_{2\mu}$
$\displaystyle=(E^{\prime},-{\mathbf{p}}^{\prime})\,,$ $\displaystyle
E^{\prime}$ $\displaystyle=\sqrt{{\mathbf{p}}^{\prime 2}+m_{p}^{2}}\,,$ (2.11)

where $2E=2E^{\prime}+2m_{e}$. Free two-nucleon states with vanishing total
three-momentum and individual three-momenta given by $\pm{\mathbf{q}}$ will be
denoted by $|{\mathbf{q}}\rangle$ [23], so for example for the initial and
final state we will have $|i_{0}\rangle=|{\mathbf{p}}\rangle$ and
$|f_{0}\rangle=|{\mathbf{p}}^{\prime}\rangle$, respectively.

###  2.2 Matching strategy

The amplitude for the process $nn\to pp$ is given in Eq. (2.9) as the integral
of the product of a massless propagator (remnant of the Majorana neutrino
propagator) with the contracted hadronic tensor ${\cal T}(k,p_{\rm
ext})=g_{\mu\nu}T^{\mu\nu}$. The neutrino four-momentum regions relevant for
the integration over $d^{4}k$ are schematically depicted in Fig. 1. Denoting
the Euclidean four-momentum squared by
$k_{E}^{2}\equiv(k^{0})^{2}+{\mathbf{k}}^{2}$, one can introduce hard
($k_{E}^{2}>\Lambda^{2}$), intermediate
($\Lambda_{\chi}^{2}<k_{E}^{2}<\Lambda^{2}$), and low-energy
($k_{E}^{2}<\Lambda_{\chi}^{2}$) regions, separated by $\Lambda_{\chi}$ (the
breakdown scale of the low-energy hadronic EFT) and $\Lambda$ (scale at which
the OPE becomes reliable). The low-energy region further includes the soft
($|k^{0}|\sim|{\mathbf{k}}|\sim k_{F}$), potential ($|k^{0}|\sim
k_{F}^{2}/m_{N}$, $|{\mathbf{k}}|\sim k_{F}$), and ultrasoft
($|k^{0}|\sim|{\mathbf{k}}|\ll k_{F}$) regions, essential to reproduce the IR
behavior of the amplitude.

The basic idea behind our approach is that model-independent representations
of the integrand in Eq. (2.9) can be constructed in the low-energy region (via
pionless and chiral EFT) and in the hard region (via the OPE). Given this, a
model for the full amplitude can be constructed by interpolating between these
two regions. This approach uses model-dependent input for the intermediate
momentum region, which we anchor to known constraints from QCD at low and high
momenta.

Figure 1: Schematic representation of the regions of neutrino virtuality
contributing to the amplitude in Eq. (2.9). The boundaries between various
regions are given by $k_{F}\sim 100\,\text{MeV}$, $\Lambda_{\chi}\lesssim
1\,\text{GeV}$, and $\Lambda\gtrsim 1.5\,\text{GeV}$.

In practice, given the non-relativistic nature of the process of interest, we
will not use $k_{E}^{2}$ as matching variable. Instead, we decompose
$d^{4}k=dk^{0}d^{3}{\mathbf{k}}$, first perform the $k^{0}$ integral in the
appropriate regions via Cauchy’s theorem, and then carry out the angular
integrations in $d^{3}{\mathbf{k}}$ to reduce the amplitude to an integral
over $d|{\mathbf{k}}|$. To LO in the expansion in external momenta
$|\mathbf{p}|,|\mathbf{p}^{\prime}|\sim Q$ (we denote by $\mu_{\chi}\sim Q$
the soft scale in the EFT), we write the full amplitude as an integral over
the internal neutrino three-momentum ${\mathbf{k}}$, which we split into a
low- plus intermediate-momentum region and a high-momentum region

$\displaystyle{\cal A}_{\nu}^{\rm full}$
$\displaystyle=\int_{0}^{\infty}\,d|{\mathbf{k}}|\ a^{\rm
full}(|{\mathbf{k}}|)={\cal A}^{<}+{\cal A}^{>}\,,$ $\displaystyle{\cal
A}^{<}$ $\displaystyle=\int_{0}^{\Lambda}\,d|{\mathbf{k}}|\
a_{<}(|{\mathbf{k}}|)\,,$ $\displaystyle{\cal A}^{>}$
$\displaystyle=\int_{\Lambda}^{\infty}\,d|{\mathbf{k}}|\
a_{>}(|{\mathbf{k}}|)\,,$ (2.12)

separated by the scale $\Lambda$ that represents the onset of the asymptotic
behavior for the current–current correlator, controlled by the OPE. This
representation introduces model dependence through: (i) The choices made to
extend the model-independent integrand $a_{\chi}(|{\mathbf{k}}|)$ dictated by
chiral EFT in the region $|{\mathbf{k}}|<\Lambda_{\chi}$ to the function
$a_{<}(|{\mathbf{k}}|)$ valid up to $|{\mathbf{k}}|\sim\Lambda$. We will
provide a simple parameterization of $a_{<}(|{\mathbf{k}}|)$ that reduces to
$a_{\chi}(|{\mathbf{k}}|)$ for $|{\mathbf{k}}|<M_{\pi}$ and incorporates
phenomenological input such as nucleon form factors of the weak current and
resonance contributions to the strong-interaction potential. (ii) The choice
of $\Lambda$ that determines the boundary of integration regions in the
variable $|{\mathbf{k}}|$. Once a representation for ${\cal A}_{\nu}^{\rm
full}$ is obtained, along with an estimate of the associated uncertainties, we
will estimate the LEC appearing in ${\cal A}_{\nu}^{\chi{\rm EFT}}$ by
enforcing the matching condition

${\cal A}_{\nu}^{\chi{\rm EFT}}={\cal A}^{<}+{\cal A}^{>}\,.$ (2.13)

In the following sections, we will describe the construction of ${\cal
A}^{<,>}$ and the matching to ${\cal A}_{\nu}^{\rm\chi{\rm EFT}}$, starting
with the spectral representation in Section 2.3.

###  2.3 Spectral representation

In this section we provide a spectral representation for the $nn\to pp$
amplitude, which will prove very useful in identifying and organizing the
various intermediate-state contributions to ${\cal A}^{<}$ and carrying out
the analysis in analogy to the Cottingham formula [60, 61].

We begin by recalling some elements of the formal theory of scattering that we
will use in various parts of the discussion. We denote by
$\hat{H}=\hat{H}_{0}+\hat{V}$ the total Hamiltonian, split into a free and
interaction term ($\hat{V}$, not to be confused with the potential). Retarded
and advanced Green’s functions in the interacting theory are given by

$\hat{G}_{\pm}(E)=\frac{1}{E-\hat{H}\pm i\epsilon}=\int
dt\,e^{iEt}\hat{G}_{\pm}(t)~{},\qquad\hat{G}_{\pm}(t)=\mp i\theta(\pm
t)\,e^{-i\hat{H}t}\,,$ (2.14)

with analogous definitions for the free-theory ones, denoted by
$\hat{G}_{\pm}^{(0)}(E)$, with the replacement $\hat{H}\to\hat{H}_{0}$. The
scattering operator $\hat{T}(E)$ is formally given by

$\hat{T}(E)=\hat{V}\ \left(I-\hat{G}^{(0)}_{+}(E)\,\hat{V}\right)^{-1}$ (2.15)

and satisfies $\hat{G}_{+}(E)\hat{V}=\hat{G}^{(0)}_{+}(E)\hat{T}(E)$. The
scattering states $\langle f_{-}|$ and $|i_{+}\rangle$ are related to the free
states $\langle f_{0}|$ and $|i_{0}\rangle$ by

$\displaystyle|i_{+}\rangle$
$\displaystyle=\left(I+\hat{G}_{+}^{(0)}(E)\,\hat{T}(E)\right)|i_{0}\rangle\,,$
$\displaystyle\langle f_{-}|$ $\displaystyle=\langle
f_{0}|\left(\hat{T}(E^{\prime})\hat{G}_{+}^{(0)}(E^{\prime})+I\right)\,.$
(2.16)

In terms of the scattering states, the amplitude for $nn\to pp$ can be written
as

${\cal A}_{\nu}=\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\ \langle
f_{-}|\,\hat{O}^{LL}({\mathbf{k}})\,|i_{+}\rangle\,,$ (2.17)

with the weak transition operator

$\hat{O}^{LL}(\mathbf{k})\equiv
2\int\frac{dk^{0}}{2\pi}\,\frac{g^{\mu\nu}\hat{\Pi}_{\mu\nu}^{LL}(k)}{k^{2}+i\epsilon}\,.$
(2.18)

From the definition of the correlator in Eq. (2.6) one obtains the following
representation for $\hat{\Pi}_{\mu\nu}^{LL}(k)$ in terms of Green’s functions:

$g^{\mu\nu}\hat{\Pi}_{\mu\nu}^{LL}(k)=i(2\pi)^{3}\,J_{\mu}^{L}(0)\,\left[\hat{G}_{+}(k_{+}^{0})\,\delta^{(3)}(\hat{{\mathbf{P}}}-{\mathbf{k}}_{+})+\hat{G}_{+}(k_{-}^{0})\,\delta^{(3)}(\hat{{\mathbf{P}}}-{\mathbf{k}}_{-})\right]\,J^{L\mu}(0)\,,$
(2.19)

where $\hat{{\mathbf{P}}}$ is the total three-momentum operator and we have
introduced the four-vectors

$k_{\pm}^{\mu}=\tilde{p}^{\mu}\pm
k^{\mu}\,,\qquad\tilde{p}^{\mu}=\frac{1}{2}(p_{i}+p_{f})^{\mu}=(\tilde{E},\tilde{{\mathbf{p}}})\,.$
(2.20)

The labels in $p_{i,f}^{\mu}$ refer to the initial and final states between
which $\hat{\Pi}^{LL}_{\mu\nu}(k)$ is evaluated. Since we are considering two-
nucleon external states with vanishing total three-momentum (and total
momentum is conserved at each vertex) we have $\tilde{{\mathbf{p}}}=0$ and
hence ${\mathbf{k}}_{\pm}=\pm{\mathbf{k}}$. In Eq. (2.19) the dependence on
$k^{0}$ is very simple, as $k^{0}$ appears only through the energy
denominators of $\hat{G}_{+}(k_{\pm}^{0})$. Performing the integration over
$k^{0}$ in Eq. (2.18) with Cauchy’s theorem,444 For each term in Eq. (2.19),
one can close the contour in the upper or lower $k^{0}$ plane so that the
integral is given by the residue at the $k^{0}$ pole from the neutrino
propagator in Eq. (2.18).  one arrives at

$\hat{O}^{LL}({\mathbf{k}})=\frac{1}{|{\mathbf{k}}|}\
J_{\mu}^{L}(0)\,\hat{G}_{+}(\tilde{E}-|{\mathbf{k}}|)\
(2\pi)^{3}\left[\delta^{(3)}(\hat{{\mathbf{P}}}-{\mathbf{k}}_{+})+\delta^{(3)}(\hat{{\mathbf{P}}}-{\mathbf{k}}_{-})\right]\,J^{L\mu}(0)\,.$
(2.21)

Further inserting a complete set of states between the current operators in
Eq. (2.21) leads to the spectral representation for the amplitude555The
summation is over intermediate states $|n(\mathbf{k}_{\pm})\rangle$ of total
three-momentum $\mathbf{k}_{\pm}$, enforced by the $\delta$-functions in Eq.
(2.21). Therefore, for an $N$-particle intermediate state
$\scriptsize{\sum_{n}}$ involves phase space integrals over the $N-1$ internal
momenta (the total momentum being fixed to $\mathbf{k}_{\pm}$) and carries
non-zero mass dimension. For example, for two-nucleon intermediate states,
using non-relativistic normalizations for the states
$\langle\mathbf{p}_{n}|\mathbf{p}_{n}^{\prime}\rangle=(2\pi)^{3}\delta^{(3)}(\mathbf{p}_{n}-\mathbf{p}_{n}^{\prime})$
one has $\sum_{n}\to\int d^{3}\mathbf{p}_{n}/(2\pi)^{3}$, where
$\mathbf{p}_{n}$ is the relative momentum of the two-nucleon pair. In general
the summation $\sum_{n}\ |n(\mathbf{k}_{\pm})\rangle\langle
n(\mathbf{k}_{\pm})|$ carries mass dimension $-3$.

${\cal
A}_{\nu}=-\sum_{n}\int\frac{d^{3}\bf{k}}{(2\pi)^{3}}\,\frac{1}{|{\mathbf{k}}|}\,\left[\frac{\langle
f_{-}|J^{L}_{\mu}|n({\mathbf{k}}_{+})\rangle\langle
n({\mathbf{k}}_{+})|J^{L\mu}|i_{+}\rangle}{|{\mathbf{k}}|+(E_{n}({\mathbf{k}}_{+})-\tilde{E})-i\epsilon}\
+\ \frac{\langle f_{-}|J^{L}_{\mu}|n({\mathbf{k}}_{-})\rangle\langle
n({\mathbf{k}}_{-})|J^{L\mu}|i_{+}\rangle}{|{\mathbf{k}}|+(E_{n}({\mathbf{k}}_{-})-\tilde{E})-i\epsilon}\right]\,.$
(2.22)

The representations (2.17) and (2.22) are quite general. The asymptotic
behavior of the integrand in Eq. (2.22) at large $|{\mathbf{k}}|$ is dictated
by the OPE for $\hat{\Pi}^{LL}_{\mu\nu}(k)$ or, equivalently,
$\hat{O}^{LL}({\mathbf{k}})$. An explicit calculation to be described below
shows the behavior $d^{3}{\mathbf{k}}/|{\mathbf{k}}|^{5}$, so the amplitude in
the full theory is finite. Moreover, Eq. (2.22) shows that once
$|{\mathbf{k}}|>k_{F}$, so that ${\mathbf{k}}^{2}/m_{N}$ is above the typical
nuclear binding energies, one expects
$(E_{n}({\mathbf{k}}_{\pm})-\tilde{E})>0$ even for bound intermediate states
(such as the deuteron), and therefore the energy denominators in Eq. (2.22)
will not lead to any singular behavior in the variable $|{\mathbf{k}}|$. The
matrix elements in the numerator are also expected to have a smooth behavior
in $|{\mathbf{k}}|$, dictated by single- and multi-hadron form factors, as
shown by explicit EFT calculations. Based on these considerations, we conclude
that a smooth interpolation between the calculable regimes of
$|{\mathbf{k}}|\lesssim\Lambda_{\chi}$ and $|{\mathbf{k}}|\gtrsim\Lambda$ is
adequate.

In order to make the integrand in Eqs. (2.17) and (2.22) more explicit, we use
the expression for the scattering states (2.16) in Eq. (2.17) and arrive at

$\displaystyle{\cal A}_{\nu}$
$\displaystyle=\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\ \langle
f_{0}|\left(\hat{T}(E^{\prime})\hat{G}_{+}^{(0)}(E^{\prime})+I\right)\,\hat{O}^{LL}({\mathbf{k}})\,\left(I+\hat{G}_{+}^{(0)}(E)\,\hat{T}(E)\right)|i_{0}\rangle$
(2.23) $\displaystyle=\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\
\Bigg{\\{}\langle f_{0}|\,\hat{O}^{LL}({\mathbf{k}})\,|i_{0}\rangle$
$\displaystyle\qquad+\sum_{m}\ \langle
f_{0}|\hat{T}(E^{\prime})|m\rangle\left[G^{(0)}_{+}(E^{\prime})\right]_{m}\langle
m|\,\hat{O}^{LL}({\mathbf{k}})\,|i_{0}\rangle$ $\displaystyle\qquad+\sum_{m}\
\langle
f_{0}|\,\hat{O}^{LL}({\mathbf{k}})\,|m\rangle\left[G^{(0)}_{+}(E)\right]_{m}\langle
m|\hat{T}(E)|i_{0}\rangle$ $\displaystyle\qquad+\sum_{m,m^{\prime}}\ \langle
f_{0}|\hat{T}(E^{\prime})|m^{\prime}\rangle\left[G^{(0)}_{+}(E^{\prime})\right]_{m^{\prime}}\langle
m^{\prime}|\,\hat{O}^{LL}({\mathbf{k}})\,|m\rangle\left[G^{(0)}_{+}(E)\right]_{m}\langle
m|\ \hat{T}(E)\ |i_{0}\rangle\Bigg{\\}}\,,$

where

$\left[G^{(0)}_{+}(E)\right]_{m}=\frac{1}{E-E^{(0)}_{m}+i\epsilon}\,,$ (2.24)

$E^{(0)}_{m}$ denotes the energy associated with the free Hamiltonian
$\hat{H}_{0}$, and

$\langle
a|\,\hat{O}^{LL}({\mathbf{k}})\,|b\rangle=-\sum_{n}\,\frac{1}{|\mathbf{k}|}\,\left[\frac{\langle
a|J^{L}_{\mu}|n({\mathbf{k}})\rangle\langle
n({\mathbf{k}})|J^{L\mu}|b\rangle}{|{\mathbf{k}}|+(E_{n}({\mathbf{k}})-\tilde{E}_{ab})-i\epsilon}\
+\ \frac{\langle a|J^{L}_{\mu}|n(-{\mathbf{k}})\rangle\langle
n(-{\mathbf{k}})|J^{L\mu}|b\rangle}{|{\mathbf{k}}|+(E_{n}(-{\mathbf{k}})-\tilde{E}_{ab})-i\epsilon}\right]\,,$
(2.25)

with $\tilde{E}_{ab}=(E_{a}+E_{b})/2$. In general the sum over hadronic
intermediate states in Eqs. (2.23) and (2.25) involves $|m\rangle$,
$|m^{\prime}\rangle$, $|n\rangle$, $\in\\{N\\!N,N\\!N\pi,...\\}$, i.e., both
elastic $|N\\!N\rangle$ contributions and inelastic contributions
$|m\rangle,|m^{\prime}\rangle,|n\rangle\neq|N\\!N\rangle$. Equations
(2.23)–(2.25) make it explicit which dynamical input is needed for the
evaluation of the $nn\to pp$ amplitude:

  1. 1.

One needs the matrix elements of the current–current operator
$\hat{O}^{LL}({\mathbf{k}})$ among two-nucleon and possibly other intermediate
states, namely $\langle m^{\prime}|\hat{O}^{LL}({\mathbf{k}})|m\rangle$, that
can be further decomposed according to Eq. (2.25).

  2. 2.

One needs the $T$-matrix elements $\langle m|\hat{T}(E)|i_{0}\rangle$ and
$\langle f_{0}|\hat{T}(E^{\prime})|m^{\prime}\rangle$ involving arbitrary
intermediate states $\langle m|$, $|m^{\prime}\rangle$ and on-shell two-
nucleon states $|i_{0}\rangle=|\mathbf{p}\rangle$ (with
$E=\mathbf{p}^{2}/m_{N}$) and $\langle f_{0}|=\langle\mathbf{p}^{\prime}|$
(with $E^{\prime}=\mathbf{p}^{\prime 2}/m_{N}$). When considering the elastic
contributions, these reduce to the so-called half-off-shell (HOS) $T$-matrix
elements $\langle\mathbf{p}_{m}|\hat{T}(E)|\mathbf{p}\rangle$ and
$\langle\mathbf{p}^{\prime}|\hat{T}(E^{\prime})|\mathbf{p}_{m^{\prime}}\rangle$,
involving loop momenta $\mathbf{p}_{m}$ and $\mathbf{p}_{m^{\prime}}$. While
it is well known that the HOS $T$-matrix elements by themselves are not
physical quantities (see for example the discussion in Ref. [79]), they enter
Eq. (2.23) in such a way that the full physical amplitude ${\cal A}_{\nu}$ is
free of off-shell ambiguities (see Appendix A for an explicit check of this
point).

To LO in chiral EFT the amplitude $\mathcal{A}_{\nu}$ is saturated by elastic
contributions, with all inputs in Eqs. (2.23)–(2.25) given to leading chiral
order. The LO chiral input provides a good representation of the low-momentum
part of the integrand but misrepresents the high-momentum component. In this
language, the ultraviolet (UV) divergence and the need for a LO contact term
arises from the $1/|\mathbf{k}|$ behavior of the integrand, as discussed in
Section 3.

On the other hand, in our estimate of the full amplitude to be described in
Section 4, we will start from Eqs. (2.23)–(2.25) and use representations of
$\langle m^{\prime}|\hat{O}^{LL}({\mathbf{k}})|m\rangle$,
$\langle\mathbf{p}_{m}|\hat{T}(E)|\mathbf{p}\rangle$, and
$\langle\mathbf{p}^{\prime}|\hat{T}(E^{\prime})|\mathbf{p}_{m^{\prime}}\rangle$
that go beyond leading chiral order to construct a UV convergent integrand.
Motivated by the leading chiral EFT analysis and the analogy with the
Cottingham approach to the pion and nucleon electromagnetic mass splitting, we
expect the elastic two-nucleon intermediate state to provide the dominant
contribution. While we will mostly focus on the elastic channel, we will also
estimate the effect of the leading $N\\!N\pi$ inelastic channel as we expect
this to be one of the dominant sources of uncertainty in our final result.

##  3 Chiral EFT result to leading order

In this section we briefly revisit the chiral EFT result of Refs. [40, 43] in
light of the representation given in Eq. (2.23). This will serve two purposes:
setting up the notation and pointing to a useful way of organizing the
integrand in the full theory amplitude.

The power counting for ${\cal A}_{\nu}$ in chiral EFT is described in Refs.
[40, 43] and we recall here some of its elements as needed. Denoting by
$\mu_{\chi}\sim Q$ the soft scale in the EFT, to LO in chiral counting, i.e.,
$1/Q^{2}$, only elastic $N\\!N$ intermediate states are relevant. The
corresponding diagrams are reported in Fig. 2. For concreteness, we regulate
all the integrals dimensionally and perform $\overline{\rm MS}$ subtraction of
the divergences when needed. The LO chiral EFT results correspond to replacing
in Eq. (2.23) the LO form for the $\hat{O}^{LL}({\mathbf{k}})$ and
$\hat{T}(E)$ operators, denoted by $\hat{O}_{\chi}^{LL}({\mathbf{k}})$ and
$\hat{T}^{\chi}(E)$, respectively, and using the non-relativistic form of the
free two-nucleon Green’s function,
$[G^{(0)}_{+}(E)]_{n}=1/(E-{\mathbf{p}}_{n}^{2}/m_{N}+i\epsilon)$.

The current–current correlator to LO in chiral EFT is given by

${}_{{}^{1}S_{0}}\langle{\mathbf{p}}^{\prime}|\hat{O}_{\chi}^{LL}({\mathbf{k}})|{\mathbf{p}}\rangle_{{}^{1}S_{0}}\
=\
-\frac{1+2g_{A}^{2}+\frac{g_{A}^{2}M_{\pi}^{4}}{({\mathbf{k}}^{2}+M_{\pi}^{2})^{2}}}{{\mathbf{k}}^{2}}\
\frac{(2\pi)^{3}}{2}\left[\delta^{(3)}({\mathbf{k}}+{\mathbf{p}}-{\mathbf{p}}^{\prime})+\delta^{(3)}({\mathbf{k}}-{\mathbf{p}}+{\mathbf{p}}^{\prime})\right]\,.\qquad$
(3.1)

The $\hat{T}$ operator is determined by the LO interaction Hamiltonian
$\hat{V}=\hat{V}_{\pi}+\hat{V}_{S}$, which contains the one-pion-exchange and
a short-range contribution, parameterized to LO by the LEC $C$:

$\displaystyle\langle{\mathbf{p}}^{\prime}|\hat{V}|{\mathbf{p}}\rangle$
$\displaystyle=V_{S}({\mathbf{p}}^{\prime},{\mathbf{p}})+V_{\pi}({\mathbf{p}}^{\prime},{\mathbf{p}})\,,$
$\displaystyle V_{S}({\mathbf{p}}^{\prime},{\mathbf{p}})$ $\displaystyle=C\,,$
$\displaystyle V_{\pi}({\mathbf{p}}^{\prime},{\mathbf{p}})$
$\displaystyle=-\frac{4\pi\alpha_{\pi}}{({\mathbf{p}}-{\mathbf{p}}^{\prime})^{2}+M_{\pi}^{2}}\,,\qquad\alpha_{\pi}=\frac{g_{A}^{2}M_{\pi}^{2}}{16\pi
F_{\pi}^{2}}\,.$ (3.2)

Here $F_{\pi}$ is the pion decay constant and $M_{\pi}$ denotes the pion mass.
The LO ${}^{1}S_{0}$ $N\\!N$ contact coupling $C$ contains contributions from
one-pion exchange as well as a contact interaction and scales as $C\sim
4\pi/(m_{N}Q)$. The split $\hat{V}=\hat{V}_{\pi}+\hat{V}_{S}$ implies that the
$\hat{T}$ operator can be similarly separated into a pion-range and short-
range contribution as follows [80]

$\displaystyle\hat{T}(E)$ $\displaystyle=\hat{T}_{\pi}(E)+\hat{T}_{S}(E)\,,$
$\displaystyle\hat{T}_{\pi}(E)$
$\displaystyle=\hat{V}_{\pi}\left(I+\hat{G}_{+}^{(\pi)}(E)\,\hat{V}_{\pi}\right)=\hat{V}_{\pi}\
\left(I-\hat{G}^{(0)}_{+}(E)\,\hat{V}_{\pi}\right)^{-1}\,,$
$\displaystyle\hat{T}_{S}(E)$
$\displaystyle=\left(I+\hat{V}_{\pi}\hat{G}_{+}^{(\pi)}(E)\right)\left[\hat{V}_{S}\,\left(I-\hat{G}^{(\pi)}_{+}(E)\,\hat{V}_{S}\right)^{-1}\,\right]\left(I+\hat{G}_{+}^{(\pi)}(E)\,\hat{V}_{\pi}\right)\,,$
(3.3)

where $\hat{G}_{+}^{(\pi)}(E)=1/(E-\hat{H}_{0}-\hat{V}_{\pi}+i\epsilon)$ is
the Green’s function associated with the pion-exchange interaction.

Using Eq. (2.23) and separating out the contributions with zero, one, and two
insertions of $\hat{T}_{S}^{\chi}$, the LO chiral EFT amplitude can be written
as

${\cal A}_{\nu}^{\chi{\rm EFT}}={\cal M}_{A}^{\chi}+{\cal
M}_{B+\bar{B}}^{\chi}+{\cal M}_{C}^{\chi}\,,$ (3.4)

with

$\displaystyle{\cal M}_{A}^{\chi}$
$\displaystyle=\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\ \langle f_{0}|\
\left(\hat{T}^{\chi}_{\pi}(E^{\prime})\hat{G}_{+}^{(0)}(E^{\prime})+I\right)\,\hat{O}_{\chi}^{LL}({\mathbf{k}})\,\left(I+\hat{G}_{+}^{(0)}(E)\,\hat{T}^{\chi}_{\pi}(E)\right)|i_{0}\rangle\,,$
$\displaystyle{\cal M}_{B+\bar{B}}^{\chi}$
$\displaystyle=\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\ \langle
f_{0}|\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)|i_{0}\rangle$
$\displaystyle+\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\ \langle
f_{0}|\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)|i_{0}\rangle\,,$
$\displaystyle{\cal M}_{C}^{\chi}$
$\displaystyle=\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\ \langle
f_{0}|\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)|i_{0}\rangle\,.$
(3.5)

Diagrammatically, ${\cal M}_{A}^{\chi}$, ${\cal M}_{B+\bar{B}}^{\chi}$, and
${\cal M}_{C}^{\chi}$ correspond to the first, second, and third row in Fig.
2. Using the definition of $\hat{T}_{S}$ in Eq. (3.3) and the fact that to LO
$\hat{V}_{S}$ is a momentum-independent contact interaction, in dimensional
regularization one can show that the rescattering factors
$K_{E},K_{E^{\prime}}$

$K_{E}=\frac{C}{1-CG^{+}_{E}({\mathbf{0}},{\mathbf{0}})}\,,\qquad
G^{\pm}_{E}({\mathbf{r}},{\mathbf{r}}^{\prime})=\int\frac{d^{3}{\mathbf{q}}}{(2\pi)^{3}}\int\frac{d^{3}{\mathbf{q}}^{\prime}}{(2\pi)^{3}}\,e^{i{\mathbf{q}}\cdot
r}e^{-i{\mathbf{q}}^{\prime}\cdot{\mathbf{r}}^{\prime}}\langle{\mathbf{q}}|\hat{G}_{\pm}^{(\pi)}(E)|{\mathbf{q}}^{\prime}\rangle\,,$
(3.6)

and Yukawa wave functions at the origin $\chi_{\mathbf{p}}^{+}({\mathbf{0}})$,
$\chi_{{\mathbf{p}}^{\prime}}^{+}({\mathbf{0}})$

$\chi^{\pm}_{{\mathbf{p}}}({\mathbf{r}})=\int\frac{d^{3}{\mathbf{q}}}{(2\pi)^{3}}\,e^{i{\mathbf{k}}\cdot{\mathbf{r}}}\langle{\mathbf{k}}|(1+\hat{G}^{(\pi)}_{\pm}(E)\,\hat{V}_{\pi})|{\mathbf{p}}\rangle\,,\\\
$ (3.7)

can be factored out of the $d^{3}{\mathbf{k}}$ integrals in ${\cal
M}_{B+\bar{B}}^{\chi}$, and ${\cal M}_{C}^{\chi}$, thus reproducing the chiral
EFT results of Refs. [40, 43]

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

The divergence in $G^{+}_{E}({\mathbf{0}},{\mathbf{0}})$ is absorbed by
$C^{-1}$, so that $K_{E}$ is well defined and independent of the chosen scheme
and scale [23].

Figure 2: Diagrammatic representation of LO contributions to $nn\rightarrow
ppee$ in chiral EFT. Double, dashed, and plain lines denote nucleons, pions,
and leptons, respectively. The blue ellipse represents iteration of
$\hat{V}_{\pi}$. Gray circles denote the nucleon axial and vector currents. In
the first three lines the black square represents an insertion of
$m_{\beta\beta}$, while in the fourth line it represents an insertion of
${\cal C}_{1}=g^{NN}_{\nu}$. The ellipses in the second to fourth lines denote
diagrams with an arbitrary number of $N\\!N$ bubble insertions.

The reduced amplitudes ${\cal A}_{A,B,C}$ correspond to the left-most diagram
in the first, second, and third rows of Fig. 2 respectively (without the
Yukawa iteration in the outer legs after the $\hat{V}_{S}$ insertion).
Neglecting Yukawa interactions ($\alpha_{\pi}\to 0$), i.e., in the pionless
EFT (${/\pi}$EFT), ${\cal A}_{A,B,C}$ have simple expressions [40, 43]. In
general, however, they are non-perturbative objects that involve the sum of
infinitely many Feynman diagrams. The amplitude ${\cal A}_{C}$ in Eqs. (3.5)
and (3.8) contains a UV-divergent term at the two-loop level, which we denote
by ${\cal A}_{C}^{\rm sing}$, as well as UV-finite terms induced by pion
exchange, which we denote by $\delta{\cal A}_{C}$, leading to the
decomposition:

${\cal A}_{C}={\cal A}_{C}^{\rm sing}+\delta{\cal A}_{C}\,.$ (3.9)

The UV-convergent term $\delta{\cal A}_{C}$ arises from (i) the iteration of
the pion-induced potential (see Fig. 2 and Eq. (3.5), as well as Eq. (3.3) for
the definition of $\hat{T}_{S}^{\chi}$); (ii) the term proportional to
$g_{A}^{2}M_{\pi}^{4}$ in $\hat{O}^{LL}_{\chi}({\mathbf{k}})$, see Eq. (3.1),
which is one of two manifestations of the induced pseudoscalar form factor of
the axial current (the other is the change $3g_{A}^{2}\to 2g_{A}^{2}$). Using
dimensional regularization with scale $\mu_{\chi}$ and the $\overline{\rm MS}$
scheme for renormalization, we thus identify the singular (UV-divergent) term
with

${\cal A}_{C}^{\rm sing}(\mu_{\chi})=-m_{N}^{2}\
\mu_{\chi}^{4-d}\,\int\frac{d^{d-1}k}{(2\pi)^{d-1}}\,\frac{1+2g_{A}^{2}}{{\bf
k}^{2}-i\epsilon}\ {\cal I}_{C}({\bf k}^{2},{\bf p}^{2},{\bf p}^{\prime
2})\,,$ (3.10)

where

$\displaystyle{\cal I}_{C}({\bf k}^{2},{\bf p}^{2},{\bf p}^{\prime 2})$
$\displaystyle=\int\frac{d^{3}{\mathbf{q}}}{(2\pi)^{3}}\,\frac{1}{{\bf
p}^{\prime 2}-({\bf q}+{\bf k})^{2}+i\epsilon}\ \frac{1}{{\bf p}^{2}-{\bf
q}^{2}+i\epsilon}\,,$ $\displaystyle{\cal I}_{C}({\bf k}^{2},{\bf p}^{2},{\bf
p}^{2})$ $\displaystyle=\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|\,.$ (3.11)

In case of $\bf p=\bf p^{\prime}$ the resulting expression becomes

${\cal A}_{C}^{\rm
sing}(\mu_{\chi})=-\frac{(1+2g_{A}^{2})}{2}\,\frac{m_{N}^{2}}{(4\pi)^{2}}\left[\log\frac{\mu_{\chi}^{2}}{-4|{\bf
p}|^{2}-i\epsilon}+1\right]\,.$ (3.12)

It will prove useful to write the real part of this result in an integral
representation

$\mathrm{Re}\,{\cal A}_{C}^{\rm
sing}(\mu_{\chi})=-\frac{(1+2g_{A}^{2})}{2}\,\frac{m_{N}^{2}}{(4\pi)^{2}}\left[1+2\int_{0}^{\mu_{\chi}}\frac{d{|\bf
k}|}{|{\bf k}|}\theta(|{\bf k}|-2|{\bf p}|)\right]\,.$ (3.13)

The UV divergence is removed by introducing a new contact term

$\mathcal{L}_{|\Delta L|=2}^{NN}=\left(2\,G_{F}^{2}V_{ud}^{2}\right)\
m_{\beta\beta}\,{\cal
C}_{1}\times\bar{e}_{L}e_{L}^{c}\,\bar{N}\mathcal{\tau}^{+}N\,\bar{N}\tau^{+}N\,,$
(3.14)

with a LNV coupling proportional to $m_{\beta\beta}$ denoted by ${\cal C}_{1}$
[43], to be identified with $g_{\nu}^{NN}$ also used in the recent literature,
i.e., ${\cal C}_{1}\equiv g_{\nu}^{NN}$. A chiral-covariant form of the
contact operator in Eq. (3.14) will be given in Section 6.

Including the contact term, the finite, renormalized amplitude is given by

$\displaystyle{\cal A}_{\nu}^{\chi{\rm EFT}}$
$\displaystyle=\mathcal{A}_{A}+\chi^{+}_{{\mathbf{p}}^{\prime}}({\mathbf{0}})\,K_{E^{\prime}}\,\mathcal{A}_{B}+\mathcal{\bar{A}}_{B}\,K_{E}\,\chi^{+}_{{\mathbf{p}}}({\mathbf{0}})$
$\displaystyle+\chi^{+}_{{\mathbf{p}}^{\prime}}({\mathbf{0}})\,K_{E^{\prime}}\,\left(\mathcal{A}_{C}^{\rm
sing}(\mu_{\chi})+\frac{2\,{\cal C}_{1}(\mu_{\chi})}{C^{2}}+\delta{\cal
A}_{C}\right)\,K_{E}\,\chi^{+}_{{\mathbf{p}}}({\mathbf{0}})\,.$ (3.15)

The matching analysis will provide a representation for the combination ${\cal
A}_{C}^{\rm sing}(\mu_{\chi})+2\,{\cal C}_{1}(\mu_{\chi})/C^{2}$.

##  4 Full theory parameterization

The starting point for parameterizing the full-theory amplitude is provided by
Eq. (2.23), in which we split the integral into a low- plus intermediate-
momentum region and high-momentum region, according to Eq. (2.12). In order to
mimic the structure of the chiral EFT amplitude and facilitate matching, it is
also convenient to separate the full theory $T$ matrix into a term induced by
one-pion exchange and a short-range contribution, i.e.,
$\hat{T}=\hat{T}_{\pi}+\hat{T}_{S}$, as done in Eq. (3.3).

###  4.1 Low- and intermediate-momentum region: $\boldsymbol{{\cal A}^{<}}$

Following Eqs. (2.23) and (3.4), we write

${\cal A}^{<}={\cal M}_{A}^{<}+{\cal M}_{B+\bar{B}}^{<}+{\cal M}_{C}^{<}\,,$
(4.1)

with

$\displaystyle{\cal M}^{<}_{A}$
$\displaystyle=\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\ \langle f_{0}|\
\left(\hat{T}^{<}_{\pi}(E^{\prime})\hat{G}_{+}^{(0)}(E^{\prime})+I\right)\,\hat{O}_{<}^{LL}({\mathbf{k}})\,\left(I+\hat{G}_{+}^{(0)}(E)\,\hat{T}^{<}_{\pi}(E)\right)|i_{0}\rangle\,,$
$\displaystyle{\cal M}_{B+\bar{B}}^{<}$
$\displaystyle=\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\ \langle
f_{0}|\left(\hat{T}^{<}_{S}(E^{\prime})\,\hat{G}_{+}^{(0)}(E^{\prime})\right)\,\hat{O}_{<}^{LL}({\mathbf{k}})\,\left(I+\hat{G}_{+}^{(0)}(E)\,\hat{T}^{<}_{\pi}(E)\right)|i_{0}\rangle$
$\displaystyle+\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\ \langle
f_{0}|\left(I+\hat{T}^{<}_{\pi}(E^{\prime})\hat{G}_{+}^{(0)}(E)\right)\,\hat{O}_{<}^{LL}({\mathbf{k}})\,\left(\hat{G}_{+}^{(0)}(E)\hat{T}^{<}_{S}(E)\right)|i_{0}\rangle\,,$
$\displaystyle{\cal M}_{C}^{<}$
$\displaystyle=\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\ \langle
f_{0}|\left(\hat{T}^{<}_{S}(E^{\prime})\hat{G}_{+}^{(0)}(E^{\prime})\right)\,\hat{O}_{<}^{LL}({\mathbf{k}})\,\left(\hat{G}_{+}^{(0)}(E)\,\hat{T}^{<}_{S}(E)\right)|i_{0}\rangle\,.$
(4.2)

We wish to provide a representation of the integrand in ${\cal M}^{<}_{A,B,C}$
that is valid up to a scale $|{\mathbf{k}}|\sim\Lambda$. We build the
integrand starting from the low-momentum region. For $|{\mathbf{k}}|\sim Q$
the integrand is controlled by chiral EFT. Using the LO chiral representation
for the two-current operator ($\hat{O}_{<}^{LL}=\hat{O}_{\chi}^{LL}$) and the
$T$ matrix ($\hat{T}^{<}_{\pi,S}=\hat{T}^{\chi}_{\pi,S}$) in ${\cal M}_{A}$
and ${\cal M}_{B+\bar{B}}$ leads to convergent integrals, and one recovers the
chiral EFT results of Eqs. (3.5), up to terms of $\mathcal{O}(Q/\Lambda)$ that
are irrelevant at LO

$\displaystyle{\cal M}_{A}^{<}$ $\displaystyle={\cal
M}_{A}^{\chi}\times\left(1+\mathcal{O}\left(\frac{Q}{\Lambda}\right)\right)\,,$
$\displaystyle{\cal M}_{B+\bar{B}}^{<}$ $\displaystyle={\cal
M}_{B+\bar{B}}^{\chi}\times\left(1+\mathcal{O}\left(\frac{Q}{\Lambda}\right)\right)\,.$
(4.3)

We will also show that ${\cal M}^{>}_{A,B+\bar{B}}$ does not contribute to the
amplitude at LO. Therefore, ${\cal M}_{A}$ and ${\cal M}_{B+\bar{B}}$ drop out
of the LO matching condition (2.13), which will therefore only involve ${\cal
M}_{C}^{<}$. In order to have a representation of the integrand valid up to
momenta $|{\mathbf{k}}|\sim\Lambda$, in Eq. (4.2) we need to model both
$\hat{O}^{LL}({\mathbf{k}})$ and $\hat{T}_{S}$ beyond the LO chiral EFT
expressions. We will do so by introducing appropriate single-nucleon weak form
factors and HOS form factors, represented by the red and blue circles in Fig.
3, respectively.

Figure 3:  Diagrammatic representation of the forward Compton amplitude
corresponding to the integrand of ${\cal A}_{C}^{\rm sing}$ in Eq. (3.10) and
${\cal A}_{C}^{<,\rm sing}$ in Eq. (4.18). Solid lines represent nucleons and
wavy lines the currents. Blue circles represent short-range $N\\!N$
interactions and red circles represent single-nucleon weak form factors.

For the current–current operator we make the replacement
$\hat{O}^{LL}_{\chi}({\mathbf{k}})\to\hat{O}^{LL}_{<}({\mathbf{k}})$, by
introducing single-nucleon form factors in the weak current vertices,

${}_{{}^{1}S_{0}}\langle{\mathbf{p}}^{\prime}|\hat{O}_{<}^{LL}({\mathbf{k}})|{\mathbf{p}}\rangle_{{}^{1}S_{0}}\
=\
-\frac{g_{V}^{2}({\mathbf{k}}^{2})+3h_{GT}^{LL}({\mathbf{k}}^{2})}{{\mathbf{k}}^{2}}\,\frac{(2\pi)^{3}}{2}\left[\delta^{(3)}({\mathbf{k}}+{\mathbf{p}}-{\mathbf{p}}^{\prime})+\delta^{(3)}({\mathbf{k}}-{\mathbf{p}}+{\mathbf{p}}^{\prime})\right]\,,$
(4.4)

where

$h^{LL}_{GT}({\mathbf{k}}^{2})=g_{A}^{2}({\mathbf{k}}^{2})+\frac{{\mathbf{k}}^{4}g_{P}^{2}({\mathbf{k}}^{2})}{12m_{N}^{2}}+\frac{{\mathbf{k}}^{2}g_{A}({\mathbf{k}}^{2})g_{P}({\mathbf{k}}^{2})}{3m_{N}}+\frac{{\mathbf{k}}^{2}g_{M}^{2}({\mathbf{k}}^{2})}{6m_{N}^{2}}\,.$
(4.5)

In our baseline analysis, we will use the simple dipole parameterization for
the vector, magnetic, and axial-vector form factors

$\displaystyle g_{V}({\mathbf{k}}^{2})$
$\displaystyle=\frac{1}{\left(1+\frac{{\mathbf{k}}^{2}}{\Lambda_{V}^{2}}\right)^{2}}\,,\qquad
g_{M}({\mathbf{k}}^{2})=(1+\kappa_{1})\,g_{V}({\mathbf{k}}^{2})\,,\qquad\kappa_{1}\simeq
3.7\,,$ $\displaystyle g_{A}({\mathbf{k}}^{2})$
$\displaystyle=\frac{g_{A}}{\left(1+\frac{{\mathbf{k}}^{2}}{\Lambda_{A}^{2}}\right)^{2}}\,,\qquad
g_{P}({\mathbf{k}}^{2})=-\frac{2m_{N}\,g_{A}({\mathbf{k}}^{2})}{{\mathbf{k}}^{2}+M_{\pi}^{2}}\,,$
(4.6)

with $\Lambda_{V}=0.84\,\text{GeV}$ and $\Lambda_{A}=1.0\,\text{GeV}$. While
we expect this to be a good first approximation, we have explored the
dependence of the contact term on the form factor input. In the case of the
axial-vector form factor, we have taken the range
$\Lambda_{A}\in[0.8,1.2]\,\text{GeV}$, which corresponds to an axial radius of
$0.47(19)\,\text{fm}^{2}$, in good agreement with the result
$0.46(16)\,\text{fm}^{2}$ quoted in Ref. [81]. For the vector form factor, we
have explored the continued fraction expansion of Ref. [82]. Since input on
the form factors turns out to induce a by far subdominant uncertainty in the
matching result (see discussion below), we will not perform a more
sophisticated error analysis based on state-of-the-art information on the
nucleon form factors [83, 84, 85, 81].

The parameterization of $\hat{O}^{LL}_{<}({\mathbf{k}})$ in terms of the weak
nucleon form factors then leads to

$\displaystyle{\cal M}_{C}^{<}$
$\displaystyle=-\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\
\left[\frac{g_{V}^{2}({\mathbf{k}}^{2})+3h_{GT}^{LL}({\mathbf{k}}^{2})}{{\mathbf{k}}^{2}}\,\right]$
$\displaystyle\times\int\frac{d^{3}{\mathbf{q}}}{(2\pi)^{3}}\
T_{S}({\mathbf{p}}^{\prime},{\mathbf{q}}+{\mathbf{k}})\
\frac{1}{E^{\prime}-({\mathbf{q}}+{\mathbf{k}})^{2}/m_{N}+i\epsilon}\,\frac{1}{E-{\mathbf{q}}^{2}/m_{N}+i\epsilon}\
T_{S}({\mathbf{q}},{\mathbf{p}})\,,$ (4.7)

where the notation for the HOS $\hat{T}_{S}$-matrix elements is666These are
$S$-wave projected matrix elements and depend only on $|{\mathbf{p}}|$ and
$|{\mathbf{q}}|$.

$\displaystyle T_{S}({\mathbf{q}},{\mathbf{p}})$
$\displaystyle=\langle{\mathbf{q}}|\hat{T}_{S}(E)|{\mathbf{p}}\rangle\,,$
$\displaystyle E$ $\displaystyle=\frac{{\mathbf{p}}^{2}}{m_{N}}\,,$
$\displaystyle T_{S}({\mathbf{p}}^{\prime},{\mathbf{q}})$
$\displaystyle=\langle{\mathbf{p}}^{\prime}|\hat{T}_{S}(E^{\prime})|{\mathbf{q}}\rangle\,,$
$\displaystyle E^{\prime}$ $\displaystyle=\frac{{\mathbf{p}}^{\prime
2}}{m_{N}}\,.$ (4.8)

We further write the HOS $\hat{T}_{S}$-matrix element as the product of the
on-shell $\hat{T}_{S}$-matrix element times the HOS factors
$\tilde{f}_{S}({\mathbf{q}},{\mathbf{p}})$ [86]

$T_{S}({\mathbf{q}},{\mathbf{p}})\equiv
T_{S}({\mathbf{p}},{\mathbf{p}})\times\tilde{f}_{S}({\mathbf{q}},{\mathbf{p}})\,.$
(4.9)

Using the LO on-shell result [23], see Eq. (3.3),

$T_{S}({\mathbf{p}},{\mathbf{p}})=\left(\chi^{+}_{\mathbf{p}}({\mathbf{0}})\right)^{2}\
K_{E}$ (4.10)

and defining

$\bar{f}_{S}({\mathbf{q}},{\mathbf{p}})\equiv\tilde{f}_{S}({\mathbf{q}},{\mathbf{p}})\
\chi^{+}_{\mathbf{p}}({\mathbf{0}})=\frac{T_{S}({\mathbf{q}},{\mathbf{p}})}{T_{S}({\mathbf{p}},{\mathbf{p}})}\
\chi^{+}_{\mathbf{p}}({\mathbf{0}})\,,$ (4.11)

we arrive at

$\displaystyle{\cal M}_{C}^{<}$
$\displaystyle=\chi^{+}_{{\mathbf{p}}^{\prime}}({\mathbf{0}})\,K_{E^{\prime}}\,\Big{[}\,{\cal
A}^{<}_{C}\,\Big{]}\,K_{E}\,\chi^{+}_{{\mathbf{p}}}({\mathbf{0}})\,,$ (4.12)
$\displaystyle{\cal A}_{C}^{<}$
$\displaystyle=-\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\
\left[\frac{g_{V}^{2}({\mathbf{k}}^{2})+3h_{GT}^{LL}({\mathbf{k}}^{2})}{{\mathbf{k}}^{2}}\,\right]\,\times\,\tilde{\cal
I}_{C}^{<}(|{\mathbf{k}}|)\,,$ (4.13) $\displaystyle\tilde{\cal
I}_{C}^{<}(|{\mathbf{k}}|)$
$\displaystyle=\int\frac{d^{3}{\mathbf{q}}}{(2\pi)^{3}}\
\bar{f}_{S}({\mathbf{p}}^{\prime},{\mathbf{q}}+{\mathbf{k}})\
\frac{1}{E^{\prime}-({\mathbf{q}}+{\mathbf{k}})^{2}/m_{N}+i\epsilon}\,\frac{1}{E-{\mathbf{q}}^{2}/m_{N}+i\epsilon}\
\bar{f}_{S}({\mathbf{q}},{\mathbf{p}})\,.$ (4.14)

Equation (4.12) shows that ${\cal M}_{C}^{<}$ has the correct IR behavior to
LO in the external momenta ${\mathbf{p}}$ and ${\mathbf{p}}^{\prime}$ and that
the quantity ${\cal A}_{C}^{<}$ corresponds to the amplitude ${\cal A}_{C}$ in
chiral EFT. As argued in Section 3, ${\cal A}_{C}$ contains a singular UV-
divergent term (${\cal A}_{C}^{\rm sing}$) and UV-finite contributions induced
by pion exchange ($\delta{\cal A}_{C}$). In analogy to Eq. (3.9) we can write

${\cal A}_{C}^{<}={\cal A}_{C}^{<,\rm sing}+\delta{\cal A}_{C}^{<}\,,$ (4.15)

where $\delta{\cal A}_{C}^{<}$ denotes the convergent pion-exchange
contributions to ${\cal A}_{C}^{<}$, satisfying

$\delta{\cal A}^{<}_{C}-\delta{\cal A}_{C}\sim{\cal
A}_{C}\times\mathcal{O}\left(\frac{Q}{\Lambda_{\chi}}\right)\,.$ (4.16)

Therefore, these finite contributions drop out of the matching relation (2.13)
to LO and for the purpose of matching one only needs to identify the singular
component ${\cal A}_{C}^{<}$.

${\cal A}_{C}^{<,\rm sing}$ is obtained from Eq. (4.13) by systematically
removing the convergent pion contributions. This requires:

  1. 1.

Discarding the convergent terms due to the induced pseudoscalar form factor in
$h_{GT}^{LL}({\mathbf{k}}^{2})$. In practice this means replacing
$3h_{GT}^{LL}({\mathbf{k}}^{2})\to
2g_{A}({\mathbf{k}}^{2})+\frac{{\mathbf{k}}^{2}g_{M}^{2}({\mathbf{k}}^{2})}{2m_{N}^{2}}$.

  2. 2.

Setting $\alpha_{\pi}\to 0$ in the evaluation of
$\bar{f}_{S}({\mathbf{q}},{\mathbf{p}})$ in Eq. (4.14), which implies

$\bar{f}_{S}({\mathbf{q}},{\mathbf{p}})\bigg{|}_{\alpha_{\pi}=0}=\frac{T_{S}({\mathbf{q}},{\mathbf{p}})}{T_{S}({\mathbf{p}},{\mathbf{p}})}\bigg{|}_{\alpha_{\pi}=0}\equiv
f_{S}({\mathbf{q}},{\mathbf{p}})\,.$ (4.17)

The quantity $f_{S}({\mathbf{q}},{\mathbf{p}})$ is real-valued and can be
interpreted as a form factor. In light of this, from Eqs. (4.12)–(4.14) we
arrive at

$\displaystyle{\cal A}_{C}^{<,{\rm sing}}$
$\displaystyle=-m_{N}^{2}\,\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\
\left[\frac{g_{V}^{2}({\mathbf{k}}^{2})+2g^{2}_{A}({\mathbf{k}}^{2})+({\mathbf{k}}^{2}g_{M}^{2}({\mathbf{k}}^{2}))/(2m_{N}^{2})}{{\mathbf{k}}^{2}}\right]\times{\cal
I}_{C}^{<}(|{\mathbf{k}}|)\,,$ $\displaystyle{\cal I}_{C}^{<}(|{\mathbf{k}}|)$
$\displaystyle=\int\frac{d^{3}{\mathbf{q}}}{(2\pi)^{3}}\
f_{S}({\mathbf{p}}^{\prime},{\mathbf{q}}+{\mathbf{k}})\ \frac{1}{{\bf
p}^{\prime 2}-({\bf q}+{\bf k})^{2}+i\epsilon}\ \frac{1}{{\bf p}^{2}-{\bf
q}^{2}+i\epsilon}~{}\ f_{S}({\mathbf{q}},{\mathbf{p}})\,.$ (4.18)

In the end, the difference between ${\cal A}_{C}^{\rm sing}(\mu_{\chi})$ and
${\cal A}_{C}^{<,{\rm sing}}$ is quite intuitive. As illustrated in Fig. 3,
the integrand of ${\cal A}_{C}^{\rm sing}$ in Eq. (3.10) contains the LO
nucleon weak current vertices and the leading contact $N\\!N$ interaction. On
the other hand, the full-theory integrand in ${\cal A}_{C}^{<,{\rm sing}}$
contains appropriate form factors in the nucleon weak current vertices and in
the $N\\!N$ vertices ($f_{S}({\mathbf{q}},{\mathbf{p}})$). The latter
parameterizes the off-shell behavior of the short-range component of the
$N\\!N$ amplitude and amounts to changing the LO chiral EFT function ${\cal
I}_{C}(|{\mathbf{k}}|)$ defined in Eq. (3.11) to ${\cal
I}_{C}^{<}(|{\mathbf{k}}|)$ given in Eq. (4.18), valid for a larger range of
momenta $|{\mathbf{k}}|$ compared to the LO result. It will prove convenient
to present results in terms of the ratio

$r(|{\mathbf{k}}|)\equiv\frac{{\rm Re}\,{\cal I}_{C}^{<}(|{\mathbf{k}}|)}{{\rm
Re}\,{\cal I}_{C}(|{\mathbf{k}}|)}\,.$ (4.19)

Figure 4:  Ratio $r(|{\mathbf{k}}|)$ defined in Eq. (4.19). The various lines
correspond to the NLO ${/\pi}$EFT result (blue), the NLO chiral EFT result
(red), and the three-Yukawa potential from Ref. [87] (green), compared to the
LO chiral EFT result $r(|{\mathbf{k}}|)=1$.

We now describe our construction of ${\cal I}_{C}^{<}(|{\mathbf{k}}|)$ and
$r(|{\mathbf{k}}|)$, based on the first two orders in chiral EFT and then
extended to larger values of $|{\mathbf{k}}|$ using models for the short-range
$N\\!N$ interaction in the ${}^{1}S_{0}$ channel. For concreteness we work at
the kinematic point $|{\mathbf{p}}|=|{\mathbf{p}}^{\prime}|=1\,\text{MeV}$.

  1. 1.

At LO in chiral EFT one has $f_{S}({\mathbf{q}},{\mathbf{p}})=1$ and the
integrals simplify to
${\rm Re}\,{\cal I}_{C}^{<}(|{\mathbf{k}}|)\to{\rm Re}\,{\cal
I}_{C}(|{\mathbf{k}}|)=\theta(|{\mathbf{k}}|-2|{\mathbf{p}}|)/(8|{\mathbf{k}}|)$.

  2. 2.

For $|{\mathbf{k}}|\leq M_{\pi}$ and $|{\mathbf{k}}|\leq\Lambda_{\chi}$ we can
evaluate the corrections to ${\rm Re}\ {\cal I}_{C}^{<}(|{\mathbf{k}}|)$ in
${/\pi}$EFT and chiral EFT, respectively, using
$V_{S}({\mathbf{p}},{\mathbf{p}}^{\prime})=C+(C_{2}/2)({\mathbf{p}}^{2}+{\mathbf{p}}^{\prime
2})$. As discussed in Appendix A, the results at next-to-leading order (NLO)
in chiral EFT and ${/\pi}$EFT are formally identical

$\displaystyle f_{S}({\mathbf{q}},{\mathbf{p}})$
$\displaystyle=1-\frac{C_{2}}{2C}({\mathbf{p}}^{2}-{\mathbf{q}}^{2})\,,$
$\displaystyle f_{S}({\mathbf{p}}^{\prime},{\mathbf{q}}+{\mathbf{k}})$
$\displaystyle=1-\frac{C_{2}}{2C}\left[{\mathbf{p}}^{\prime\,2}-({\mathbf{q}}+{\mathbf{k}})^{2}\right]\,,$
(4.20)

with the identifications $C=C_{0}$ in ${/\pi}$EFT and
$C=C_{0}+g_{A}^{2}/(4F_{\pi}^{2})$ in chiral EFT. Note that these HOS form
factors are different depending on whether they involve the initial- or final-
state on-shell momenta, so that $f_{S}({\mathbf{p}}^{\prime},{\mathbf{q}})$
cannot be obtained from the expression of $f_{S}({\mathbf{q}},{\mathbf{p}})$
by interchanging ${\mathbf{q}}\leftrightarrow{\mathbf{p}}$ and subsequently
replacing ${\mathbf{p}}\to{\mathbf{p}^{\prime}}$. This leads to (see Appendix
A for details)

${\rm Re}\ {\cal
I}_{C}^{<}(|{\mathbf{k}}|)=\frac{1}{8|{\mathbf{k}}|}-\frac{C_{2}}{m_{N}C^{2}}\,.$
(4.21)

In terms of the ratio $r(|{\mathbf{k}}|)$ introduced in Eq. (4.19) the NLO
analysis gives

$r^{\rm NLO}(|{\mathbf{k}}|)=1-\frac{8C_{2}}{m_{N}C^{2}}|{\mathbf{k}}|\,.$
(4.22)

In both EFTs the ratio $C_{2}/C^{2}$ is renormalization-scale independent. In
${/\pi}$EFT $C_{2}/C^{2}=m_{N}r_{0}/(8\pi)$ is linked to the effective range,
leading to

$r^{\rm NLO}_{/\pi}(|{\mathbf{k}}|)=1-\frac{r_{0}}{\pi}|{\mathbf{k}}|\,.$
(4.23)

Since $r_{0}\simeq 1/(72\,\text{MeV})$ (corresponding to $C_{2}/C^{2}\simeq
0.52$), this produces sizable fractional deviations in ${\rm Re}\ {\cal
I}_{C}^{<}(|{\mathbf{k}}|)$ already for $|{\mathbf{k}}|\leq M_{\pi}$. In
chiral EFT one finds a slightly reduced suppression $C_{2}/C^{2}\simeq 0.38$
[23], as part of the effective range is already captured by the pion-exchange
contribution. In summary, the NLO analysis shows that there is a sizable
suppression of the LO result already for relatively low $|{\mathbf{k}}|$,
linked to the large effective range in the ${}^{1}S_{0}$ channel.

  3. 3.

To extend ${\cal I}_{C}^{<}(|{\mathbf{k}}|)$ or equivalently
$r(|{\mathbf{k}}|)$ to higher values of $|{\mathbf{k}}|$, we have computed
$f_{S}({\mathbf{q}},{\mathbf{p}})$ using models of $V_{S}$, i.e., the short-
range $N\\!N$ interaction in the ${}^{1}S_{0}$ channel. We have used the
three-Yukawa potential from Ref. [87, 88] and the AV18 potential [89], all of
which reproduce the ${}^{1}S_{0}$ phase shifts for $|{\mathbf{p}}|$ up to
several hundred MeV. Using the Kaplan–Steel three-Yukawa potential [87] as our
baseline model, we have solved numerically the Lippmann–Schwinger (LS)
equation to obtain $f_{S}({\mathbf{q}},{\mathbf{p}})$, which was then used to
evaluate ${\cal I}_{C}^{<}(|{\mathbf{k}}|)$ (see Appendix A for details).
Reassuringly, as illustrated in Fig. 4, the behavior of
$r_{V_{S}}(|{\mathbf{k}}|)$ in the model calculation very closely tracks the
NLO chiral EFT result for $|{\mathbf{k}}|<200\,\text{MeV}$, i.e., the region
in which NLO chiral EFT is expected to provide model-independent and accurate
results. Using other $N\\!N$ potentials [88, 89] does not change the
qualitative picture, but induces small changes in the slope of
$r(|{\mathbf{k}}|)$ for $|{\mathbf{k}}|<400\,\text{MeV}$ as well as the
location and depth of the minimum for $|{\mathbf{k}}|\sim 600\,\text{MeV}$.
This makes it clear that using $r_{V_{S}}(|{\mathbf{k}}|)$ all the way to
$|{\mathbf{k}}|\sim\Lambda\sim 1\,\text{GeV}$ introduces model dependence.
However, in the integral (4.18) the region in which
$r_{V_{S}}(|{\mathbf{k}}|)$ has the largest model dependence is weighted
considerably less than the model-independent low-$|{\mathbf{k}}|$ region,
because the integration kernel involves $1/|{\mathbf{k}}|$ and the nucleon
form factors, both rapidly decreasing functions of $|{\mathbf{k}}|$. In
Section 5 we will quantify how this model dependence affects the extraction of
the effective coupling.

###  4.2 High-momentum region: $\boldsymbol{{\cal A}^{>}}$

At large Euclidean virtualities ($k_{E}^{2}\gg\Lambda^{2}_{\text{QCD}}$) the
time-ordered product of currents $\hat{\Pi}_{\mu\nu}^{LL}(k,0)$ defined in Eq.
(2.6) admits an OPE. An analysis of the leading diagrams, depicted in Fig. 5,
shows that

$\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}_{\alpha\beta}(0)$
$\displaystyle+g_{\mu\nu}k_{\alpha}k_{\beta}\hat{O}^{\alpha\beta}(0)+k^{2}\hat{O}_{\nu\mu}(0)-k_{\nu}k^{\alpha}\,\hat{O}_{\alpha\mu}(0)-k_{\mu}k^{\alpha}\,\hat{O}_{\nu\alpha}(0)\Big{]}\,,$
(4.24)

where

$\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)=\hat{O}_{\beta\alpha}(0)\,,$ (4.25)

and $T^{a}$ are $SU(3)$ color generators.

Figure 5: Quark-level diagrams that determine the asymptotic behavior of
$\hat{\Pi}^{LL}_{\mu\nu}$. Solid lines represent quarks, curly lines represent
gluons, and the black dots represent insertions of the current
$\bar{u}_{L}\gamma_{\mu}d_{L}$.

In the actual matching calculation, since we contract $T_{\mu\nu}$ with
$g_{\mu\nu}$, we need

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

where the last step holds under symmetric integration ($k_{\alpha}k_{\beta}\to
k^{2}g_{\alpha\beta}/4$). Finally, we note that by using Fierz identities in
color and then in flavor one obtains

$\hat{O}_{\alpha\beta}=\frac{1}{3}\bar{u}_{L}\gamma_{\alpha}d_{L}\
\bar{u}_{L}\gamma_{\beta}d_{L}\ \to\
\hat{O}_{\alpha}^{\alpha}=\frac{1}{3}O_{1}\,,$ (4.27)

where $O_{1}$ is defined in Eq. (2.2).

Performing the $k^{0}$ integration in order to go from
$g^{\mu\nu}\,\hat{\Pi}_{\mu\nu}^{LL}(k,0)$ to $\hat{O}^{LL}({\mathbf{k}})$ we
obtain

$\hat{O}^{LL}_{>}({\mathbf{k}})=\frac{3g_{s}^{2}}{4}\frac{1}{|{\mathbf{k}}|^{5}}\
O_{1}(0)\,,$ (4.28)

which is the expression to be used to evaluate the ${\cal A}^{>}$ component in
Eq. (2.23), to obtain

${\cal A}^{>}=\frac{3\alpha_{s}}{2\pi}\,\langle
f_{-}|\,O_{1}(0)\,|i_{+}\rangle\,\int_{\Lambda}^{\infty}\
d|{\mathbf{k}}|\,\frac{1}{|{\mathbf{k}}|^{3}}\,.$ (4.29)

This result shows a factorization of the hard and soft contributions. The
final step requires the evaluation of $\langle
f_{-}|\,O_{1}(0)\,|i_{+}\rangle$, which can be done in chiral EFT to LO.

As discussed in Refs. [35, 42] the four-quark local operator $O_{1}$ admits a
low-energy realization consistent with chiral symmetry and its breaking. To LO
(and neglecting higher powers of the pion fields) one has

$O_{1}\ \to\ g_{1}^{NN}\,\bar{p}n\bar{p}n\ +\
\frac{5}{6}g_{1}^{\pi\pi}\,F_{\pi}^{2}\,\partial_{\mu}\pi^{-}\partial^{\mu}\pi^{-}\
+\
\sqrt{2}g_{A}g_{1}^{\pi{N}}F_{\pi}\,\bar{p}(S^{\mu}\partial_{\mu}\pi^{-})n+\ldots$
(4.30)

The scaling of the non-perturbative parameters is
$g_{1}^{NN},g_{1}^{\pi{N}}\sim\mathcal{O}(1)$, while a precise lattice
calculation gives $g_{1}^{\pi\pi}=0.36(2)$ at the renormalization scale
$\mu=2\,\text{GeV}$ in the $\overline{\rm MS}$ scheme [90].777See Refs. [91,
92] for a determination of $g_{1}^{\pi\pi}$ based on chiral $SU(3)$,
consistent with the direct lattice result. We have suppressed the dependence
of the short-distance couplings on the QCD renormalization scale $\mu$, needed
to cancel the scale dependence in $\alpha_{s}(\mu)$. The $\pi\pi$ and $\pi N$
couplings induce a pion-range transition operator, while the $N\\!N$ couplings
gives a short-range transition operator. The final result is given by

$\langle f_{-}|\,O_{1}(0)\,|i_{+}\rangle=\tilde{\cal
A}_{A}+\chi^{+}_{{\mathbf{p}}^{\prime}}({\mathbf{0}})\,K_{E^{\prime}}\,\tilde{\cal
A}_{B}+\ \tilde{\cal
A}_{\bar{B}}\,K_{E}\,\chi^{+}_{{\mathbf{p}}}({\mathbf{0}})+\chi^{+}_{{\mathbf{p}}^{\prime}}({\mathbf{0}})\,K_{E^{\prime}}\,\left(\tilde{\cal
A}_{C}+\frac{2\hat{g}_{1}^{NN}}{C^{2}}\right)\,K_{E}\,\chi^{+}_{{\mathbf{p}}}({\mathbf{0}})\,,$
(4.31)

where

$\hat{g}_{1}^{NN}=g_{1}^{NN}+\frac{g_{A}^{2}}{2}\left(\frac{5}{6}g_{1}^{\pi\pi}-g_{1}^{\pi
N}\right)$ (4.32)

and $\tilde{\cal A}_{A,B,C}$ have the same formal expression as ${\cal
A}_{A,B,C}$ [43] with the replacement [42]888Note that the mismatch in
dimensions between $V^{{}^{1}S_{0}}_{\nu\,{\rm L}}({\mathbf{q}})$ and
$V^{{}^{1}S_{0}}_{O_{1}}({\mathbf{q}})$ is compensated by the additional
powers of $1/|{\mathbf{k}}|$ in Eq. (4.29).

$\displaystyle V^{{}^{1}S_{0}}_{\nu\,{\rm L}}({\mathbf{q}})\ \to\
V^{{}^{1}S_{0}}_{O_{1}}({\mathbf{q}})=\
g_{A}^{2}\,\left[\frac{5g_{1}^{\pi\pi}}{6}\,\frac{M_{\pi}^{2}\,{\mathbf{q}}^{2}}{({\mathbf{q}}^{2}+M_{\pi}^{2})^{2}}-\left(g_{1}^{\pi{N}}-\frac{5g_{1}^{\pi\pi}}{6}\right)\,\frac{M_{\pi}^{2}}{{\mathbf{q}}^{2}+M_{\pi}^{2}}\right]\,.$
(4.33)

These relations allow us to write ${\cal A}^{>}$ in a form very useful for the
matching to chiral EFT:

${\cal
A}^{>}=\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}})\,,$
(4.34)

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}}\ \tilde{\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{2\hat{g}_{1}^{NN}}{C^{2}}+\tilde{\cal A}_{C}\right)\,.$ (4.35)

The quantities ${\cal A}_{A,B,\bar{B}}^{>}$ are sub-leading compared to their
finite chiral EFT counterparts ${\cal A}_{A,B,\bar{B}}$, e.g., ${\cal
A}_{A}^{>}/{\cal A}_{A}\sim Q^{2}/\Lambda^{2}$ and similarly for ${\cal
A}_{B}$. Therefore, as anticipated in the previous section, these quantities
do not enter the LO matching relation.

On the other hand, ${\cal A}_{C}^{>}$ has to be retained because its integrand
provides the QCD asymptotic behavior to which the integrand in ${\cal
A}_{C}^{<,\text{sing}}$ given in Eq. (4.18) has to tend for large
$|\mathbf{k}|$. Given our limited knowledge of $g_{1}^{{NN},\pi{N}}$ and the
fact that the terms involving $g_{1}^{{NN}}$, $g_{1}^{\pi{N}}$, and
$g_{1}^{\pi\pi}$ contribute at the same order to $\langle
f_{-}|\,O_{1}(0)\,|i_{+}\rangle$, we will retain for simplicity only the term
proportional to $g_{1}^{NN}$. For convenience, in the matching analysis we
express $g_{1}^{NN}$ as follows

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

In numerical estimates we will allow for $\bar{g}_{1}^{NN}$ as large as
$\mathcal{O}(10)$, given that the $\pi\pi$ short-distance coupling
$\bar{g}_{LR}^{\pi\pi}=8.23$ [90] relevant for the vector–vector correlator by
far exceeds its $\mathcal{O}(1)$ expectation. However, the numerical impact of
these poorly known short-distance parameters on the overall analysis is very
minor, as the main role of ${\cal A}_{C}^{>}$ is to enforce the UV finiteness
of the $nn\to pp$ amplitude.

##  5 Matching and extraction of $\boldsymbol{{\cal C}_{1}}$

###  5.1 Results

The matching condition (2.13), taking into account the form of ${\cal
A}_{\nu}^{\chi{\rm EFT}}$ in Eq. (3.15) and the structure of ${\cal A}^{<}$
(see Eq. (4.1) and following discussion) and ${\cal A}^{>}$ (see Eq. (4.34)
and following discussion) reduces to

${\cal A}_{C}^{<,{\rm sing}}\ +\ {\cal A}_{C}^{>}\ =\ {\cal A}_{C}^{\rm
sing}(\mu_{\chi})\ +\ \frac{2{\cal C}_{1}(\mu_{\chi})}{C^{2}}\,.$ (5.1)

The amplitudes ${\cal A}_{C}^{\rm sing}(\mu_{\chi})$ and ${\cal A}_{C}^{<,{\rm
sing}}$ are given in Eqs. (3.10) and (4.18), respectively, and ${\cal
A}_{C}^{>}$ is given in Eq. (4.35).

Defining the dimensionless coupling $\tilde{\cal C}_{1}$ via

${\cal C}_{1}=\left(\frac{m_{N}}{4\pi}C\right)^{2}\ \tilde{\cal C}_{1}$ (5.2)

and the dimensionless quantities

$\bar{\cal A}_{C}^{X}\equiv\left(\frac{4\pi}{m_{N}}\right)^{2}\ {\cal
A}^{X}_{C}\,,$ (5.3)

the matching condition reads

$\tilde{\cal C}_{1}(\mu_{\chi})=\frac{1}{2}\left[\left(\bar{\cal
A}_{C}^{<,{\rm sing}}\ -\bar{\cal A}_{C}^{\rm
sing}(\mu_{\chi})\right)+\bar{\cal A}_{C}^{>}\right]\,.$ (5.4)

Using our previous results, we can write the various quantities entering the
matching condition in terms of the following integrals:

$\displaystyle\bar{\cal A}^{\rm sing}_{C}(\mu_{\chi})$
$\displaystyle=-\frac{1+2g_{A}^{2}}{2}+\int_{0}^{\mu_{\chi}}\ d|{\mathbf{k}}|\
a_{\chi}(|{\mathbf{k}}|)\,,$ $\displaystyle\bar{\cal A}^{<.{\rm sing}}_{C}$
$\displaystyle=\int_{0}^{\Lambda}\ d|{\mathbf{k}}|\ a_{<}(|{\mathbf{k}}|)\,,$
$\displaystyle\bar{\cal A}^{>}_{C}$ $\displaystyle=\int_{\Lambda}^{\infty}\
d|{\mathbf{k}}|\ a_{>}(|{\mathbf{k}}|)\,,$ (5.5)

with999For the matching condition only the real part of $a_{\chi}$ matters.

$\displaystyle a_{\chi}(|{\mathbf{k}}|)$ $\displaystyle=-(1+2g_{A}^{2})\
\frac{1}{|{\bf k}|}\,\theta(|{\bf k}|-2|{\bf p}|)\,,$ $\displaystyle
a_{<}(|{\mathbf{k}}|)$
$\displaystyle=-\left(g_{V}^{2}({\mathbf{k}}^{2})+2g_{A}({\mathbf{k}}^{2})+\frac{{\mathbf{k}}^{2}g_{M}^{2}({\mathbf{k}}^{2})}{2m_{N}^{2}}\right)\,\frac{r(|{\mathbf{k}}|)}{|{\mathbf{k}}|}\theta(|{\bf
k}|-2|{\bf p}|)\,,$ $\displaystyle a_{>}(|{\mathbf{k}}|)$
$\displaystyle=\frac{3\alpha_{s}}{\pi}\ \bar{g}_{1}^{NN}\
\frac{F_{\pi}^{2}}{|{\mathbf{k}}|^{3}}\,,$ (5.6)

with $r(|{\mathbf{k}}|)$ given in Eq. (4.19). Putting all pieces together we
can write

$2\,\tilde{\cal
C}_{1}(\mu_{\chi})=\frac{1+2g_{A}^{2}}{2}-\int_{0}^{\mu_{\chi}}\
d|{\mathbf{k}}|\ a_{\chi}(|{\mathbf{k}}|)+\int_{0}^{\Lambda}\ d|{\mathbf{k}}|\
a_{<}(|{\mathbf{k}}|)+\int_{\Lambda}^{\infty}\ d|{\mathbf{k}}|\
a_{>}(|{\mathbf{k}}|)\,.$ (5.7)

By construction the $\mu_{\chi}$ dependence is consistent with the
renormalization group equation (RGE) for the rescaled coupling $\tilde{\cal
C}_{1}(\mu_{\chi})$ [43]. This representation shows that up to an additive
constant the LEC $\tilde{\cal C}_{1}(\mu_{\chi})$ can be thought of as the
difference between two integrals in $|{\mathbf{k}}|$, one in the full theory
extending all the way to $|{\mathbf{k}}|\to\infty$, and one in the EFT
extending to $|{\mathbf{k}}|=\mu_{\chi}$. Therefore, the LEC $\tilde{\cal
C}_{1}(\mu_{\chi})$ corresponds to (i) a possible mismatch between the LO
chiral EFT and the full amplitude at $|{\mathbf{k}}|<\mu_{\chi}$; (ii) the
component of the full amplitude arising from $|{\mathbf{k}}|>\mu_{\chi}$.

Figure 6:  Left panel: $a_{\chi}(|{\mathbf{k}}|)$ (yellow),
$a_{<}(|{\mathbf{k}}|)$ (blue), and $a_{>}(|{\mathbf{k}}|)$ (green band)
assuming $\bar{g}_{1}^{NN}\in[-10,10]$. Right panel: same plot for the
absolute values in logarithmic scale.

In Fig. 6 we show $a_{\chi}(|{\mathbf{k}}|)$ (yellow), $a_{<}(|{\mathbf{k}}|)$
(blue), and $a_{>}(|{\mathbf{k}}|)$ (green) assuming
$\bar{g}_{1}^{NN}\in[-10,+10]$. The behavior of the integrand indicates that
the LEC is dominated by the low- and intermediate-momentum regions. In the
left panel of Fig. 7 we show the dependence of $\tilde{\cal
C}_{1}(\mu_{\chi}=M_{\pi})$ on the separation scale $\Lambda$, which proves to
be relatively mild. The right panel shows the dependence of $\tilde{\cal
C}_{1}$ on the chiral renormalization scale $\mu_{\chi}$. The impact of
varying $\bar{g}_{1}^{NN}\in[-10,+10]$ is illustrated by the two different
curves in the left panel of Fig. 7 and is very small.

Figure 7:  Left panel: dependence of $\tilde{\cal C}_{1}(\mu_{\chi}=M_{\pi})$
on the matching scale $\Lambda$ using $\bar{g}_{1}^{NN}=\pm 10$. Right panel:
dependence of $\tilde{\cal C}_{1}$ on the chiral renormalization scale
$\mu_{\chi}$.

Overall, the result for the LEC $\tilde{\cal
C}_{1}(\mu_{\chi}=M_{\pi})\simeq+1.32$ is relatively stable and confirms the
expected scaling of the contact term, i.e., $\tilde{\cal
C}_{1}\sim\mathcal{O}(1)$ [40]. This matching analysis implies that the
contribution from $\tilde{\cal C}_{1}$ and the long-range neutrino exchange
encoded in ${\cal A}_{C}$ interfere destructively, at least when one uses
dimensional regularization with minimal subtraction. Another way to see this
is as follows: the ${}^{1}S_{0}$ long-range neutrino potential $V_{\nu{\rm
L}}=(1+2g_{A}^{2})/{\mathbf{q}}^{2}>0$ has the opposite sign compared to the
short-range potential $V_{\nu{\rm S}}=-2{\cal C}_{1}<0$, since ${\cal
C}_{1}(\mu_{\chi})>0$ for values of $\mu_{\chi}$ appropriate for the chiral
EFT analysis (see the right panel of Fig. 7 for a quantitative statement).

The origin of the sign can be traced back to the mismatch between the LO
chiral EFT and the full amplitude for $|{\mathbf{k}}|<\mu_{\chi}$, an effect
quite visible in the doubly-logarithmic plot in the right panel of Fig. 6. In
turn, the sizable deviation of $a_{<}(|{\mathbf{k}}|)$ (blue) from
$a_{\chi}(|{\mathbf{k}}|)$ (yellow) at $|{\mathbf{k}}|<\mu_{\chi}=M_{\pi}$ is
due to the behavior of $r(|{\mathbf{k}}|)$ (see Fig. 4 and discussion in
Section 4.1), which enters multiplicatively in $a_{<}(|{\mathbf{k}}|)$, as per
Eq. (5.6), and effectively encodes higher-order corrections to the forward
Compton amplitude. As discussed in Section 4.1, the all-important negative
slope of $r(|{\mathbf{k}}|)$ for $|{\mathbf{k}}|<300\,\text{MeV}$ is
controlled by the large ${}^{1}S_{0}$ effective range $r_{0}\simeq
2.7\,\text{fm}$ and is not a model artifact.

While these arguments allow us to understand the origin of the sign of the
contact term at a given scale in the $\overline{\text{MS}}$ scheme, we stress
that no general statement is possible, see Section 7 for more details. For
instance, in the cutoff schemes discussed there, both destructive and
constructive effects are possible depending on the scale of the regulator.

###  5.2 Discussion of uncertainties

Our baseline result $\tilde{\cal C}_{1}(\mu_{\chi}=M_{\pi})\simeq+1.32$ is
obtained using the Kaplan–Steele [87] $V_{S}$ potential. This result is
relatively stable with respect to input parameters in the form factors, the
local matrix element $\bar{g}_{1}^{NN}$ controlling the high-$|{\mathbf{k}}|$
tail, and the matching scale $\Lambda$, as long as $\Lambda\geq
1\,\text{GeV}$. The coupling $\bar{g}_{1}^{NN}(\mu)$, expected to be
$\mathcal{O}(1)$, is presently unknown, but in view of the large value of the
similar two-pion matrix element $\bar{g}_{LR}^{\pi\pi}=8.23$ (see Appendix B
for details), we take the range $\bar{g}_{1}^{NN}\in[-10,+10]$, with minor
impact on the final result. These parametric uncertainties do not exceed
$\delta\tilde{\cal C}_{1}\simeq\pm 0.05_{\rm par}$ and are dominated by the
input $\Lambda_{A}=1.0(2)\,\text{GeV}$.

The main systematics of our approach are due to the following effects:

  1. 1.

We have effectively truncated the spectral representation (2.22) to keep only
the elastic $N\\!N$ intermediate state. In order to assess the size of the
neglected terms, we provide below a rough estimate of the simplest inelastic
channel, namely $N\\!N\pi$,101010In the discussion of higher-order EFT
contributions to the $nn\to pp$ amplitude given in Ref. [38], these
contributions were called “non-factorizable.”  see Fig. 8. As discussed in
Appendix C, we find contributions to $\tilde{\cal C}_{1}$ from the $N\\!N\pi$
intermediate states on the order of $|\delta\tilde{\cal
C}_{1}|=0.1\text{--}0.35$ and due to this we assign an uncertainty of
$\delta\tilde{\cal C}_{1}|_{\rm inelastic}=\pm 0.5$ .

  2. 2.

Within the $N\\!N$ channel, we have built our integrand in $|{\mathbf{k}}|$
starting from the low-$|{\mathbf{k}}|$ end. To extend the form of the
integrand to $|{\mathbf{k}}|\sim\Lambda\sim 1\,\text{GeV}$, we have taken two
main steps: (i) We have included nucleon form factors of the weak currents.
This is known to saturate the elastic contribution to the electromagnetic mass
difference of both the nucleon and pion in the Cottingham approach. (ii) We
have included HOS form factors $f_{S}({\mathbf{q}},{\mathbf{p}})$ that encode
the higher-momentum behavior of the $N\\!N$ scattering amplitude. This effect,
parameterized by the ratio $r(|{\mathbf{k}}|)$ defined in Eq. (4.19), has no
analog in the Cottingham literature, which usually deals with two electroweak
current insertions on a single hadron. As discussed in Section 4.1, we have
performed the modeling of $r(|{\mathbf{k}}|)$ using a three-Yukawa potential
and tested against the NLO chiral EFT analysis. Comparison with the chiral EFT
analysis shows that inclusion of this effect is essential to reproduce the
correct Compton amplitude, and hence the integrand, already at
$|{\mathbf{k}}|\leq M_{\pi}$. On the other hand, the extension of
$r(|{\mathbf{k}}|)$ to $|{\mathbf{k}}|\sim\Lambda\sim 1\,\text{GeV}$
introduces model dependence. To quantify this model dependence we have
employed both the simple three-resonance model and the AV18 model for the
short-range $N\\!N$ interactions in the ${}^{1}S_{0}$ channel. In both cases
the behavior of $r(|{\mathbf{k}}|)$ shows the same features, see Fig. 17: a
rapid drop controlled by the large effective range, a zero for
$|{\mathbf{k}}|$ between 400 and 500 MeV, and then a minimum around
$(600\text{--}700)\,\text{MeV}$. However, the key point is that because of the
fall-off of the integrand multiplying $r(|{\mathbf{k}}|)$, any feature above
500 MeV is washed out, and the differences between AV18 and Kaplan–Steele
potentials have little impact on the extracted LEC. Numerically, compared to
our baseline Kaplan–Steele model, we find a variation on $\tilde{\cal C}_{1}$
of $\approx+0.20$ and $\approx-0.20$ when using the Reid or AV18 potential,
respectively. This difference correlates with the different slopes of
$r(|{\mathbf{k}})$ for $|{\mathbf{k}}|<300\,\text{MeV}$, see Fig. 17. Based on
this, we will assign an error of $\delta\tilde{\cal C}_{1}$ of $\pm 0.2$ due
to the choice of the short-range potential $V_{S}$.

  3. 3.

The above discussion and Fig. 6 both point to the fact that the bulk of
$\tilde{\cal C}_{1}$ is controlled by the behavior of $a_{<}(|{\mathbf{k}}|)$
for $|{\mathbf{k}}|<300\,\text{MeV}$, where our integrand can be linked to the
model-independent chiral EFT behavior. Quantitatively, we find that the region
$|{\mathbf{k}}|\in[0.4,1.5]\,\text{GeV}$, entailing the largest unknowns,
contributes $\Delta|\tilde{\cal C}_{1}|\leq 0.05$, which is a reassuring
result.

  4. 4.

Finally, we note that our approach can be used to estimate the combination of
LECs $\tilde{\cal C}_{1}+\tilde{\cal C}_{2}$, which corresponds to two
insertions of the electromagnetic current and can be extracted from data. As
discussed in Section 6, our matching calculation for $\tilde{\cal
C}_{1}+\tilde{\cal C}_{2}$ compares within uncertainties fairly well with the
value extracted from the CIB combination of $N\\!N$ scattering lengths. This
is a non-trivial validation of the method and gives us confidence that our
results provide a realistic estimate of these couplings and their
uncertainties.

Figure 8:  Representative inelastic contribution to the generalized forward
Compton amplitude. Solid lines represent nucleons, dashed lines pions, and
wavy lines the current insertions.

In light of the above discussion, we quote the following estimate for the LEC

$\tilde{\cal C}_{1}(\mu_{\chi}=M_{\pi})\simeq 1.32(50)_{\rm
inel}(20)_{V_{S}}(5)_{\rm par}=1.3(6)\,,$ (5.8)

which becomes

$\tilde{\cal C}_{1}(\mu_{\chi}=4M_{\pi})=4.2(6)$ (5.9)

at a renormalization scale that corresponds more closely to cutoffs used in
ab-initio many-body calculations (albeit in a different scheme and thus not in
direct correspondence, see Section 7 for more details).

The final uncertainty is dominated by the missing inelastic contributions and
implies a relative precision of (20–30)% on the renormalized singular
amplitude $\bar{\mathcal{A}}_{C}^{\text{sing}}+2\,\tilde{\cal C}_{1}$ at
$|\mathbf{p}|\sim(20\text{--}30)\,\text{MeV}$—in line with the expectation
from the Cottingham analyses of pion and nucleon masses. Note that this
translates into a smaller relative error on the total amplitude
$\mathcal{A}_{\nu}^{\chi\text{EFT}}$, as we will discuss in Section 7.

##  6 Vector–vector amplitude and $\boldsymbol{{\cal C}_{1}+{\cal C}_{2}}$

In analogy to the matching for the purely left-handed coupling ${\cal
C}_{1}=g_{\nu}^{NN}$, one can determine the left–right coupling ${\cal C}_{2}$
and the vector combination ${\cal C}_{1}+{\cal C}_{2}$. Since ${\cal
C}_{1}+{\cal C}_{2}$ can be extracted from a fit to the $I=2$ $N\\!N$
scattering amplitudes, this offers a way to validate the matching approach
through experimental data. The matching calculation is most easily done by
considering the $nn\to pp$ transition in the unphysical theory in which the
$W$ boson couples to both the left-handed quark current
$\bar{u}_{L}\gamma_{\mu}d_{L}$ (as in the Standard Model) and the right-handed
current $\bar{u}_{R}\gamma_{\mu}d_{R}$ with equal couplings.

###  6.1 Effective Lagrangians

Let us briefly recall the basic elements of this extended analysis [40, 43].
The starting point is the quark-level electromagnetic and weak Lagrangian in
the Standard Model

$\mathcal{L}=\bar{q}_{L}\gamma^{\mu}\left(l_{\mu}+\hat{l}_{\mu}\right)q_{L}+\bar{q}_{R}\gamma^{\mu}\left(r_{\mu}+\hat{r}_{\mu}\right)q_{R}\,,$
(6.1)

where $q$ denotes the quark doublet $q=(u,\;d)^{T}$ and we defined

$\displaystyle l_{\mu}$
$\displaystyle=\frac{e}{2}A_{\mu}\,\tau^{3}-2\sqrt{2}\,G_{F}\left[V_{ud}\,\bar{e}_{L}\gamma_{\mu}\nu_{L}\,\tau^{+}+{\rm
h.c.}\right]\,,$ $\displaystyle\hat{l}_{\mu}$
$\displaystyle=\frac{e}{6}A_{\mu}\,,$ $\displaystyle r_{\mu}$
$\displaystyle=\frac{e}{2}A_{\mu}\,\tau^{3}\,,$ $\displaystyle\hat{r}_{\mu}$
$\displaystyle=\frac{e}{6}A_{\mu}\,,$ (6.2)

where $\tau^{+}=(\tau_{1}+i\tau_{2})/2$. We are interested in the isovector
components of the interactions, involving $l_{\mu}$ and $r_{\mu}$. We will
consider the unphysical theory in which the $W$ boson has vector-like
couplings and hence couples with same strength to both
$\bar{u}_{L}\gamma_{\mu}d_{L}$ and $\bar{u}_{R}\gamma_{\mu}d_{R}$. This
amounts to $r_{\mu}=l_{\mu}$, namely

$r_{\mu}=\frac{e}{2}A_{\mu}\,\tau^{3}-2\sqrt{2}\,G_{F}\left[V_{ud}\,\bar{e}_{L}\gamma_{\mu}\nu_{L}\,\tau^{+}+{\rm
h.c.}\right]\,.$ (6.3)

Double insertions of the isovector component of Eq. (6.1) give rise to
$0\nu\beta\beta$ amplitudes and $I=2$ electromagnetic effects through the
effective actions:

$S_{\rm LNV}=S^{++}\,,\qquad S_{\rm EM,I=2}=S^{33}\,,$ (6.4)

with

$S^{aa}=F^{aa}\,\int d^{4}xd^{4}y\ g^{\mu\nu}S(x-y)\ T\Big{(}J^{a}_{\mu
L}(x)\,J^{a}_{\nu L}(y)+J^{a}_{\mu R}(x)\,J^{a}_{\nu R}(y)+2J^{a}_{\mu
L}(x)\,J^{a}_{\nu R}(y)\Big{)}\,,$ (6.5)

and

$J^{a}_{\mu L}(x)=\bar{q}_{L}(x)\gamma_{\mu}t^{a}q_{L}(x)\,,\qquad J^{a}_{\mu
R}(x)=\bar{q}_{R}(x)\gamma_{\mu}t^{a}q_{R}(x)\,.$ (6.6)

The interaction-specific factors are111111For simplicity we neglect the
electron four-momenta, so we can replace $e_{L}(x),e_{L}(y)\to e_{L}(0)$ and
factor the electron fields out of the convolution integral.

$\displaystyle F^{++}$
$\displaystyle=\frac{8G_{F}^{2}V_{ud}^{2}m_{\beta\beta}}{2!}\
\bar{e}_{L}e_{L}^{c}\,,$ $\displaystyle t^{+}$
$\displaystyle=\tau^{+}=\frac{\tau_{1}+i\tau_{2}}{2}\,,$ $\displaystyle
F^{33}$ $\displaystyle=\frac{e^{2}}{2!}\,,$ $\displaystyle t^{3}$
$\displaystyle=\frac{\tau^{3}}{2}\,,$ (6.7)

and $S(x-y)$ is the massless scalar propagator defined in Eq. (2.4), which
arises from the neutrino or photon propagator in the $0\nu\beta\beta$ and
$I=2$ cases, respectively.

At low energy, the above effective actions manifest themselves through both
(i) long-distance effects with exchange of soft, potential, and ultrasoft
neutrinos (photons) between the pion and nucleon realization of the
electroweak currents; (ii) local interactions, which can be thought of as
arising from the exchange of hard neutrino (photon) modes. Since the
amplitudes transform according to the same irreducible representation of the
isospin group, in the isospin symmetry limit the exchange of hard neutrinos
leads to identical contributions as photon exchange, up to the overall factors
$F^{++}$ and $F^{33}$.

The above considerations imply the following relations between low-energy
chiral Lagrangians. In the nucleon sector we have:

$\displaystyle\mathcal{L}_{e^{2}}^{NN}$ $\displaystyle=\frac{e^{2}}{4}\
\Big{[}{\cal C}_{1}\,O^{33}_{1}+{\cal C}_{2}\,O^{33}_{2}\Big{]}\,,$
$\displaystyle\mathcal{L}_{|\Delta L|=2}^{NN}$
$\displaystyle=2\,G_{F}^{2}V_{ud}^{2}m_{\beta\beta}\ \bar{e}_{L}e_{L}^{c}\
\Big{[}{\cal C}_{1}\,O^{++}_{1}+{\cal C}_{2}\,O^{++}_{2}\Big{]}\,,$ (6.8)

with

$\displaystyle O_{1}^{aa}$
$\displaystyle=\bar{N}\mathcal{Q}^{a}_{L}N\,\bar{N}\mathcal{Q}^{a}_{L}N-\frac{\textrm{Tr}[\mathcal{Q}^{a}_{L}\mathcal{Q}^{a}_{L}]}{6}\bar{N}\mbox{\boldmath$\tau$}N\cdot\bar{N}\mbox{\boldmath$\tau$}N+\\{L\leftrightarrow
R\\}\,,$ $\displaystyle O_{2}^{aa}$
$\displaystyle=2\left(\bar{N}\mathcal{Q}^{a}_{L}N\bar{N}\mathcal{Q}^{a}_{R}N-\frac{\textrm{Tr}[\mathcal{Q}^{a}_{L}\mathcal{Q}^{a}_{R}]}{6}\bar{N}\mbox{\boldmath$\tau$}N\cdot\bar{N}\mbox{\boldmath$\tau$}N\right)\,,$
(6.9)

where

$\mathcal{Q}_{L}^{a}=u^{\dagger}t^{a}u\,,\qquad\mathcal{Q}_{R}^{a}=ut^{a}u^{\dagger}\,,$
(6.10)

and $u^{2}=U=\exp(i\mbox{\boldmath$\tau$}\cdot\boldsymbol{\pi}/F_{\pi})$
incorporates the pion fields.

In the pion sector we have non-derivative LO local operators arising from the
LR current correlator:

$\displaystyle\mathcal{L}_{e^{2}}^{\pi\pi}$ $\displaystyle=e^{2}ZF_{\pi}^{4}\
\textrm{Tr}[\mathcal{Q}^{3}_{L}\mathcal{Q}^{3}_{R}]\,,$
$\displaystyle\mathcal{L}_{|\Delta L|=2}^{\pi\pi}$
$\displaystyle=\left(8G_{F}^{2}V_{ud}^{2}m_{\beta\beta}\
\bar{e}_{L}e_{L}^{c}\right)\ ZF_{\pi}^{4}\
\textrm{Tr}[\mathcal{Q}^{+}_{L}\mathcal{Q}^{+}_{R}]\,,$ (6.11)

where, at LO in chiral perturbation theory (ChPT), $Z$ is related to the pion-
mass splitting by121212Note that the operator in the first line of Eq. (6.11)
differs from the usual definition in ChPT [93] by an inessential constant.

$Ze^{2}F_{\pi}^{2}=\frac{1}{2}\delta
M^{2}_{\pi}=\frac{1}{2}\left(M_{\pi^{\pm}}^{2}-M_{\pi^{0}}^{2}\right).$ (6.12)

###  6.2 Pion two-point function and the low-energy constant $\boldsymbol{Z}$

Before discussing the $nn\to pp$ amplitude, it is instructive to consider the
simpler $\pi^{-}\to\pi^{+}$ transition. In complete analogy with Eq. (2.6) one
can define a vector–vector correlator via the replacement
$J_{\mu}^{L}=\bar{u}_{L}\gamma_{\mu}d_{L}\to
V_{\mu}=J_{\mu}^{L}+J_{\mu}^{R}=\bar{u}\gamma_{\mu}d$. The object of interest
for our matching calculation is

$T_{\mu\nu}^{\pi\pi}(k,p)=\langle\pi^{+}(p)|\ \hat{\Pi}^{VV}_{\mu\nu}(k,0)\
|\pi^{-}(p)\rangle\,,\qquad{\cal T}^{\pi\pi}(k,p)\equiv
g^{\mu\nu}T_{\mu\nu}^{\pi\pi}(k,p)\,.$ (6.13)

${\cal T}^{\pi\pi}(k,0)$ determines the LEC $Z$ defined in Eq. (6.11), while
${\cal T}^{\pi\pi}(k,p)$ (where $p$ is an off-shell four-momentum with
potential scaling) contributes to $nn\to pp$ as discussed below.

In the small $p,k$ regime ${\cal T}^{\pi\pi}(k,p)$ can be computed reliably in
ChPT. At LO one has131313Up to a proportionality factor due to isospin, this
expression agrees with the forward Compton amplitude, see for example Ref.
[66].

${\cal
T}_{\chi}^{\pi\pi}(k,p)=2i\left[2g_{\mu\nu}-\frac{(2p+k)_{\mu}(2p+k)_{\nu}}{(p+k)^{2}-M_{\pi}^{2}+i\epsilon}-\frac{(2p-k)_{\mu}(2p-k)_{\nu}}{(p-k)^{2}-M_{\pi}^{2}+i\epsilon}\right]\
g^{\mu\nu}\,,$ (6.14)

which in the chiral limit and for $p\to 0$ reduces to

${\cal
T}_{\chi}^{\pi\pi}(k,0)\Big{|}_{M_{\pi}=0}=4i\left[g_{\mu\nu}-\frac{k_{\mu}k_{\nu}}{k^{2}}\right]\,g^{\mu\nu}=4i(d-1)\,.$
(6.15)

The expression for ${\cal T}^{\pi\pi}(k,p)$ can be extended to the
intermediate $k$ momentum region by using resonance models that work
reasonably well in the meson sector [62] and automatically respect the
constraints of chiral symmetry beyond the LO amplitude (6.14). We will not
follow this route but we will rather use the following prescription

${\cal T}_{<}^{\pi\pi}(k,p)={\cal
T}_{\chi}^{\pi\pi}(k,p)\times\left(F_{\pi}^{V}(k^{2})\right)^{2}\,,$ (6.16)

which captures the full pion-pole contribution to the on-shell Compton
amplitude [94, 95, 96] in terms of the pion vector form factor
$F_{\pi}^{V}(k^{2})$, for which we will take the simple monopole form

$F_{\pi}^{V}({k}^{2})=\frac{M_{V}^{2}}{M_{V}^{2}-k^{2}}\,,$ (6.17)

with $M_{V}=M_{\rho}$. While the form (6.16) does not capture the
contributions form intermediate-state axial-vector mesons [63, 62, 64, 65,
66], the resulting numerical impact on the LEC $Z$ is below the 20% level
[64].

In the large $k$ regime, the form of ${\cal T}^{\pi\pi}(k,p)$ is dictated by
the OPE, see Section 4.2 and Appendix B. It is given by

${\cal
T}_{>}^{\pi\pi}(k,p)=\frac{4ig_{s}^{2}}{\left(k^{2}+i\epsilon\right)^{2}}\,\langle\pi^{+}(p)|\
\hat{O}^{VV}\
|\pi^{-}(p)\rangle\,,\qquad\hat{O}^{VV}=O_{1}+\frac{1}{2}\left(O_{4}-3O_{5}\right)\,.$
(6.18)

To LO in ChPT, the pion matrix elements are given by

$\langle\pi^{+}(p)|\ \hat{O}^{VV}\ |\pi^{-}(p)\rangle=F_{0}^{2}\
\left[\frac{5}{3}g_{1}^{\pi\pi}p^{2}+\frac{1}{2}\left(g_{4}^{\pi\pi}-3g_{5}^{\pi\pi}\right)\right]\,,$
(6.19)

where the LECs $g^{\pi\pi}_{1,4,5}$ are known from lattice QCD [92, 90] and
scale as $g_{1}^{\pi\pi}\sim\mathcal{O}(1)$ and
$g_{4,5}^{\pi\pi}\sim\mathcal{O}(\Lambda_{\chi}^{2})$. From this scaling it is
clear that the dominant contribution to the two-pion matrix element is
proportional to the combination
$g_{LR}^{\pi\pi}\equiv(1/2)(g_{4}^{\pi\pi}-3g_{5}^{\pi\pi})$.

We are now in a position to write down a matching relation for the LEC $Z$
introduced in Eq. (6.11). Writing the pion two-point function at zero momentum
in ChPT and in full QCD, one derives the relation

$8ZF_{0}^{2}=\int\frac{d^{4}k}{(2\pi)^{4}}\ \frac{1}{k^{2}+i\epsilon}\,{\cal
T}^{\pi\pi}(k,0)\,.$ (6.20)

In the spirit of the matching strategy used in this work, we split up the
integration over three-momentum $|{\mathbf{k}}|$ in two regions and by using
${\cal T}_{<}^{\pi\pi}(k,0)$ from Eq. (6.15) and ${\cal T}_{>}^{\pi\pi}(k,0)$
from Eqs. (6.18)–(6.19) we arrive at (see Appendix D for details)

$\displaystyle Z$ $\displaystyle=Z^{<}+Z^{>}\,,$ $\displaystyle Z^{<}$
$\displaystyle=\frac{6}{(4\pi
F_{0})^{2}}\,\int_{0}^{\Lambda}\,d|{\mathbf{k}}|\,|{\mathbf{k}}|\times\left(1-\frac{|{\mathbf{k}}|}{\omega_{V}}\left(1+\frac{M_{V}^{2}}{2\omega_{V}^{2}}\right)\right)\,,$
$\displaystyle Z^{>}$
$\displaystyle=\frac{3\alpha_{s}(\mu)g_{LR}^{\pi\pi}(\mu)}{16\pi}\,\int_{\Lambda}^{\infty}d|{\mathbf{k}}|\,\frac{1}{|{\mathbf{k}}|^{3}}\,,$
(6.21)

where $\omega_{V}=\sqrt{{\mathbf{k}}^{2}+M_{V}^{2}}$ and $\mu$ is the QCD
renormalization scale. The above result is in agreement with estimates of the
pion electromagnetic mass splitting [63, 62, 64, 65, 66]. Taking
$M_{V}=775\,\text{MeV}$, $F_{\pi}=92.28\,\text{MeV}$, and
${g}_{LR}^{\pi\pi}=8.23$ from lattice QCD [90], the sum of low- and high-
momentum components $Z^{<}(\Lambda)+Z^{>}(\Lambda)$ equals 0.63 at
$\Lambda=2\,\text{GeV}$ and reaches the asymptotic value
$Z^{<}(\Lambda\to\infty)=0.67$. The deficit compared to the experimental value
$Z\simeq 0.8$ is understood in terms of the neglected inelastic contributions
from the axial-vector resonances [63, 62, 64], thus providing another estimate
of the error due to neglecting inelastic corrections.

###  6.3 $\boldsymbol{nn\to pp}$ vector-like amplitude in chiral EFT and full
theory

Let us now consider the amplitude ${\cal A}_{VV}$ in the unphysical theory in
which the $W$ boson has vector-like couplings to quarks. Denoting by ${\cal
A}_{LL,LR}$ the $nn\to pp$ amplitudes generated by $W$-exchange between two
left-handed quark currents and a left-handed and a right-handed quark current,
respectively (see Eqs. (2.6) and (2.9)), we have

${\cal A}_{VV}=2{\cal A}_{LL}+2{\cal A}_{LR}\,.$ (6.22)

Note that ${\cal A}_{LL}\equiv{\cal A}_{\nu}$, i.e., the amplitude in the
physical theory, defined in Eq. (2.9), and the factor of 2 in front of ${\cal
A}_{LL}$ arises from the fact that $LL$ and $RR$ products, both present due to
vector-like couplings, give the same result by parity.

####  6.3.1 $\boldsymbol{{\cal A}_{VV}}$ in chiral EFT

At LO in chiral EFT, the $nn\to pp$ vector amplitude takes the form

${\cal A}_{VV}^{\chi{\rm EFT}}=2{\cal A}_{LL}^{NN}+2{\cal A}_{LR}^{NN}+{\cal
A}_{VV}^{\pi\pi}\,,$ (6.23)

where ${\cal A}_{LL}^{NN}$ is given in Eq. (3.15). ${\cal A}^{NN}_{LR}$ can be
obtained from ${\cal A}^{NN}_{LL}$ by flipping the sign of the “$A\times A$”
axial contribution, which in practice amounts to setting $g_{A}^{2}\to-
g_{A}^{2}$ everywhere and replacing ${\cal C}_{1}\to{\cal C}_{2}$ in the
counter term amplitude.

In the vector theory, the main new effect compared to the physical amplitude
is the presence of a LO contribution induced by pion exchange, denoted by
${\cal A}_{VV}^{\pi\pi}$ in Eq. (6.23), induced by the non-derivative operator
in Eq. (6.11) with coupling constant $Z$. The amplitude ${\cal
A}_{VV}^{\pi\pi}$ has the same form of ${\cal A}_{\nu}^{\chi{\rm EFT}}$ in Eq.
(3.4), except for the fact that in each diagram the neutrino propagator is
replaced by a pion propagator with one insertion of $Z$ that converts a
$\pi^{-}$ into a $\pi^{+}$, see Fig. 9. This implies that in the computation
of ${\cal A}_{A,B,C}$ in Eqs. (3.8)–(3.15) one needs to make the replacement:

$V^{{}^{1}S_{0}}_{\nu\,{\rm
L}}({\mathbf{k}})=\frac{1}{{\mathbf{k}}^{2}}\left[1+2g_{A}^{2}+\frac{g_{A}^{2}M_{\pi}^{4}}{({\mathbf{k}}^{2}+M_{\pi}^{2})^{2}}\right]\
\to\ V^{{}^{1}S_{0}}_{Z}({\mathbf{k}})\ =\frac{8Zg_{A}^{2}}{{\mathbf{k}}^{2}}\
\left(1-\frac{M_{\pi}^{2}}{{\mathbf{k}}^{2}+M_{\pi}^{2}}\right)^{2}\,.$ (6.24)

As discussed in Ref. [43], the UV behavior of this potential induces
additional divergences reabsorbed by the coupling ${\cal C}_{2}$.

Finally, we note that in complete analogy to the physical $LL$ case, the
matching condition for the $VV$ amplitude will involve only the real part of
the singular component of the amplitude ${\cal A}_{C}$ associated with Eq.
(6.23). Combining the $NN$ and $\pi\pi$ contributions, we write it in the
following way that will prove useful in the matching procedure:

$\displaystyle{\cal A}_{C}^{\rm sing}(\mu_{\chi})\Big{|}_{VV}$
$\displaystyle=-2\frac{m_{N}^{2}}{(4\pi)^{2}}\,\left(1+2Zg_{A}^{2}\right)\,\left[\log\frac{\mu_{\chi}^{2}}{4|{\bf
p}|^{2}}+1\right]$
$\displaystyle=2\frac{m_{N}^{2}}{(4\pi)^{2}}\,\left[-(1+2Zg_{A}^{2})-4Zg_{A}^{2}\,\log\frac{\mu_{\chi}}{2|{\bf
p}|}+\int_{0}^{\mu_{\chi}}\ d|{\mathbf{k}}|\
a^{VV}_{\chi}(|{\mathbf{k}}|)\right]\,,$ $\displaystyle
a^{VV}_{\chi}(|{\mathbf{k}}|)$ $\displaystyle=-2\ \frac{1}{|{\bf
k}|}\,\theta(|{\bf k}|-2|{\bf p}|)\,.$ (6.25)

Figure 9: Diagrammatic representation of ${\cal A}_{VV}^{\pi\pi}$, the pion
contributions to the LO $nn\rightarrow ppee$ amplitude in the vector-like
unphysical theory. Double, dashed, and plain lines denote nucleons, pions, and
electrons, respectively. The black square denotes an insertion of the non-
derivative coupling $Z$, see Eq. (6.11). The blue ellipse represents iteration
of $\hat{V}_{\pi}$, as in Fig. 2. The ellipses in the second and third lines
denote diagrams with arbitrary numbers of $N\\!N$ bubble insertions.

####  6.3.2 $\boldsymbol{{\cal A}_{VV}}$ in the full theory

As before, we split the full amplitude into two terms, capturing the
contributions from small plus intermediate and hard neutrino momenta: ${\cal
A}_{VV}^{\rm full}={\cal A}_{VV}^{<}+{\cal A}_{VV}^{>}$. One can again
organize the full theory calculation of ${\cal A}_{VV}^{<}$ in close analogy
to the EFT expression given in Eq. (6.23). ${\cal A}_{LL}^{NN,<}$ is discussed
in Section 4.1 and the corresponding ${\cal A}_{LR}^{NN,<}$ is obtained by
defining $\hat{O}^{LR}({\mathbf{k}})$ in complete analogy to
$\hat{O}^{LL}({\mathbf{k}})$ in Eq. (2.18) and then following through. In
practice, this again amounts to flipping the sign of the axial contribution in
${\cal A}_{LL}^{NN,<}$. As discussed in Section 5, for matching purposes one
only needs the singular parts of the amplitude ${\cal A}_{C}^{<}$, which can
be written as

${\cal A}_{C}^{<,{\rm sing}}\Big{|}_{VV}=2\left({\cal A}_{C}^{<,{\rm
sing}}\Big{|}_{LL}+{\cal A}_{C}^{<,{\rm sing}}\Big{|}_{LR}\right)+{\cal
A}_{C}^{\pi\pi,<,{\rm sing}}\,.$ (6.26)

${\cal A}_{C}^{<,{\rm sing}}|_{LL}$ is given in Eqs. (4.18), (5.5), and (5.6).
${\cal A}_{C}^{<,{\rm sing}}|_{LR}$ is obtained by flipping the sign of the
axial terms in ${\cal A}_{C}^{<,{\rm sing}}|_{LL}$, and the sum of the $N\\!N$
intermediate-state contributions is

$\displaystyle 2\left({\cal A}_{C}^{<,{\rm sing}}\Big{|}_{LL}+{\cal
A}_{C}^{<,{\rm sing}}\Big{|}_{LR}\right)$
$\displaystyle=2\,\frac{m_{N}^{2}}{(4\pi)^{2}}\,\int_{0}^{\Lambda}\
d|{\mathbf{k}}|\ a^{NN}_{<}(|{\mathbf{k}}|)\,,$ $\displaystyle
a^{NN}_{<}(|{\mathbf{k}}|)$
$\displaystyle=-\left(2g_{V}^{2}({\mathbf{k}}^{2})+\frac{{\mathbf{k}}^{2}g_{M}^{2}({\mathbf{k}}^{2})}{m_{N}^{2}}\right)\,8\,{\cal
I}_{C}^{<}(|{\mathbf{k}}|)\,.$ (6.27)

Figure 10:  Contribution to ${\cal A}_{VV}$ due to the insertion of two vector
currents on a potential pion line. In the small $|{\mathbf{k}}|$ region this
diagram represents ${\cal A}_{C}^{\pi\pi,<,{\rm sing}}$. In the large
$|{\mathbf{k}}|$ region it represents $\hat{\cal A}_{C}^{\rm sing}$. See main
text for notation.

In the full theory analysis, the counterpart of the EFT terms proportional to
$Z$ (called ${\cal A}_{VV}^{\pi\pi}$ in Eq. (6.23)) is due to the insertion of
$\hat{O}^{VV}({\mathbf{k}})\equiv(\hat{O}^{LL}+\hat{O}^{RR}+2\hat{O}^{LR})({\mathbf{k}})$
on a pion exchanged between the two nucleons, as depicted in Fig. 10. The
corresponding ${\cal A}_{C}^{\pi\pi,<,{\rm sing}}$ is given by

${\cal A}_{C}^{\pi\pi,<,{\rm
sing}}=-m_{N}^{2}\,\frac{g_{A}^{2}}{F_{0}^{2}}\int\frac{d^{4}k}{(2\pi)^{4}}\
\frac{1}{k^{2}+i\epsilon}\,\int\frac{d^{3}{\mathbf{q}}}{(2\pi)^{3}}\
\frac{1}{{\mathbf{q}}^{2}}\,{\cal I}_{C}^{<}(|{\mathbf{q}}|)\ {\cal
T}_{<}^{\pi\pi}(k,q)\,.$ (6.28)

The integral in $d^{3}{\mathbf{q}}$ is actually convergent. This can be seen
by rewriting ${\cal T}_{<}^{\pi\pi}(k,p)$ as follows:

$\displaystyle{\cal T}_{<}^{\pi\pi}(k,q)$
$\displaystyle=2i\bigg{[}2(d-4)-F_{+}(k,q)-F_{-}(k,q)\bigg{]}\times\left(F_{\pi}^{V}(k^{2})\right)^{2}\,,$
$\displaystyle F_{\pm}(k,q)$ $\displaystyle=\frac{4M_{\pi}^{2}-3k^{2}\mp
4k\cdot q}{(q\pm k)^{2}-M_{\pi}^{2}+i\epsilon}\,.$ (6.29)

Note that terms involving $q^{0}\sim|{\mathbf{q}}|^{2}/m_{N}$ are neglected at
the order at which we work. Up to a finite piece coming from the term
proportional to $d-4$ in ${\cal T}_{<}^{\pi\pi}$, the integral in Eq. (6.28)
can be evaluated as follows: first performing the integration in $k^{0}$ using
the residues’ theorem; then performing the angular integral in
$d^{3}{\mathbf{q}}$; finally evaluating the integral in $d|{\mathbf{q}}|$. We
will present our results in order of complexity:

  1. 1.

First, we will take $F_{\pi}^{V}(k^{2})=1$ and ${\cal
I}_{C}^{<}(|{\mathbf{q}}|)=\theta(|{\mathbf{q}}|-2|{\mathbf{p}}|)/(8|{\mathbf{q}}|)$.

  2. 2.

Next, we will use the realistic pion electromagnetic form factor
$F_{\pi}^{V}(k^{2})$ from Eq. (6.17).

  3. 3.

Finally, we will include higher-momentum effects in the $N\\!N$ interactions
by writing

${\cal
I}_{C}^{<}(|{\mathbf{q}}|)=\theta(|{\mathbf{q}}|-2|{\mathbf{p}}|)\frac{r(|{\mathbf{q}}|)}{8|{\mathbf{q}}|}\,,$
(6.30)

with $r(|{\mathbf{q}}|)$ defined in Eq. (4.19) and illustrated in Fig. 4.

Throughout, we work in the chiral limit ($M_{\pi}\to 0$).

The result is quite simple when using $F_{\pi}^{V}(k^{2})=1$ (corresponding to
the limit $M_{V}\to\infty$) and ${\cal
I}_{C}^{<}(|{\mathbf{q}}|)=\theta(|{\mathbf{q}}|-2|{\mathbf{p}}|)/(8|{\mathbf{q}}|)$:

${\cal A}_{C}^{\pi\pi,<,{\rm
sing}}=-\frac{m_{N}^{2}}{(4\pi)^{2}}\,\frac{16g_{A}^{2}}{(4\pi F_{0})^{2}}\
\int_{0}^{\Lambda}d|{\mathbf{k}}||{\mathbf{k}}|\left(4\log
2-2+3\log\frac{|{\mathbf{k}}|}{2|{\mathbf{p}}|}\right)\,.$ (6.31)

This expression can be recast into the following useful form using Eq. (6.21)

$\displaystyle{\cal A}_{C}^{\pi\pi,<,{\rm sing}}$
$\displaystyle=2\,\frac{m_{N}^{2}}{(4\pi)^{2}}\,\left[-4g_{A}^{2}Z^{<}\,\log\frac{\nu_{\chi}}{{2|{\mathbf{p}}|}}\
+\ \ \int_{0}^{\Lambda}d|{\mathbf{k}}|\
a_{<}^{\pi\pi}(|{\mathbf{k}}|,\nu_{\chi})\right]\,,$ $\displaystyle
a_{<}^{\pi\pi}(|{\mathbf{k}}|,\nu_{\chi})$ $\displaystyle=-8g_{A}^{2}\
\frac{|{\mathbf{k}}|}{(4\pi F_{0})^{2}}\left(4\log
2-2+3\log\frac{|{\mathbf{k}}|}{\nu_{\chi}}\right)\,,$ (6.32)

where we have extracted a factor of two for later convenience, have introduced
the arbitrary scale $\nu_{\chi}$, and isolated the dependence on the IR
physics (external momentum $|{\mathbf{p}}|$), which must reproduce the same
dependence in the EFT amplitude (cf. the term in Eq. (6.25) proportional to
$Z$).

For a realistic pion form factor (6.17), the residues from the $\rho$
propagators need to be added, which changes the bracket in Eq. (6.31) to

$\displaystyle 4\log
2-\frac{(\omega_{V}-|{\mathbf{k}}|)\big{(}2\omega_{V}^{2}-3|{\mathbf{k}}|\omega_{V}-|{\mathbf{k}}|^{2}\big{)}}{2\omega_{V}^{3}}+\frac{(\omega_{V}-|{\mathbf{k}}|)^{2}(2\omega_{V}+|{\mathbf{k}}|)}{2\omega_{V}^{3}}\bigg{(}3\log\frac{\omega_{V}+|{\mathbf{k}}|}{2|{\mathbf{p}}|}-1\bigg{)}$
$\displaystyle-3\log\bigg{(}1+\frac{\omega_{V}}{|{\mathbf{k}}|}\bigg{)}-\log\bigg{(}1+\frac{|{\mathbf{k}}|}{\omega_{V}}\bigg{)}\,.$
(6.33)

In the limit $\omega_{V}\to\infty$ this reduces to Eq. (6.31), and for
$|{\mathbf{k}}|\to\infty$ this expression displays the $1/|{\mathbf{k}}|^{4}$
fall-off that makes the integrand in Eq. (6.31) consistent with the OPE
behavior at large $|\mathbf{k}|$. To extract the new form of the integrand, we
first express $Z^{<}$, see Eq. (6.21), as

$\displaystyle
Z^{<}=\frac{3}{16\pi^{2}F_{0}^{2}}\int_{0}^{\Lambda}d|{\mathbf{k}}||{\mathbf{k}}|\frac{(\omega_{V}-|{\mathbf{k}}|)^{2}(2\omega_{V}+|{\mathbf{k}}|)}{\omega_{V}^{3}}\,.$
(6.34)

In this way, Eq. (6.32) still applies upon the replacement

$\displaystyle a_{<}^{\pi\pi}(|{\mathbf{k}}|,\nu_{\chi})$
$\displaystyle=-4g_{A}^{2}\ \frac{|{\mathbf{k}}|}{(4\pi
F_{0})^{2}}\Bigg{\\{}8\log
2-\frac{(\omega_{V}-|{\mathbf{k}}|)\big{(}2\omega_{V}^{2}-3|{\mathbf{k}}|\omega_{V}-|{\mathbf{k}}|^{2}\big{)}}{\omega_{V}^{3}}$
$\displaystyle+\frac{(\omega_{V}-|{\mathbf{k}}|)^{2}(2\omega_{V}+|{\mathbf{k}}|)}{\omega_{V}^{3}}\bigg{(}3\log\frac{\omega_{V}+|{\mathbf{k}}|}{\nu_{\chi}}-1\bigg{)}$
$\displaystyle-6\log\bigg{(}1+\frac{\omega_{V}}{|{\mathbf{k}}|}\bigg{)}-2\log\bigg{(}1+\frac{|{\mathbf{k}}|}{\omega_{V}}\bigg{)}\Bigg{\\}}\,.$
(6.35)

Finally, we include the HOS form factor effects by rewriting in Eq. (6.28)

${\cal
I}_{C}^{<}(|{\mathbf{q}}|)=\theta(|{\mathbf{q}}|-2|{\mathbf{p}}|)\times\left[\frac{1}{8|{\mathbf{q}}|}+\frac{r(|{\mathbf{q}}|)-1}{8|{\mathbf{q}}|}\right]\,.$
(6.36)

The first term gives the analytic result in Eqs. (6.32) and (6.3.2). The
second term gives a correction to $a_{<}^{\pi\pi}(|{\mathbf{k}}|,\nu_{\chi})$,
which we denote by $\delta a_{<}^{\pi\pi}(|{\mathbf{k}}|)$,

$\delta a_{<}^{\pi\pi}(|{\mathbf{k}}|)=4g_{A}^{2}\ \frac{|{\mathbf{k}}|}{(4\pi
F_{0})^{2}}\
\int_{2|{\mathbf{p}}|}^{\infty}d|{\mathbf{q}}|\,\frac{1-r(|{\mathbf{q}}|)}{|{\mathbf{q}}|}\,g(|{\mathbf{k}}|,|{\mathbf{q}}|)\,,$
(6.37)

where $g(|{\mathbf{k}}|,|{\mathbf{q}}|)$ is given by

$\displaystyle g(|{\mathbf{k}}|,|{\mathbf{q}}|)$
$\displaystyle=\begin{cases}g_{-}(|{\mathbf{k}}|,|{\mathbf{q}}|)\qquad\text{if}\qquad|{\mathbf{q}}|<|{\mathbf{k}}|\\\
g_{+}(|{\mathbf{k}}|,|{\mathbf{q}}|)\qquad\text{if}\qquad|{\mathbf{q}}|>|{\mathbf{k}}|\end{cases},\qquad\omega=\omega_{V}=\sqrt{\mathbf{k}^{2}+M_{V}^{2}}\,,$
(6.38) $\displaystyle g_{-}(|{\mathbf{k}}|,|{\mathbf{q}}|)$
$\displaystyle=\frac{\omega-|{\mathbf{k}}|}{\omega^{2}}\bigg{[}\frac{(\omega+|{\mathbf{k}}|)^{2}(|{\mathbf{k}}|+3\omega)}{(\omega+|{\mathbf{k}}|)^{2}-|{\mathbf{q}}|^{2}}-\frac{2(|{\mathbf{k}}|^{2}+2|{\mathbf{k}}|\omega-\omega^{2})}{\omega}\bigg{]}$
$\displaystyle-\frac{2|{\mathbf{q}}|}{|{\mathbf{k}}|}\log\frac{2|{\mathbf{k}}|+|{\mathbf{q}}|}{2|{\mathbf{k}}|-|{\mathbf{q}}|}-\frac{|{\mathbf{k}}|^{4}+6|{\mathbf{q}}|^{2}\omega^{2}+\omega^{4}-2|{\mathbf{k}}|^{2}(|{\mathbf{q}}|^{2}+\omega^{2})}{2|{\mathbf{q}}|\omega^{3}}\log\frac{\omega+|{\mathbf{k}}|-|{\mathbf{q}}|}{\omega+|{\mathbf{k}}|+|{\mathbf{q}}|}\,,$
$\displaystyle g_{+}(|{\mathbf{k}}|,|{\mathbf{q}}|)$
$\displaystyle=\frac{\omega-|{\mathbf{k}}|}{\omega^{2}}\bigg{[}\frac{(|{\mathbf{q}}|-2\omega)(|{\mathbf{k}}|\omega+(|{\mathbf{q}}|+\omega)^{2})}{(\omega+|{\mathbf{q}}|)^{2}-|{\mathbf{k}}|^{2}}-\frac{\omega|{\mathbf{q}}|(|{\mathbf{q}}|-6\omega)+|{\mathbf{k}}|(\omega+|{\mathbf{k}}|)(2|{\mathbf{q}}|-\omega)}{\omega|{\mathbf{q}}|}\bigg{]}$
$\displaystyle-\frac{2|{\mathbf{q}}|}{|{\mathbf{k}}|}\log\frac{2|{\mathbf{k}}|+|{\mathbf{q}}|}{|{\mathbf{q}}|}-\frac{|{\mathbf{k}}|^{4}+6|{\mathbf{q}}|^{2}\omega^{2}+\omega^{4}-2|{\mathbf{k}}|^{2}(|{\mathbf{q}}|^{2}+\omega^{2})}{2|{\mathbf{q}}|\omega^{3}}\log\frac{\omega-|{\mathbf{k}}|+|{\mathbf{q}}|}{\omega+|{\mathbf{k}}|+|{\mathbf{q}}|}\,.$

The integration in Eq. (6.37) will be done numerically.

Putting together the $N\\!N$ and $\pi\pi$ contributions, the full low- and
intermediate-$|{\mathbf{k}}|$ integrand in ${\cal A}_{C}^{<,\rm sing}|_{VV}$
is given by

$a^{VV}_{<}(|{\mathbf{k}}|,\nu_{\chi})=a^{NN}_{<}(|{\mathbf{k}}|)+a^{\pi\pi}_{<}(|{\mathbf{k}}|,\nu_{\chi})+\delta
a^{\pi\pi}_{<}(|{\mathbf{k}}|)\,,$ (6.39)

where the input functions are given in Eqs. (6.27), (6.3.2), and (6.37).

Concerning ${\cal A}_{VV}^{>}$, we use the OPE results given in Section 4.2
and Appendix B, assembled to produce the vector–vector linear combination, as
done in Section 6.2. This leads to

${\cal A}_{C}^{>,{\rm
sing}}\Big{|}_{VV}=\frac{3\alpha_{s}}{\pi}\,\int_{\Lambda}^{\infty}\
d|{\mathbf{k}}|\ \frac{1}{|{\mathbf{k}}|^{3}}\
\left(\frac{2(g_{1}^{NN}+g_{LR}^{NN})(\nu_{\chi})}{C^{2}}+\tilde{\cal
A}_{C}^{\rm sing}(\nu_{\chi})+\hat{\cal A}_{C}^{\rm
sing}(\nu_{\chi})\right)\,,$ (6.40)

where $\tilde{\cal A}_{C}$ is defined in Section 4.2,
$g_{LR}^{NN}(\nu_{\chi})$ and $\hat{\cal A}_{C}$ are defined in Appendix B,
and $\nu_{\chi}$ is the chiral EFT renormalization scale used in evaluating
$\langle pp|\hat{O}^{VV}|nn\rangle$. We chose to keep
$\nu_{\chi}\neq\mu_{\chi}$ for clarity. However, to simplify the analysis of
the IR divergences, we use here the same scale $\nu_{\chi}$ introduced in Eq.
(6.32) to separate out the $\log|\mathbf{p}|$ term from the low-$|\mathbf{k}|$
amplitude ${\cal A}_{C}^{\pi\pi,<,{\rm sing}}$. Given the scaling of the
$N\\!N$ couplings and the potentials that determine $\tilde{\cal A}_{C}$ and
$\hat{\cal A}_{C}$, discussed in Appendix B, in the above expression we keep
only the leading terms $g_{LR}^{NN}$ and $\hat{\cal A}_{C}$. The real part of
the singular component of $\hat{\cal A}_{C}$ is generated by pion exchange and
is given by

$\hat{\cal A}_{C}^{\rm
sing}(\nu_{\chi})=-\frac{m_{N}^{2}}{(4\pi)^{2}}\,\frac{g_{LR}^{\pi\pi}\,g_{A}^{2}}{4}\,\left[\log\frac{\nu_{\chi}^{2}}{4|{\bf
p}|^{2}}+1\right]\,,$ (6.41)

so that we can write141414Here and below the dependence on the QCD short-
distance renormalization scale $\mu$ in $\alpha_{s}$, $g_{LR}^{\pi\pi}$, and
$g_{LR}^{NN}$ is suppressed for simplicity. Note that, in contrast to
$g_{1}^{NN}$, $g_{LR}^{NN}$ carries an additional dependence on the chiral EFT
scale $\nu_{\chi}$.

$\displaystyle{\cal A}_{C}^{>,{\rm sing}}\Big{|}_{VV}$
$\displaystyle=\frac{3\alpha_{s}}{\pi}\,\int_{\Lambda}^{\infty}\
d|{\mathbf{k}}|\ \frac{1}{|{\mathbf{k}}|^{3}}\
\left(\frac{2g_{LR}^{NN}(\nu_{\chi})}{C^{2}}-\frac{m_{N}^{2}}{(4\pi)^{2}}\,\frac{g_{LR}^{\pi\pi}\,g_{A}^{2}}{4}\,\right)$
$\displaystyle-\frac{m_{N}^{2}}{(4\pi)^{2}}\,8g_{A}^{2}\,\left[\frac{3\,\alpha_{s}\,g_{LR}^{\pi\pi}}{16\pi}\,\int_{\Lambda}^{\infty}d|{\mathbf{k}}|\,\frac{1}{|{\mathbf{k}}|^{3}}\right]\,\log\frac{\nu_{\chi}}{2|{\mathbf{p}}|}\,.$
(6.42)

In the second line of the above equation, the expression within square
brackets is identified as $Z^{>}$, cf. Eq. (6.21). Moreover, using the
definitions (B.14) we can further simplify this expression to obtain:

$\displaystyle{\cal A}_{C}^{>,{\rm sing}}\Big{|}_{VV}$
$\displaystyle=2\,\frac{m_{N}^{2}}{(4\pi)^{2}}\,\left[-4g_{A}^{2}\ Z^{>}\
\log\frac{\nu_{\chi}}{2|{\mathbf{p}}|}+\
\int_{\Lambda}^{\infty}d|{\mathbf{k}}|\
a_{>}^{VV}(|{\mathbf{k}}|,\nu_{\chi})\right]\,,$ $\displaystyle
a_{>}^{VV}(|{\mathbf{k}}|,\nu_{\chi})$
$\displaystyle=\frac{3\alpha_{s}}{2\pi}\,\ \frac{(4\pi
F_{\pi})^{2}}{|{\mathbf{k}}|^{3}}\
\left(2\,\bar{g}_{LR}^{NN}(\nu_{\chi})-\frac{\bar{g}_{LR}^{\pi\pi}\,g_{A}^{2}}{4}\right)\,.$
(6.43)

Using the lattice QCD results of Ref. [90], we find
$\bar{g}_{LR}^{\pi\pi}=8.23$ in the $\overline{\rm MS}$ scheme at
$\mu=2\,\text{GeV}$. While both $\bar{g}_{LR}^{\pi\pi}$ and
$\bar{g}_{LR}^{NN}$ are expected to be $\mathcal{O}(1)$, given the large
numerical value of $\bar{g}_{LR}^{\pi\pi}$, for the nucleon coupling we will
assume $\bar{g}_{LR}^{NN}\in[-10,+10]$.

###  6.4 Matching

The matching condition for the vector amplitude ${\cal A}_{VV}$ leads to
(compare to Eq. (5.1) for the $LL$ amplitude)

${\cal A}_{C}^{<,{\rm sing}}\Big{|}_{VV}\ +\ {\cal A}_{C}^{>,{\rm
sing}}\Big{|}_{VV}\ =\ {\cal A}_{C}^{\rm sing}(\mu_{\chi})\Big{|}_{VV}\ +\
\frac{4({\cal C}_{1}(\mu_{\chi})+{\cal C}_{2}(\mu_{\chi}))}{C^{2}}\,.$ (6.44)

The input needed in Eq. (6.44) can be found in Eqs. (6.25), (6.26), (6.27),
(6.32), (6.3.2), and (6.43). Rescaling the couplings and integrals as in
Section 5 and using $Z=Z^{<}+Z^{>}$, one arrives at

$\displaystyle 2\,\bigg{(}\tilde{C}_{1}(\mu_{\chi})$
$\displaystyle+\tilde{C}_{2}(\mu_{\chi})\bigg{)}=1+2\,Z\,g_{A}^{2}-4Zg_{A}^{2}\,\log\frac{\nu_{\chi}}{\mu_{\chi}}-\int_{0}^{\mu_{\chi}}d|{\mathbf{k}}|\
a^{VV}_{\chi}(|{\mathbf{k}}|)$
$\displaystyle+\int_{0}^{\Lambda}d|{\mathbf{k}}|\
a^{VV}_{<}(|{\mathbf{k}}|,\nu_{\chi})+\int_{\Lambda}^{\infty}d|{\mathbf{k}}|\
a^{VV}_{>}(|{\mathbf{k}}|,\nu_{\chi})\,,$ (6.45)

with $a^{VV}_{<}(|{\mathbf{k}}|,\nu_{\chi})$ given in Eq. (6.39). Note that
the logarithmic IR dependence on the external momentum $|\mathbf{p}|$ has
disappeared in the matching relation, providing a strong consistency check on
the calculation. Except for $\delta a_{<}^{\pi\pi}(|{\mathbf{k}}|)$, which is
given in Eq. (6.37) and has to be obtained via numerical integration, the
quantities relevant to evaluate the matching condition (6.45) are

$\displaystyle a^{VV}_{\chi}(|{\mathbf{k}}|)$ $\displaystyle=-2\
\frac{1}{|{\bf k}|}\,\theta(|{\bf k}|-2|{\bf p}|)\,,$ $\displaystyle
a^{NN}_{<}(|{\mathbf{k}}|)$
$\displaystyle=-\left(2g_{V}^{2}({\mathbf{k}}^{2})+\frac{{\mathbf{k}}^{2}g_{M}^{2}({\mathbf{k}}^{2})}{m_{N}^{2}}\right)\,8\,{\cal
I}_{C}^{<}(|{\mathbf{k}}|)\,,$ $\displaystyle
a_{<}^{\pi\pi}(|{\mathbf{k}}|,\nu_{\chi})$ $\displaystyle=-4g_{A}^{2}\
\frac{|{\mathbf{k}}|}{(4\pi F_{0})^{2}}\Bigg{\\{}8\log
2-\frac{(\omega_{V}-|{\mathbf{k}}|)\big{(}2\omega_{V}^{2}-3|{\mathbf{k}}|\omega_{V}-|{\mathbf{k}}|^{2}\big{)}}{\omega_{V}^{3}}$
$\displaystyle+\frac{(\omega_{V}-|{\mathbf{k}}|)^{2}(2\omega_{V}+|{\mathbf{k}}|)}{\omega_{V}^{3}}\bigg{(}3\log\frac{\omega_{V}+|{\mathbf{k}}|}{\nu_{\chi}}-1\bigg{)}$
$\displaystyle-6\log\bigg{(}1+\frac{\omega_{V}}{|{\mathbf{k}}|}\bigg{)}-2\log\bigg{(}1+\frac{|{\mathbf{k}}|}{\omega_{V}}\bigg{)}\Bigg{\\}}\,,$
$\displaystyle a_{>}^{VV}(|{\mathbf{k}}|,\nu_{\chi})$
$\displaystyle=\frac{3\alpha_{s}}{2\pi}\,\ \frac{1}{|{\mathbf{k}}|^{3}}\
\left(2\bar{g}_{LR}^{NN}(\nu_{\chi})-\frac{\bar{g}_{LR}^{\pi\pi}\,g_{A}^{2}}{4}\right)\,.$
(6.46)

Figure 11:  Left panel: $a^{VV}_{\chi}(|{\mathbf{k}}|)$ (yellow),
$a^{VV}_{<}(|{\mathbf{k}}|,M_{\pi})$ (blue), and
$a^{VV}_{>}(|{\mathbf{k}}|,M_{\pi})$ (green) assuming
$\bar{g}_{LR}^{NN}\in[-10,10]$ (lower to upper green curves and shaded band).
Right panel: same plot for the absolute values in logarithmic scale.

Without loss of generality, in the numerical results we set
$\nu_{\chi}=M_{\pi}$. The various components of the integrand are shown in
Figs. 11 and 12. Fig. 11 displays the chiral, low- and high-$|{\mathbf{k}}|$
components of the integrand, namely $a^{VV}_{\chi}(|{\mathbf{k}}|)$ (yellow),
$a^{VV}_{<}(|{\mathbf{k}}|,M_{\pi})$ (blue), and
$a^{VV}_{>}(|{\mathbf{k}}|,M_{\pi})$, in both linear scale (left panel) and
logarithmic scale (right panel). Fig. 12 shows
$a^{\pi\pi}_{<}(|{\mathbf{k}}|)$ (blue), $\delta
a^{\pi\pi}_{<}(|{\mathbf{k}}|)$ (yellow), and their sum in green. Fig. 13
shows the dependence of $(\tilde{\cal C}_{1}+\tilde{\cal C}_{2})$ on the
matching scale $\Lambda$ (left panel) and the chiral renormalization scale
$\mu_{\chi}$ (right panel).

Figure 12:  $a^{\pi\pi}_{<}(|{\mathbf{k}}|,M_{\pi})$ (blue), $\delta
a^{\pi\pi}_{<}(|{\mathbf{k}}|)$ (yellow), and their sum (green).

We now discuss the uncertainties associated with our estimate of $(\tilde{\cal
C}_{1}+\tilde{\cal C}_{2})$: (i) Concerning the parametric uncertainties, the
dependence on the vector form factor parameterization is mild: we find an
upward shift of $\pm 0.03$ when using the Arrington–Sick fit [82], compared to
our baseline dipole parameterization. On the other hand, the uncertainty due
to the choice of $\Lambda\in[2\,\text{GeV},4\,\text{GeV}]$, coupled to the
range $\bar{g}_{LR}^{NN}|\in[-10,+10]$ (motivated by the large value of the
corresponding two-pion matrix element, see Appendix B), is at the level of
$\pm 0.3$, larger than in the case of $\tilde{\cal C}_{1}$, see Fig. 13 (left
panel). (ii) To estimate the systematic uncertainty due to the choice of the
short-range $N\\!N$ interaction $V_{S}$, we repeat the analysis using the Reid
and AV18 potentials, finding shifts in $(\tilde{\cal C}_{1}+\tilde{\cal
C}_{2})$ of $\approx+0.25$ and $\approx-0.3$, respectively. (iii) For the
systematic effect due to the neglected inelastic channels we take a range
$\delta(\tilde{\cal C}_{1}+\tilde{\cal C}_{2})\approx\pm 1.1$, which leads to
a relative error of about 50% in the singular $N\\!N$ electromagnetic
amplitude at $|\mathbf{p}|\sim 25\,\text{MeV}$, larger than the 30% in the
weak amplitude controlled by $\tilde{\cal C}_{1}$. In addition to the new
class of pion-exchange diagrams, this larger assignment is motivated as
follows: in contrast to the $\tilde{\cal C}_{1}$ analysis, the parametric
error is now more sizable, mainly driven by the $N\\!N$ short-distance
coupling $\bar{g}_{LR}^{NN}$. In fact, reducing $\Lambda$ to values as low as
$1\,\text{GeV}$ and thus into the energy region where the applicability of the
OPE becomes questionable and inelastic effects important, leads to a variation
$\delta(\tilde{\cal C}_{1}+\tilde{\cal C}_{2})\simeq\pm 1.0$. Since this
effect concerns the intermediate-momentum region which is most uncertain in
our analysis, the resulting variation could be either booked as an uncertainty
obtained by extrapolating the OPE expression beyond its region of validity, or
in terms of neglected intermediate states. We prefer to keep the OPE scale
$\Lambda\gtrsim 2\,\text{GeV}$, and thus increase the estimate of inelastic
contributions accordingly to account for the more prominent effect of the
short-distance coupling for $\tilde{\cal C}_{1}+\tilde{\cal C}_{2}$.

Altogether, our final result is

$(\tilde{\cal C}_{1}+\tilde{\cal C}_{2})(\mu_{\chi}=M_{\pi})\simeq
2.9(1.1)_{\rm inel}(0.3)_{V_{S}}(0.3)_{\rm par}=2.9(1.2)\,.$ (6.47)

Concerning the central value, the pion contributions ($2Zg_{A}^{2}$ plus the
integral of $a^{\pi\pi}_{<}(|{\mathbf{k}}|,M_{\pi})+\delta
a^{\pi\pi}_{<}(|{\mathbf{k}}|)$) amount to $+2.4$. Moreover, the most
uncertain intermediate momentum region
$|{\mathbf{k}}|\in[0.4,1.5]\,\text{GeV}$ contributes $\Delta(\tilde{\cal
C}_{1}+\tilde{\cal C}_{2})\simeq 0.55$, safely below our uncertainty estimate.
At the renormalization scale $\mu=4M_{\pi}$ our result becomes

$(\tilde{\cal C}_{1}+\tilde{\cal C}_{2})(\mu_{\chi}=4M_{\pi})=7.8(1.2)\,.$
(6.48)

Figure 13:  Left panel: dependence of $(\tilde{\cal C}_{1}+\tilde{\cal
C}_{2})(\mu_{\chi}=M_{\pi})$ on the matching scale $\Lambda$ using
$\bar{g}_{LR}^{NN}=\pm 10$. Right panel: dependence of $\tilde{\cal
C}_{1}+\tilde{\cal C}_{2}$ on the chiral renormalization scale $\mu_{\chi}$.

###  6.5 Charge-independence-breaking contribution to $\boldsymbol{N\\!N}$
scattering

The result (6.47) already compares quite well to the phenomenological
determination $(\tilde{\cal C}_{1}+\tilde{\cal
C}_{2})(\mu_{\chi}=M_{\pi})=5.0$ from Ref. [43]. However, since the contact
term is scale and scheme dependent, it is more appropriate to compare directly
observables calculated based on Eq. (6.47). We therefore focus on the CIB
contribution to $N\\!N$ scattering. We summarize here our main findings and
report technical details of the analysis in Appendix E.

We first note that within LO chiral EFT the scattering lengths $a_{nn}$,
$a_{np}$, and $a_{pp}^{C}$ (the latter defined in the modified effective range
expansion to account for Coulomb effects [97, 98, 99]) can be mapped onto
contact terms for each channel

$\tilde{C}_{np}=\tilde{C}+\frac{e^{2}}{3}\left(\tilde{\cal C}_{1}+\tilde{\cal
C}_{2}\right)\,,\qquad\tilde{C}_{nn/pp}=\tilde{C}-\frac{e^{2}}{6}\left(\tilde{\cal
C}_{1}+\tilde{\cal C}_{2}\right)\pm\frac{1}{2}\tilde{\cal C}_{\text{CSB}}\,,$
(6.49)

where $\tilde{C}$ denotes the isospin-symmetric combination, $\tilde{\cal
C}_{1}+\tilde{\cal C}_{2}$ the CIB contribution, and $\tilde{\cal
C}_{\text{CSB}}$ a charge-symmetry-breaking (CSB) term [100]. To test our
prediction for $\tilde{\cal C}_{1}+\tilde{\cal C}_{2}$, we can thus use two
observables to determine $\tilde{C}$ and $\tilde{\cal C}_{\text{CSB}}$, and
then predict the third based on Eq. (6.47). Several choices are possible (see
Appendix E), of which arguably the simplest is

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

This combination of scattering lengths would isolate the CIB contribution if
$N\\!N$ scattering were perturbative and Coulomb interactions absent, but in
practice depends on all couplings in a complicated manner. To obtain the
numerical result we have used the empirical values
$a_{pp}^{C}=-7.817(4)\,\text{fm}$ [101, 102], $a_{np}=-23.74(2)\,\text{fm}$
[103, 104], $a_{nn}=-18.9(4)\,\text{fm}$ [105]. From Eq. (6.47) we find
$a_{\text{CIB}}=15.9^{+4.5}_{-4.0}\,\text{fm}$, in good agreement with Eq.
(6.50), given that additional uncertainties from higher chiral orders could be
attached. The comparison to the phenomenology of CIB in $N\\!N$ scattering
thus validates our approach at the level of (30–50)% and shows that our
uncertainty estimates are realistic.

##  7 Synthetic data for $\boldsymbol{nn\to pp}$ near threshold

Having determined the LEC $\tilde{\cal C}_{1}$ in the $\overline{\text{MS}}$
scheme, we can now compute the low-energy $nn\rightarrow ppe^{-}e^{-}$
amplitude $\mathcal{A}_{\nu}$. Our matching strategy was based on dimensional
regularization with minimal subtraction as it provides convenient and
factorized expressions for the amplitude. While dimensional regularization is
rarely used by nuclear practitioners, this is no obstacle to applying our
results in nuclear-structure calculations. Observables, such as
$\mathcal{A}_{\nu}$, are scheme independent so that the LNV contact term can
by obtained in any scheme through a fit to our synthetic data for the
amplitude

$\mathcal{A}_{\nu}(|{\mathbf{p}}|,|{\mathbf{p}}^{\prime}|)=-\langle\Psi_{pp}(|{\mathbf{p}}^{\prime}|)|\
V_{\nu\,\rm L}^{{}^{1}S_{0}}+V_{\nu\,\rm S}^{{}^{1}S_{0}}\
|\Psi_{nn}(|{\mathbf{p}}|)\rangle\,,$ (7.1)

where both initial $|\Psi_{nn}(|{\mathbf{p}}|)\rangle$ and final
$|\Psi_{pp}(|{\mathbf{p}}^{\prime}|)\rangle$ states are in the ${}^{1}S_{0}$
channel, $V_{\nu\,\rm L}^{{}^{1}S_{0}}$ denotes the usual long-range neutrino
potential, and $V_{\nu\,\rm S}^{{}^{1}S_{0}}$ is the short-range interaction
proportional to $\tilde{\cal C}_{1}$. We denote by $E={\mathbf{p}}^{2}/m_{n}$
and $E^{\prime}={\mathbf{p}}^{\prime 2}/m_{p}$ the center-of-mass energies of
the incoming neutrons and outgoing protons of masses $m_{n}$ and $m_{p}$,
respectively, and by ${\mathbf{p}}$ and ${\mathbf{p}}^{\prime}$ the
corresponding relative momenta. Assuming the outgoing electrons to be at rest
one can determine the maximum momentum carried by the outgoing protons, given
the incoming momentum of the neutrons $|{\mathbf{p}}|$, via

$E^{\prime}=E+2(m_{n}-m_{p}-m_{e})\,,\qquad|{\mathbf{p}}^{\prime}|=\sqrt{{\mathbf{p}}^{2}+2m_{N}(m_{n}-m_{p}-m_{e})}\,,$
(7.2)

with $m_{e}$ the electron mass and $2m_{N}=m_{n}+m_{p}$. If the electrons are
not at rest but carry total zero momentum in the incoming neutrons’ rest
frame, the outgoing protons fly back-to-back with momentum
$|{\mathbf{p}}^{\prime}|$ ranging between zero and the maximum value given in
Eq. (7.2).

In Ref. [59] we chose the momenta $|{\mathbf{p}}|=25\,\text{MeV}$ and
$|{\mathbf{p}}^{\prime}|=30\,\text{MeV}$, which proved an advantageous
kinematic point for which the LO chiral amplitude provides an accurate
prediction, while staying away form the sizable isospin-breaking effects at
threshold. In practice we compute the amplitude in coordinate space151515The
relation between ${\cal C}_{1}$ and $\tilde{\cal C}_{1}$ is given in Eq.
(5.2).

$\mathcal{A}_{\nu}(|{\mathbf{p}}|,|{\mathbf{p}}^{\prime}|)=-\int
d^{3}{\mathbf{r}}\,\psi^{-\,*}_{{\mathbf{p}}^{\prime}}({\mathbf{r}})\left(V^{{}^{1}S_{0}}_{\nu\,\rm
L}({\mathbf{r}})-2\,{\cal
C}_{1}\,\delta^{(3)}({\mathbf{r}})\right)\psi^{+}_{{\mathbf{p}}}({\mathbf{r}})\,,$
(7.3)

with the long-range neutrino potential given by

$V^{{}^{1}S_{0}}_{\nu{\rm L}}(r)=\frac{1+2g_{A}^{2}}{4\pi
r}+\frac{g_{A}^{2}}{4\pi
r}\left[1-e^{-M_{\pi}r}\left(1+\frac{M_{\pi}r}{2}\right)\right]\,,$ (7.4)

and wave functions for the scattering states