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

We present a method to determine the leading-order (LO) contact term contributing to the $nn \to pp e^-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 $^1S_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 $^1S_0$ $N\! N$ scattering lengths. While our analysis is performed in the $\overline{\rm MS}$ scheme, we express our final result in terms of the scheme-independent renormalized amplitude ${\cal A}_\nu(|{\bf p}|,|{\bf p}^\prime|)$ at a set of kinematic points near threshold. We illustrate for two cutoff schemes how, using our synthetic data for ${\cal 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.


Introduction
Neutrinoless double β decay (0νββ) 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νββ 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 ββ = i U 2 ei m i , where the sum runs over light neutrino masses m i and U ei are elements of the neutrino mixing matrix. 0νββ 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 β 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 ∆L = 2 effective neutrino mass m ββ . When considering the 1 S 0 channel, in analogy to the nucleon-nucleon (NN ) 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 → 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 β 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 NN 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 → 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 [46,[52][53][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νββ decays besides the nuclear-structure ones, especially given that its impact is enhanced in ∆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 → ppe − e − amplitude including this contact-term contribution [59].
The hadronic component of the light-Majorana-neutrino-exchange amplitude has the structure 1) and is controlled by the two-nucleon matrix element of the time-ordered product T {j α w (x)j β w (0)} of two weak currents. Such matrix elements with the weak current replaced by the electromagnetic current j α em (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νββ decay nn → 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 protonneutron mass difference m el p−n = 0.75 (2) MeV. Despite the tension between the estimate of the inelastic contributions in lattice QCD, m inel p−n = 0.28 (11) MeV [75][76][77], and from nucleon structure functions, m inel p−n = −0.17 (16) MeV [67,[72][73][74], also in this case the elastic estimate is accurate at the 30% level. The main complication in the generalization to 0νββ decay arises from the twoparticle 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 NN 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 NN scattering amplitude. To validate this approach, we also consider the twonucleon matrix element with two electromagnetic currents, which can be accessed experimentally in terms of charge independence breaking (CIB) in the NN 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νββ 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 → 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 NN scattering lengths. This section is fairly technical and can be omitted by readers primarily interested in 0νββ decay. In Section 7 we return to the LNV amplitude nn → 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 NN 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 NN 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 µ Λ χ ∼ 1 GeV is given by 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 ββ = i U 2 ei m i in terms of mass eigenstates and elements of the neutrino mixing matrix, 1 and a dimension-nine ∆L = 2 operator generated at the electroweak threshold: with e c L = Cē T L . Since C L = (8V 2 ud G 2 F m ββ )/M 2 W × (1 + O(α s /π)), the effect of the latter term on the 0νββ amplitude is suppressed by (k F /M W ) 2 (where k F ∼ O(100) MeV is the typical Fermi momentum of nucleons in a nucleus) compared to light-neutrino exchange and can be safely 1 The effective mass probed in 0νββ decay is often defined as m ββ = i U 2 ei mi , but to simplify the notation at the Lagrangian level we formally keep its phase. 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 ∆L = 2 transitions (such as π − π − → e − e − , nn → ppe − e − , 76 Ge → 76 Se e − e − , 136 Xe → 136 Ba e − e − , . . . ) through the non-local effective action obtained by contracting the neutrino fields in the two weak vertices, is the scalar massless propagator, a remnant of the neutrino propagator. Computing matrix elements of S ∆L=2 eff 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 2 and p ext denotes generically the hadronic external momenta p f and p i . Using translational invariance one has one arrives at The hadronic amplitude A ν 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 µν (k, p ext ) will be thought of as the generalized forward Compton amplitude (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 where 2E = 2E + 2m e . Free two-nucleon states with vanishing total three-momentum and individual three-momenta given by ±q will be denoted by |q [23], so for example for the initial and final state we will have |i 0 = |p and |f 0 = |p , respectively.

Matching strategy
The amplitude for the process nn → 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 T (k, p ext ) = g µν T µν . 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 2 , and low-energy (k 2 E < Λ 2 χ ) regions, separated by Λ χ (the breakdown scale of the low-energy hadronic EFT) and Λ (scale at which the OPE becomes reliable). The low-energy region further includes the soft (|k 0 | ∼ |k| ∼ k F ), potential (|k 0 | ∼ k 2 F /m N , |k| ∼ k F ), and ultrasoft (|k 0 | ∼ |k| 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.
In practice, given the non-relativistic nature of the process of interest, we will not use k 2 E as matching variable. Instead, we decompose d 4 k = dk 0 d 3 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 k to reduce the amplitude to an integral over d|k|. To LO in the expansion in external momenta |p|, |p | ∼ Q (we denote by µ χ ∼ Q the soft scale in the EFT), we write the full amplitude as an integral over the internal neutrino three-momentum k, which we split into a low-plus

intermediate-momentum region and a high-momentum region
separated by the scale Λ that represents the onset of the asymptotic behavior for the currentcurrent correlator, controlled by the OPE. This representation introduces model dependence through: (i) The choices made to extend the model-independent integrand a χ (|k|) dictated by chiral EFT in the region |k| < Λ χ to the function a < (|k|) valid up to |k| ∼ Λ. We will provide a simple parameterization of a < (|k|) that reduces to a χ (|k|) for |k| < M π 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 Λ that determines the boundary of integration regions in the variable |k|. Once a representation for A full ν is obtained, along with an estimate of the associated uncertainties, we will estimate the LEC appearing in A χEFT ν by enforcing the matching condition (2.13) In the following sections, we will describe the construction of A <,> and the matching to A χEFT ν , starting with the spectral representation in Section 2.3.

Spectral representation
In this section we provide a spectral representation for the nn → pp amplitude, which will prove very useful in identifying and organizing the various intermediate-state contributions to 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Ĥ =Ĥ 0 +V the total Hamiltonian, split into a free and interaction term (V , not to be confused with the potential). Retarded and advanced Green's functions in the interacting theory are given bŷ with analogous definitions for the free-theory ones, denoted byĜ ± (E), with the replacement H →Ĥ 0 . The scattering operatorT (E) is formally given bŷ The scattering states f − | and |i + are related to the free states f 0 | and |i 0 by In terms of the scattering states, the amplitude for nn → pp can be written as with the weak transition operator From the definition of the correlator in Eq. (2.6) one obtains the following representation for Π LL µν (k) in terms of Green's functions: whereP is the total three-momentum operator and we have introduced the four-vectors The labels in p µ i,f refer to the initial and final states between whichΠ LL µν (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 havep = 0 and hence k ± = ±k. In Eq. (2.19) the dependence on k 0 is very simple, as k 0 appears only through the energy denominators of Ĝ + (k 0 ± ). Performing the integration over k 0 in Eq. (2.18) with Cauchy's theorem, 4 one arrives atÔ Further inserting a complete set of states between the current operators in Eq. (2.21) leads to the spectral representation for the amplitude 5 (2.22) The representations (2.17) and (2.22) are quite general. The asymptotic behavior of the integrand in Eq. (2.22) at large |k| is dictated by the OPE forΠ LL µν (k) or, equivalently,Ô LL (k). An explicit calculation to be described below shows the behavior d 3 k/|k| 5 , so the amplitude in the full theory is finite. Moreover, Eq. (2.22) shows that once |k| > k F , so that k 2 /m N is above the typical nuclear binding energies, one expects (E n (k ± ) −Ẽ) > 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 |k|. The matrix elements in the numerator are also expected to have a smooth behavior in |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 |k| Λ χ and |k| Λ 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 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). 5 The summation is over intermediate states |n(k±) of total three-momentum k±, enforced by the δ-functions in Eq. (2.21). Therefore, for an N -particle intermediate state n involves phase space integrals over the N − 1 internal momenta (the total momentum being fixed to k±) and carries non-zero mass dimension. For example, for two-nucleon intermediate states, using non-relativistic normalizations for the states pn|p n = (2π) 3 δ (3) (pn − p n ) one has n → d 3 pn/(2π) 3 , where pn is the relative momentum of the two-nucleon pair. In general the summation n |n(k±) n(k±)| carries mass dimension −3.
(0) m denotes the energy associated with the free HamiltonianĤ 0 , and To LO in chiral EFT the amplitude A ν 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/|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 m |Ô LL (k)|m , p m |T (E)|p , and p |T (E )|p m 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 NN π inelastic channel as we expect this to be one of the dominant sources of uncertainty in our final result.

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 A ν in chiral EFT is described in Refs. [40,43] and we recall here some of its elements as needed. Denoting by µ χ ∼ Q the soft scale in the EFT, to LO in chiral counting, i.e., 1/Q 2 , only elastic NN intermediate states are relevant. The corresponding diagrams are reported in Fig. 2. For concreteness, we regulate all the integrals dimensionally and perform MS subtraction of the divergences when needed. The LO chiral EFT results correspond to replacing in Eq. (2.23) the LO form for theÔ LL (k) andT (E) operators, denoted byÔ LL χ (k) andT χ (E), respectively, and using the non-relativistic form of the free two-nucleon Green's function, [G (0) . The current-current correlator to LO in chiral EFT is given by (3.1) TheT operator is determined by the LO interaction HamiltonianV =V π +V S , which contains the one-pion-exchange and a short-range contribution, parameterized to LO by the LEC C: Here F π is the pion decay constant and M π denotes the pion mass. The LO 1 S 0 NN contact coupling C contains contributions from one-pion exchange as well as a contact interaction and scales as C ∼ 4π/(m N Q). The splitV =V π +V S implies that theT operator can be similarly separated into a pion-range and short-range contribution as follows [80] T whereĜ (π) 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 T χ S , the LO chiral EFT amplitude can be written as contact interaction, in dimensional regularization one can show that the rescattering factors and Yukawa wave functions at the origin χ + p (0), χ + p (0) can be factored out of the d 3 k integrals in M χ B+B , and M χ C , thus reproducing the chiral EFT results of Refs. [40,43] The divergence in G + E (0, 0) is absorbed by C −1 , so that K E is well defined and independent of the chosen scheme and scale [23].
The reduced amplitudes 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 thê V S insertion). Neglecting Yukawa interactions (α π → 0), i.e., in the pionless EFT (/ πEFT), 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 A C in Eqs. (3.5) and (3.8) contains a UV-divergent term at the two-loop level, which we denote by A sing C , as well as UV-finite terms induced by pion exchange, which we denote by δA C , leading to the decomposition: The UV-convergent term δ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 ofT χ S ); (ii) the term proportional to g 2 A M 4 π inÔ LL χ (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 2 A → 2g 2 A ). Using dimensional regularization with scale µ χ and the MS scheme for renormalization, we thus identify the singular (UV-divergent) term with where In case of p = p the resulting expression becomes It will prove useful to write the real part of this result in an integral representation The UV divergence is removed by introducing a new contact term with a LNV coupling proportional to m ββ denoted by C 1 [43], to be identified with g N N ν also used in the recent literature, i.e., C 1 ≡ g N N ν . 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 The matching analysis will provide a representation for the combination A sing C (µ χ )+2 C 1 (µ χ )/C 2 .

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.,T =T π +T S , as done in Eq. (3.3).

Low-and intermediate-momentum region: A <
Following Eqs. (2.23) and (3.4), we write We wish to provide a representation of the integrand in M < A,B,C that is valid up to a scale |k| ∼ Λ. We build the integrand starting from the low-momentum region. For |k| ∼ Q the integrand is controlled by chiral EFT. Using the LO chiral representation for the two-current operator (Ô LL < =Ô LL χ ) and the T matrix (T < π,S =T χ π,S ) in M A and M B+B leads to convergent integrals, and one recovers the chiral EFT results of Eqs. (3.5), up to terms of O(Q/Λ) that are irrelevant at LO We will also show that M > A,B+B does not contribute to the amplitude at LO. Therefore, M A and M B+B drop out of the LO matching condition (2.13), which will therefore only involve M < C . In order to have a representation of the integrand valid up to momenta |k| ∼ Λ, in Eq. (4.2) we need to model bothÔ LL (k) andT 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.
For the current-current operator we make the replacementÔ LL χ (k) →Ô LL < (k), by introducing single-nucleon form factors in the weak current vertices, where In our baseline analysis, we will use the simple dipole parameterization for the vector, magnetic, and axial-vector form factors with Λ V = 0.84 GeV and Λ A = 1.0 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 Λ A ∈ [0.8, 1.2] GeV, which corresponds to an axial radius of 0.47 (19) fm 2 , in good agreement with the result 0.46(16) 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 [81,[83][84][85]. The parameterization ofÔ LL < (k) in terms of the weak nucleon form factors then leads to where the notation for the HOST S -matrix elements is 6 We further write the HOST S -matrix element as the product of the on-shellT S -matrix element times the HOS factorsf S (q, p) [86] Using the LO on-shell result [23], see Eq. (3.3), and definingf we arrive at Equation (4.12) shows that M < C has the correct IR behavior to LO in the external momenta p and p and that the quantity A < C corresponds to the amplitude A C in chiral EFT. As argued in Section 3, A C contains a singular UV-divergent term (A sing C ) and UV-finite contributions induced by pion exchange (δA C ). In analogy to Eq. (3.9) we can write where δA < C denotes the convergent pion-exchange contributions to A < C , satisfying 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 A < C . A <,sing C is obtained from Eq. (4.13) by systematically removing the convergent pion contributions. This requires: 1. Discarding the convergent terms due to the induced pseudoscalar form factor in h LL GT (k 2 ). In practice this means replacing 3h LL 6 These are S-wave projected matrix elements and depend only on |p| and |q|. . The various lines correspond to the NLO / π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(|k|) = 1.
2. Setting α π → 0 in the evaluation off S (q, p) in Eq. (4.14) , which implies The quantity f S (q, 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 In the end, the difference between A sing C (µ χ ) and A <,sing C is quite intuitive. As illustrated in Fig. 3, the integrand of A sing C in Eq. (3.10) contains the LO nucleon weak current vertices and the leading contact NN interaction. On the other hand, the full-theory integrand in A <,sing C contains appropriate form factors in the nucleon weak current vertices and in the NN vertices (f S (q, p)). The latter parameterizes the off-shell behavior of the short-range component of the NN amplitude and amounts to changing the LO chiral EFT function I C (|k|) defined in Eq. (3.11) to I < C (|k|) given in Eq. (4.18), valid for a larger range of momenta |k| compared to the LO result. It will prove convenient to present results in terms of the ratio .
We now describe our construction of I < C (|k|) and r(|k|), based on the first two orders in chiral EFT and then extended to larger values of |k| using models for the short-range NN interaction in the 1 S 0 channel. For concreteness we work at the kinematic point |p| = |p | = 1 MeV.
1. At LO in chiral EFT one has f S (q, p) = 1 and the integrals simplify to Re 2. For |k| ≤ M π and |k| ≤ Λ χ we can evaluate the corrections to Re I < C (|k|) in / πEFT and chiral EFT, respectively, using V S (p, p ) = C + (C 2 /2)(p 2 + p 2 ). As discussed in Appendix A, the results at next-to-leading order (NLO) in chiral EFT and / πEFT are formally identical with the identifications C = C 0 in / πEFT and C = C 0 + g 2 A /(4F 2 π ) in chiral EFT. Note that these HOS form factors are different depending on whether they involve the initialor final-state on-shell momenta, so that f S (p , q) cannot be obtained from the expression of f S (q, p) by interchanging q ↔ p and subsequently replacing p → p . This leads to (see Appendix A for details) In terms of the ratio r(|k|) introduced in Eq. (4.19) the NLO analysis gives In both EFTs the ratio Since r 0 1/(72 MeV) (corresponding to C 2 /C 2 0.52), this produces sizable fractional deviations in Re I < C (|k|) already for |k| ≤ M π . In chiral EFT one finds a slightly reduced suppression C 2 /C 2 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 |k|, linked to the large effective range in the 1 S 0 channel.
3. To extend I < C (|k|) or equivalently r(|k|) to higher values of |k|, we have computed f S (q, p) using models of V S , i.e., the short-range NN 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 |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 (q, p), which was then used to evaluate I < C (|k|) (see Appendix A for details). Reassuringly, as illustrated in Fig. 4, the behavior of r V S (|k|) in the model calculation very closely tracks the NLO chiral EFT result for |k| < 200 MeV, i.e., the region in which NLO chiral EFT is expected to provide modelindependent and accurate results. Using other NN potentials [88,89] does not change the qualitative picture, but induces small changes in the slope of r(|k|) for |k| < 400 MeV as well as the location and depth of the minimum for |k| ∼ 600 MeV. This makes it clear that using r V S (|k|) all the way to |k| ∼ Λ ∼ 1 GeV introduces model dependence. However, in the integral (4.18) the region in which r V S (|k|) has the largest model dependence is weighted Figure 5: Quark-level diagrams that determine the asymptotic behavior ofΠ LL µν . Solid lines represent quarks, curly lines represent gluons, and the black dots represent insertions of the currentū L γ µ d L .
considerably less than the model-independent low-|k| region, because the integration kernel involves 1/|k| and the nucleon form factors, both rapidly decreasing functions of |k|. In Section 5 we will quantify how this model dependence affects the extraction of the effective coupling.

High-momentum region:
and T a are SU (3) color generators.
In the actual matching calculation, since we contract T µν with g µν , we need where the last step holds under symmetric integration (k α k β → k 2 g αβ /4). Finally, we note that by using Fierz identities in color and then in flavor one obtainŝ Performing the k 0 integration in order to go from g µνΠLL µν (k, 0) toÔ LL (k) we obtain which is the expression to be used to evaluate the A > component in Eq. (2.23), to obtain This result shows a factorization of the hard and soft contributions. The final step requires the evaluation of f − | O 1 (0) |i + , 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 The scaling of the non-perturbative parameters is g N N 1 , g πN , while a precise lattice calculation gives g ππ 1 = 0.36 (2) at the renormalization scale µ = 2 GeV in the MS scheme [90]. 7 The ππ and πN couplings induce a pion-range transition operator, while the NN couplings gives a short-range transition operator. The final result is given by andÃ A,B,C have the same formal expression as A A,B,C [43] with the replacement [42] (4.33) These relations allow us to write A > in a form very useful for the matching to chiral EFT: The quantities A > A,B,B are sub-leading compared to their finite chiral EFT counterparts A A,B,B , e.g., A > A /A A ∼ Q 2 /Λ 2 and similarly for A B . Therefore, as anticipated in the previous section, these quantities do not enter the LO matching relation. 7 See Refs. [91,92] for a determination of g ππ 1 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 µ, needed to cancel the scale dependence in αs(µ). 8 Note that the mismatch in dimensions between V On the other hand, A > C has to be retained because its integrand provides the QCD asymptotic behavior to which the integrand in A <,sing C given in Eq. (4.18) has to tend for large |k|. Given our limited knowledge of g N N ,πN 1 and the fact that the terms involving g N N 1 , g πN 1 , and g ππ 1 contribute at the same order to f − | O 1 (0) |i + , we will retain for simplicity only the term proportional to g N N 1 . For convenience, in the matching analysis we express g N N 1 as follows In numerical estimates we will allow forḡ N N 1 as large as O (10), given that the ππ short-distance couplingḡ ππ LR = 8.23 [90] relevant for the vector-vector correlator by far exceeds its 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 A > C is to enforce the UV finiteness of the nn → pp amplitude.
5 Matching and extraction of C 1

Results
The matching condition (2.13), taking into account the form of A χEFT ν in Eq. (3.15) and the structure of A < (see Eq. (4.1) and following discussion) and A > (see Eq. (4.34) and following discussion) reduces to The amplitudes A sing C (µ χ ) and A <,sing C are given in Eqs.
and the dimensionless quantitiesĀ the matching condition reads Using our previous results, we can write the various quantities entering the matching condition in terms of the following integrals: , with r(|k|) given in Eq. (4.19). Putting all pieces together we can write By construction the µ χ dependence is consistent with the renormalization group equation (RGE) for the rescaled couplingC 1 (µ χ ) [43]. This representation shows that up to an additive constant the LECC 1 (µ χ ) can be thought of as the difference between two integrals in |k|, one in the full theory extending all the way to |k| → ∞, and one in the EFT extending to |k| = µ χ . Therefore, the LECC 1 (µ χ ) corresponds to (i) a possible mismatch between the LO chiral EFT and the full amplitude at |k| < µ χ ; (ii) the component of the full amplitude arising from |k| > µ χ . In Fig. 6 we show a χ (|k|) (yellow), a < (|k|) (blue), and a > (|k|) (green) assumingḡ N N 1 ∈ [−10, +10]. The behavior of the integrand indicates that the LEC is dominated by the lowand intermediate-momentum regions. In the left panel of Fig. 7 we show the dependence of C 1 (µ χ = M π ) on the separation scale Λ, which proves to be relatively mild. The right panel shows the dependence ofC 1 on the chiral renormalization scale µ χ . The impact of varyinḡ g N N 1 ∈ [−10, +10] is illustrated by the two different curves in the left panel of Fig. 7 and is very small.
Overall, the result for the LECC 1 (µ χ = M π ) +1.32 is relatively stable and confirms the expected scaling of the contact term, i.e.,C 1 ∼ O(1) [40]. This matching analysis implies that the contribution fromC 1 and the long-range neutrino exchange encoded in 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 νL = (1 + 2g 2 A )/q 2 > 0 has the opposite sign compared to the short-range potential V νS = −2C 1 < 0, since C 1 (µ χ ) > 0 for values of µ χ 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 |k| < µ χ , an effect quite visible in the doubly-logarithmic plot in the right panel of Fig. 6. In turn, the sizable deviation of a < (|k|) (blue) from a χ (|k|) (yellow) at |k| < µ χ = M π is due to the behavior of r(|k|) (see Fig. 4 and discussion in Section 4.1), which enters multiplicatively in a < (|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(|k|) for |k| < 300 MeV is controlled by the large 1 S 0 effective range r 0 2.7 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 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.

Discussion of uncertainties
Our baseline resultC 1 (µ χ = M π ) +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ḡ N N 1 controlling the high-|k| tail, and the matching scale Λ, as long as Λ ≥ 1 GeV. The couplingḡ N N 1 (µ), expected to be O(1), is presently unknown, but in view of the large value of the similar two-pion matrix elementḡ ππ LR = 8.23 (see Appendix B for details), we take the rangeḡ N N 1 ∈ [−10, +10], with minor impact on the final result. These parametric uncertainties do not exceed δC 1 ±0.05 par and are dominated by the input Λ A = 1.0(2) GeV.
The main systematics of our approach are due to the following effects: 1. We have effectively truncated the spectral representation (2.22) to keep only the elastic NN intermediate state. In order to assess the size of the neglected terms, we provide below a rough estimate of the simplest inelastic channel, namely NN π, 10 see Fig. 8. As discussed in Appendix C, we find contributions toC 1 from the NN π intermediate states on the order of |δC 1 | = 0.1-0.35 and due to this we assign an uncertainty of δC 1 | inelastic = ±0.5 .
2. Within the NN channel, we have built our integrand in |k| starting from the low-|k| end.
To extend the form of the integrand to |k| ∼ Λ ∼ 1 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 (q, p) that encode the higher-momentum behavior of the NN scattering amplitude. This effect, parameterized by the ratio r(|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(|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 |k| ≤ M π . On the other hand, the extension of r(|k|) to |k| ∼ Λ ∼ 1 GeV introduces model dependence. To quantify this model dependence we have employed both the simple three-resonance model and the AV18 model for the shortrange NN interactions in the 1 S 0 channel. In both cases the behavior of r(|k|) shows the same features, see Fig. 17: a rapid drop controlled by the large effective range, a zero for |k| between 400 and 500 MeV, and then a minimum around (600-700) MeV. However, the key point is that because of the fall-off of the integrand multiplying r(|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 onC 1 of ≈ +0.20 and ≈ −0.20 when using the Reid or AV18 potential, respectively. This difference correlates with the different slopes of r(|k) for |k| < 300 MeV, see Fig. 17. Based on this, we will assign an error of δC 1 of ±0.2 due to the choice of the short-range potential V S .
3. The above discussion and Fig. 6 both point to the fact that the bulk ofC 1 is controlled by the behavior of a < (|k|) for |k| < 300 MeV, where our integrand can be linked to the model-independent chiral EFT behavior. Quantitatively, we find that the region |k| ∈ [0.4, 1.5] GeV, entailing the largest unknowns, contributes ∆|C 1 | ≤ 0.05, which is a reassuring result.
4. Finally, we note that our approach can be used to estimate the combination of LECsC 1 +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 forC 1 +C 2 compares within uncertainties fairly well with the value extracted from the CIB combination of NN 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.
In light of the above discussion, we quote the following estimate for the LEC which becomesC 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Ā sing C +2C 1 at |p| ∼ (20-30) MeVin 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 A χEFT ν , as we will discuss in Section 7.
6 Vector-vector amplitude and C 1 + C 2 In analogy to the matching for the purely left-handed coupling C 1 = g N N ν , one can determine the left-right coupling C 2 and the vector combination C 1 + C 2 . Since C 1 + C 2 can be extracted from a fit to the I = 2 NN 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 → pp transition in the unphysical theory in which the W boson couples to both the lefthanded quark currentū L γ µ d L (as in the Standard Model) and the right-handed currentū R γ µ d R with equal couplings.

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 where q denotes the quark doublet q = (u, d) T and we defined where τ + = (τ 1 + iτ 2 )/2. We are interested in the isovector components of the interactions, involving l µ and r µ . We will consider the unphysical theory in which the W boson has vectorlike couplings and hence couples with same strength to bothū L γ µ d L andū R γ µ d R . This amounts to r µ = l µ , namely Double insertions of the isovector component of Eq. (6.1) give rise to 0νββ amplitudes and I = 2 electromagnetic effects through the effective actions: The interaction-specific factors are 11 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νββ 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: where Q a L = u † t a u , Q a R = ut a u † , (6.10) and u 2 = U = exp(iτ · π/F π ) incorporates the pion fields.
In the pion sector we have non-derivative LO local operators arising from the LR current correlator: 11) where, at LO in chiral perturbation theory (ChPT), Z is related to the pion-mass splitting by 12 6.2 Pion two-point function and the low-energy constant Z Before discussing the nn → pp amplitude, it is instructive to consider the simpler π − → π + transition. In complete analogy with Eq. (2.6) one can define a vector-vector correlator via the The object of interest for our matching calculation is T ππ (k, 0) determines the LEC Z defined in Eq. (6.11), while T ππ (k, p) (where p is an off-shell four-momentum with potential scaling) contributes to nn → pp as discussed below.
In the small p, k regime T ππ (k, p) can be computed reliably in ChPT. At LO one has 13 14) which in the chiral limit and for p → 0 reduces to The expression for T ππ (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 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 V π (k 2 ), for which we will take the simple monopole form 12 Note that the operator in the first line of Eq. (6.11) differs from the usual definition in ChPT [93] by an inessential constant. 13 Up to a proportionality factor due to isospin, this expression agrees with the forward Compton amplitude, see for example Ref. [66].
In the large k regime, the form of T ππ (k, p) is dictated by the OPE, see Section 4.2 and Appendix B. It is given by To LO in ChPT, the pion matrix elements are given by where the LECs g ππ 1,4,5 are known from lattice QCD [90,92] and scale as g ππ 1 ∼ O(1) and g ππ 4,5 ∼ O(Λ 2 χ ). From this scaling it is clear that the dominant contribution to the two-pion matrix element is proportional to the combination g ππ LR ≡ (1/2)(g ππ 4 − 3g ππ 5 ). 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 In the spirit of the matching strategy used in this work, we split up the integration over three-momentum |k| in two regions and by using T ππ < (k, 0) from Eq. (6.15) and T ππ > (k, 0) from Eqs. (6.18)-(6.19) we arrive at (see Appendix D for details) where ω V = k 2 + M 2 V and µ is the QCD renormalization scale. The above result is in agreement with estimates of the pion electromagnetic mass splitting [62][63][64][65][66]. Taking M V = 775 MeV, F π = 92.28 MeV, and g ππ LR = 8.23 from lattice QCD [90], the sum of low-and high-momentum components Z < (Λ) + Z > (Λ) equals 0.63 at Λ = 2 GeV and reaches the asymptotic value Z < (Λ → ∞) = 0.67. The deficit compared to the experimental value Z 0.8 is understood in terms of the neglected inelastic contributions from the axial-vector resonances [62][63][64], thus providing another estimate of the error due to neglecting inelastic corrections.

nn → pp vector-like amplitude in chiral EFT and full theory
Let us now consider the amplitude A V V in the unphysical theory in which the W boson has vector-like couplings to quarks. Denoting by A LL,LR the nn → pp amplitudes generated by Wexchange 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 Note that A LL ≡ A ν , i.e., the amplitude in the physical theory, defined in Eq. (2.9), and the factor of 2 in front of A LL arises from the fact that LL and RR products, both present due to vector-like couplings, give the same result by parity.

A V V in chiral EFT
At LO in chiral EFT, the nn → pp vector amplitude takes the form 23) where A N N LL is given in Eq. (3.15). A N N LR can be obtained from A N N LL by flipping the sign of the "A × A" axial contribution, which in practice amounts to setting g 2 A → −g 2 A everywhere and replacing C 1 → 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 A ππ V V in Eq. (6.23), induced by the non-derivative operator in Eq. (6.11) with coupling constant Z. The amplitude A ππ V V has the same form of A χ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 π − into a π + , see Fig. 9. This implies that in the computation of A A,B,C in Eqs. (3.8)-(3.15) one needs to make the replacement: As discussed in Ref. [43], the UV behavior of this potential induces additional divergences reabsorbed by the coupling C 2 . Finally, we note that in complete analogy to the physical LL case, the matching condition for the V V amplitude will involve only the real part of the singular component of the amplitude A C associated with Eq. (6.23). Combining the N N and ππ contributions, we write it in the following way that will prove useful in the matching procedure:

A V V 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: One can again organize the full theory calculation of A < V V in close analogy to the EFT expression given in Eq.
In the full theory analysis, the counterpart of the EFT terms proportional to Z (called A ππ V V in Eq. (6.23)) is due to the insertion ofÔ V V (k) ≡ (Ô LL +Ô RR + 2Ô LR )(k) on a pion exchanged between the two nucleons, as depicted in Fig. 10. The corresponding A ππ,<,sing C is given by The integral in d 3 q is actually convergent. This can be seen by rewriting T ππ < (k, p) as follows:  Figure 10: Contribution to A V V due to the insertion of two vector currents on a potential pion line. In the small |k| region this diagram represents A ππ,<,sing C . In the large |k| region it representsÂ sing C . See main text for notation.
Note that terms involving q 0 ∼ |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 T ππ < , 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 q; finally evaluating the integral in d|q|. We will present our results in order of complexity: 1. First, we will take F V π (k 2 ) = 1 and I < C (|q|) = θ(|q| − 2|p|)/(8|q|).
2. Next, we will use the realistic pion electromagnetic form factor F V π (k 2 ) from Eq. (6.17).
3. Finally, we will include higher-momentum effects in the NN interactions by writing with r(|q|) defined in Eq. (4.19) and illustrated in Fig. 4.
Throughout, we work in the chiral limit (M π → 0). The result is quite simple when using F V π (k 2 ) = 1 (corresponding to the limit M V → ∞) and I < C (|q|) = θ(|q| − 2|p|)/(8|q|): This expression can be recast into the following useful form using Eq. (6.21) where we have extracted a factor of two for later convenience, have introduced the arbitrary scale ν χ , and isolated the dependence on the IR physics (external momentum |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 ρ propagators need to be added, which changes the bracket in Eq. (6.31) to In the limit ω V → ∞ this reduces to Eq. (6.31), and for |k| → ∞ this expression displays the 1/|k| 4 fall-off that makes the integrand in Eq. (6.31) consistent with the OPE behavior at large |k|. To extract the new form of the integrand, we first express Z < , see Eq. (6.21), as In this way, Eq. (6.32) still applies upon the replacement Finally, we include the HOS form factor effects by rewriting in Eq. (6.28) The first term gives the analytic result in Eqs. (6.32) and (6.35). The second term gives a correction to a ππ < (|k|, ν χ ), which we denote by δa ππ < (|k|), where g(|k|, |q|) is given by The integration in Eq. (6.37) will be done numerically.
Putting together the NN and ππ contributions, the full low-and intermediate-|k| integrand in A <,sing C | V V is given by where the input functions are given in Eqs. (6.27), (6.35), and (6.37). Concerning A > V V , 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 whereÃ C is defined in Section 4.2, g N N LR (ν χ ) andÂ C are defined in Appendix B, and ν χ is the chiral EFT renormalization scale used in evaluating pp|Ô V V |nn . We chose to keep ν χ = µ χ for clarity. However, to simplify the analysis of the IR divergences, we use here the same scale ν χ introduced in Eq. (6.32) to separate out the log |p| term from the low-|k| amplitude A ππ,<,sing C . Given the scaling of the NN couplings and the potentials that determineÃ C andÂ C , discussed in Appendix B, in the above expression we keep only the leading terms g N N LR andÂ C . The real part of the singular component ofÂ C is generated by pion exchange and is given bŷ 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: Using the lattice QCD results of Ref. [90], we findḡ ππ LR = 8.23 in the MS scheme at µ = 2 GeV. While bothḡ ππ LR andḡ N N LR are expected to be O(1), given the large numerical value ofḡ ππ LR , for the nucleon coupling we will assumeḡ N N LR ∈ [−10, +10].
14 Here and below the dependence on the QCD short-distance renormalization scale µ in αs, g ππ LR , and g N N LR is suppressed for simplicity. Note that, in contrast to g N N 1 , g N N LR carries an additional dependence on the chiral EFT scale νχ.

Matching
The matching condition for the vector amplitude A V V leads to (compare to Eq. (5.1) for the LL amplitude) The input needed in Eq. (6.44) can be found in Eqs. (6.25), (6.26), (6.27), (6.32), (6.35), and (6.43). Rescaling the couplings and integrals as in Section 5 and using Z = Z < + Z > , one arrives at with a V V < (|k|, ν χ ) given in Eq. (6.39). Note that the logarithmic IR dependence on the external momentum |p| has disappeared in the matching relation, providing a strong consistency check on the calculation. Except for δa ππ < (|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 Without loss of generality, in the numerical results we set ν χ = M π . The various components of the integrand are shown in Figs. 11 and 12. Fig. 11 displays the chiral, low-and high-|k| components of the integrand, namely a V V χ (|k|) (yellow), a V V < (|k|, M π ) (blue), and a V V > (|k|, M π ), in both linear scale (left panel) and logarithmic scale (right panel). Fig. 12 shows a ππ < (|k|) (blue), δa ππ < (|k|) (yellow), and their sum in green. Fig. 13 shows the dependence of (C 1 +C 2 ) on the matching scale Λ (left panel) and the chiral renormalization scale µ χ (right panel).
We now discuss the uncertainties associated with our estimate of (C 1 +C 2 ): (i) Concerning the parametric uncertainties, the dependence on the vector form factor parameterization is mild: we find an upward shift of ±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 Λ ∈  Figure 12: a ππ < (|k|, M π ) (blue), δa ππ < (|k|) (yellow), and their sum (green).
of the corresponding two-pion matrix element, see Appendix B), is at the level of ±0.3, larger than in the case ofC 1 , see Fig. 13 (left panel). (ii) To estimate the systematic uncertainty due to the choice of the short-range NN interaction V S , we repeat the analysis using the Reid and AV18 potentials, finding shifts in (C 1 +C 2 ) of ≈ +0.25 and ≈ −0.3, respectively. (iii) For the systematic effect due to the neglected inelastic channels we take a range δ(C 1 +C 2 ) ≈ ±1.1, which leads to a relative error of about 50% in the singular NN electromagnetic amplitude at |p| ∼ 25 MeV, larger than the 30% in the weak amplitude controlled byC 1 . In addition to the new class of pion-exchange diagrams, this larger assignment is motivated as follows: in contrast to theC 1 analysis, the parametric error is now more sizable, mainly driven by the NN shortdistance couplingḡ N N LR . In fact, reducing Λ to values as low as 1 GeV and thus into the energy region where the applicability of the OPE becomes questionable and inelastic effects important, leads to a variation δ(C 1 +C 2 ) ±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 Λ 2 GeV, and thus increase the estimate of inelastic contributions accordingly to account for the more prominent Figure 13: Left panel: dependence of (C 1 +C 2 )(µ χ = M π ) on the matching scale Λ usinḡ g N N LR = ±10. Right panel: dependence ofC 1 +C 2 on the chiral renormalization scale µ χ .

Charge-independence-breaking contribution to NN scattering
The result (6.47) already compares quite well to the phenomenological determination (C 1 + C 2 )(µ χ = M π ) = 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 NN 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 C pp (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  [103,104], a nn = −18.9(4) fm [105]. From Eq. (6.47) we find a CIB = 15.9 +4.5 −4.0 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 NN scattering thus validates our approach at the level of (30-50)% and shows that our uncertainty estimates are realistic.

Synthetic data for nn → pp near threshold
Having determined the LECC 1 in the MS scheme, we can now compute the low-energy nn → ppe − e − amplitude A ν . 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 A ν , 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 where both initial |Ψ nn (|p|) and final |Ψ pp (|p |) states are in the 1 S 0 channel, V 1 S 0 ν L denotes the usual long-range neutrino potential, and V 1 S 0 ν S is the short-range interaction proportional tõ C 1 . We denote by E = p 2 /m n and E = p 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 p and p 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 |p|, via 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 |p | ranging between zero and the maximum value given in Eq. (7.2). In Ref. [59] we chose the momenta |p| = 25 MeV and |p | = 30 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 space 15 with the long-range neutrino potential given by 1 − e −Mπr 1 + M π r 2 , (7.4) 15 The relation between C1 andC1 is given in Eq. (5.2).
with C tuned to reproduce either the np or nn scattering length (the difference being negligible at the chosen kinematic point). To solve for the wave functions we have used the methods described in Ref. [23], and to compute the amplitude, isolating its singular component, we have used the approach described in detail in Ref. [43]. Our amplitude satisfies Watson's theorem [106] and at the kinematic point |p| = 25 MeV, |p | = 30 MeV we find 16 where in the second line we have used our predictionC 1 (4M π ) = 4.2 (6). At the chosen kinematic point and chiral renormalization scale we find that the contact-term contribution interferes destructively with the long-range neutrino exchange and reduces the amplitude by about 15%. Since, as discussed in Refs. [40,43], the effect of a contact term of natural size that affects ∆I = 0 transitions such as nn → pp at the (10-20)% level, becomes amplified to the level of (50-70)% in ∆I = 2 nuclear transitions due to a node in the matrix element density, it is possible that a more pronounced effect occurs in nuclei of interest for 0νββ searches, but we stress that our observation is based on the MS scheme at µ = 4M π and thus need not carry over to other scales and schemes. We advocate using the synthetic data (7.6) to determine the LECC 1 in any regularization/renormalization scheme applicable in many-body nuclear calculations [26-28, 107, 108]. As an example, we show here how to perform such a determination for often-applied momentum- 16 The amplitude Aν is related to the S-matrix element for the process n(p) n(−p) → p(p ) p(−p ) e(pe) e(−pe) by and coordinate-space regulators. In momentum space, the strong potential is regulated through an exponential cutoff of the form in terms of a momentum cutoff Λ p . The LS equation is solved numerically for different values of Λ p and n and the strong counter term, C(Λ p , n) in the 1 S 0 channel, is determined from a fit to the NN scattering lengths. We then insert the long-and short-range neutrino potential between initial-and final-state scattering states and determineC 1 (Λ p , n) by fitting to Eq. (7.6). An analogous procedure is followed in coordinate space. In this case, contact interactions are regulated by a local Gaussian regulator The pp and nn wave functions obtained for a fixed value of R S are then used to compute the amplitude (7.3).
In Fig. 14 we depict the values of the counter termC 1 in the two schemes. In the left panel, we showC 1 (Λ p , n) for a wide range of Λ p and three choices of n, showing a clear logarithmic dependence of Λ p and a mild dependence on n. In the right panel, we show the same plot in the R S scheme, which again exhibits a logarithmic dependence on R S , with power corrections becoming important at small R −1 S . The relative importance of the short-distance contribution depends on the regulator and on the kinematics. In Fig. 15 we show the amplitude A ν for |p| = 25 MeV, while varying |p | between 26 and 46 MeV. The region includes the fitting point |p | = 30 MeV, and at these small momenta the LO chiral wave functions provide a good description of the 1 S 0 phase shift. The plots in the top panel use the momentum-space regulator Λ p , while the bottom plots use the coordinate-space R S scheme. In the top-left panel, we have chosen two specific regulators Λ p = 2 fm −1 (blue) and Λ p = 20 fm −1 (red) and kept n = 2 as the dependence on n is mild. We show individually the long-distance (dashed), short-distance (dotted), and sum (solid) contributions. The sum does not depend on Λ p , but the counter term goes from constructive (+15%) to destructive (−20%) between the two choices of the regulators. This is shown more clearly in the right panel of Fig. 15, where the ratio of short-to-long-distance contributions is given for three choices of regulators Λ p = {2, 5, 20} fm −1 in blue, green, and red respectively. This confirms explicitly that the separation of long-and short-distance contributions is not physical.
The bottom panel shows the same plots in the R S scheme, using R −1 S = 1.25 fm −1 (green), R −1 S = 2 fm −1 (blue), and R −1 S = 20 fm −1 (red). Also in this case, the interference goes from destructive at large cutoff to constructive at small cutoff, with the short-range amplitude being approximately zero at 2 fm −1 . Of course, the total amplitude agrees between the Λ p and R S regularization schemes as it should. These plots show that, as already emphasized above, we cannot say whether the new short-distance contribution will add constructively or destructively in ab-initio calculations of nuclear transitions. This question can only be answered within a specific regularization scheme and choice for the strong NN potential.

Conclusions
In this work we have presented the details of the new method developed in Ref. [59] to estimate the LO contact-term contribution to the amplitude nn → ppe − e − through the exchange of light Majorana neutrinos. This is currently the missing ingredient in order to construct complete LO nn → pp transition operators in chiral EFT, to be used in ab-initio nuclear-structure calculations of matrix elements in nuclei of experimental interest for 0νββ searches. Our approach to estimate the contact term is based on a representation of the nn → ppe − e − amplitude as the momentum integral of the neutrino propagator (1/k 2 ) times the generalized forward Compton scattering amplitude n(p 1 )n(p 2 )W + (k) → p(p 1 )p(p 2 )W − (k), in close analogy to the Cottingham formula for the electromagnetic mass splittings of pions and nucleons. To extract the contact term from this integral representation, we have constructed model-independent descriptions of the integrand in the low-momentum region, using chiral EFT, and in the highmomentum region, using the OPE. In the low-and intermediate-momentum region we have kept only the elastic contribution, i.e., the effect of two-nucleon intermediate states, with the most important momentum dependence generated by single-nucleon form factors as well as the NN amplitude itself. While we do not have a strict dispersive derivation as in the case of hadron mass splittings, this approach has the potential to match the (20-30)% accuracy of the elastic approximation observed there, an expectation that we verified by studying CIB in low-energy NN scattering within the same framework and reproducing, within uncertainties, the CIB contribution to the 1 S 0 NN scattering lengths. This phenomenological success gives us confidence that our method is sound and the uncertainty estimate realistic.
The phenomenological validation is particularly important given that our method does introduce model-dependent input in the intermediate-momentum region, albeit anchored to known constraints from QCD at low and high momenta. The extraction of the LO contact term, in a given scheme, then proceeds by matching the full amplitude obtained in this way to the LO chiral EFT amplitude. We considered several sources of uncertainty, chief among them the missing contributions from inelastic intermediate states, such as NN π. The intermediate steps of our analysis have been performed in dimensional regularization with the MS scheme, while the final result can be expressed in terms of the scheme-independent renormalized amplitude A ν (|p|, |p |) at a set of kinematic points near threshold, where LO chiral EFT is expected to give an excellent approximation. Using our synthetic data, as discussed in Section 7, one can then determine the contact term in any regularization and renormalization scheme, in particular the ones employed in nuclear-structure calculations for isotopes of experimental interest for 0νββ searches. This application is timely in view of the remarkable progress in ab-initio calculations of 0νββ decay rates of light and intermediate-mass nuclei [29-31, 43, 109], ranging from 6 He to 48 Ca and 76 Ge, starting from microscopic nuclear forces obtained from chiral EFT. So far these decay rates include the long-distance neutrino-exchange contributions while omitting the contact term, which can now be remedied using our synthetic data for nn → ppe − e − , allowing, for the first time, for complete LO calculations of nuclear 0νββ decay rates. For heavier nuclei such as 136 Xe, which are still beyond the capabilities of ab-initio techniques, the impact of the short-range term could be studied indirectly, e.g., by matching ab-initio results and nuclear models for nuclei accessible to both approaches (this strategy was used in Ref. [110] for the axial-vector current).
As a benchmark point, we find for the amplitude A ν in the MS scheme at µ χ = 4M π , with initial-and final-state nucleon momenta |p| = 25 MeV and |p | = 30 MeV, that the contactterm contribution adds destructively to the neutrino exchange at the 15% level. However, as discussed in detail in Section 7, this statement depends on the renormalization scheme and scale, e.g., the two cutoff schemes studied there can lead to both constructive and destructive effects depending on the choice of scale. This illustrates that the separation into short-and long-distance contributions is unphysical, see Fig. 15 for the decomposition in the Λ p and R S schemes as a function of the final-state momentum and the cutoff scales. Further, as discussed in Refs. [40,43], while a contact term of natural size affects ∆I = 0 transitions such as nn → ppe − e − at the (10-20)% level, its effect is amplified to the level of (50-70)% in ∆I = 2 nuclear transitions due to a node in the matrix element density. The size of the effect in realistic 0νββ nuclear transitions can now be addressed, greatly reducing a crucial uncertainty in the interpretation of future experimental searches [111][112][113][114][115][116][117].
Improving the accuracy of the results presented here would require at least a thorough study of the inelastic NN π channel, but at the same time a more systematic connection to the Cottingham formula would likely become necessary as well: while the identification of the leading elastic effects is rather intuitive, the extension to subleading corrections is not. Additional experimental input could, in principle, separate C 1 from C 2 , e.g., via CIB in nuclei, but such an extraction also requires the development of a suitable theoretical framework. While our results allow for first phenomenological estimates of the impact of the contact term on 0νββ decay rates, they thus also define a benchmark for future lattice-QCD calculations [46,[52][53][54][55][56][57]. In addition to comparing the final result for A ν , there could also be aspects of the matching strategy and spectral representation, as described in Section 2, that might prove synergistic between the two approaches. Finally, having concentrated on the LO contact term for light Majorana exchange in this paper, we remark that contact-term contributions arise at LO for other operators mediating 0νββ decay as well, both at dimension 7 and dimension 9 [37], or through exchange of massive sterile neutrinos [44], and generalizations of the strategies presented here could help constrain the contact terms in these mechanisms.

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

A.1 The half-off-shell form factor in pionless EFT
In this appendix we discuss the HOS form factors in the framework of / πEFT. Since analytic results are readily available, the / πEFT analysis will provide several insights on the general problem at hand.

A.1.1 Half-off-shell T matrix in pionless EFT at NLO
Working in dimensional regularization with power divergence subtraction [25], the resummation of bubble diagrams with NN vertices involving C 0 and C 2 leads to (E = p 2 /m N , q 2 = p 2 ) where we used the definitions (3.3) withV π → 0 and This leads to As mentioned before, these HOS form factors are different depending on whether they involve the initial-or final-state on-shell momenta, cf. Eq. (4.20). By themselves these form factors are clearly not physical. First, one notes that they depend on the renormalization scale µ, as C 2 /C 0 is µ-dependent. The form factors vary significantly with µ, while having the expected qualitative features at small |k|. Second, as shown below, f S can be changed by performing field redefinitions, as expected for any off-shell quantity. However, in the actual application insertions of f S combine into a physical contribution to the nnW + (k) → ppW − (k) amplitude, or equivalently to the integrand in the sum of M < B+B and M < C in Eq. (4.2). In the end, the net effect is a modification of the function I < C (|k|) and hence a change in r(|k|).
A.1.2 Behavior of f S and I < C (|k|) under field redefinitions We consider the non-linear field redefinition given in Ref. [79], namely where η is an arbitrary parameter of mass dimension −3. This transformation induces a derivative NN interaction from the kinetic term as well as a two-nucleon current operator, schematically The resulting η-dependent four-nucleon vertex is [79] (A.6) This modifies the HOS form factors (A.3) as follows These η-dependent HOS form factors combine with new contributions to the current matrix element in such a way to leave no η-dependent term in the amplitude nnW + (k) → ppW − (k) itself. It is simple to verify that cancellations happen through groups of diagrams represented in Fig. 16. Therefore, the integrands in the amplitudes M < B+B and M < C in Eq. (4.2) (and hence I < C (|k|) via Eqs. (4.7) and (4.18)) are individually invariant under the field redefinition considered here, i.e., they do not depend on η. We have used this simple field redefinition of Ref. [79] to illustrate a more generally valid point.
For the real part of I < C (|k|), relevant for the matching, this implies In the last step we have used the / πEFT expressions of C 0,2 (µ) and taken µ ∼ M π 1/a. The size of the effective range r 0 = 2.73 fm = 1/(72 MeV) implies sizable deviations from the leading 1/(8|k|) behavior already at |k| ≤ M π .
While the analysis above is promising, there are two unsatisfactory aspects: (i) Even though it does not enter the matching, the imaginary part of I < C (|k|) depends on the scale µ through the ratio C 2 /C 0 . (ii) The real part of I < C (|k|) is only approximately independent of µ, as we have neglected 1/(µa) 1 in the expression (A.10) Both issues can be addressed by studying how f S impacts all components of the generalized Compton amplitude, i.e., the integrand of M < B+B and M < C in Eq. (4.2), and not just the integrand of M < C as done above. The final result is very simple: the contributions combine to give a real-valued and scale-independent shift to I < C (|k|) equal to δI < C (|k|) = −r 0 /(8π), which agrees with the result in the second line of Eq. (A.9).
The argument proceeds as follows. Omitting a common factor proportional to the weak current insertions, namely −(1 + 3g 2 A )/k 2 , the integrand corresponding to M < B+B + M < C in Eq. (4.2) has the following structure: where T (E) denotes the LO on-shell NN T matrix in the 1 S 0 channel. Substituting the form factors (A.3), one obtains where I 0 (E) is given in Eq. (A.2) and T (E) = C 0 /(1 − C 0 I 0 (E)). Using the relation 13) one sees that the terms proportional to C 2 /C 0 in Eq. (A.12) cancel and one is left with the scale-independent quantity T (E ) (−C 2 /C 2 0 ) T (E). Upon matching the overall factors of m N , this leads to To conclude, we discuss several lessons that can be abstracted from the / πEFT analysis: first, the large fractional correction in I < C (|k|) (and hence a < (|k|)) at relatively low |k| ∼ M π (in / πEFT) is due to the large effective range in the 1 S 0 channel. Second, in / πEFT, the inclusion of HOS form factor effects from M < B+B was needed to cancel a subleading O(1/(µa)) scale dependence in the integrand corresponding to M < C , proportional to I < C (|k|). In practice one gets the leading corrections by just including the HOS form factors in the M < C integrand. This is consistent with the previous result that M < B+B drops out of the LO matching formula and the expectation that the HOS form factor makes the integrand in M < B+B even more convergent. Accordingly, this suggests that the inclusion of the HOS form factors in the integrand of M < C captures the bulk of the physical suppression of the nnW + (k) → ppW − (k) amplitude at small |k|, while the inclusion of HOS form factors in M < B+B would lead to effects that are below the current level of accuracy.

A.2 Half-off-shell effects in chiral EFT at NLO
In this appendix we extend the analysis to chiral EFT, thus keeping pion-exchange effects. To take into account off-shell effects inT S , we start from Eq. (3.3), where q|V S |p = C + C 2 is the NLO short-range potential, so that we can writê (A. 16) HereV (C) is the LO part of the short-distance potential and the last term gives rise to off-shell effects proportional to C 2 . ExpandingT χ S to first order in C 2 we have E can be absorbed into an NLO correction to C → C + C (1) , after which K (0,1) E become scale independent and the NLO on-shell matrix E . Inserting the above off-shell matrix elements q|T χ S |k into Eq. (4.2) we can write the lowmomentum component of the M C amplitude as . The first term in curly brackets contributes to the most singular part of M C (the parts that diverge when integrating over k). In particular, the terms with χ + q,q+k (0) → 1 give rise to the part of I C in Eq. (4.18), with f S → 1. In other words, these pieces are responsible for the contributions involving the on-shell T S matrix. Instead, the contributions from the terms proportional to 1 − χ + q,q+k (0) correspond to diagrams involving pion exchanges within the bubble, which are convergent.
The last line in Eq. (A.18) arises from the off-shell piece ofT S , the second term in Eq. (A.17), and also contributes to the most divergent parts of M C . As in pionless EFT, these terms are µ dependent by themselves and have to be combined with pieces from the M B diagrams in order to obtain a µ-independent answer. This part of the amplitude can be written as where the terms within square brackets are finite, while the last line comes from the off-shell part ofT χ S , which is again µ dependent. Combining the terms that have explicit C 2 dependence with those in M C , we obtain which is now µ independent and has the factors of χ + p (0)K E we expect for the M C -type diagrams. Using the terms in Eq. (A.20) provide a shift to I < C equal to δI < C = −C 2 /(m N C 2 ), which is the same as in / πEFT, up to the fact that C 2 /C 2 cannot simply be expressed in terms of the effective range parameter r 0 .

A.3 Half-off-shell effects in NN potential models
Finally, in this appendix we describe the calculation of the HOS form factor and the ratio r(|k|) in potential models. As a benchmark resonance model for the 1 S 0 channel interaction we will use the three-Yukawa potential, known as the Reid potential [88], which has also been studied in the context of analyzing the convergence of the NN EFT [87]. This is a simplified version of the general resonance saturation model for NN interactions [118]. In coordinate space the potential takes the form with the short-range term modeled through the attractive σ-meson and repulsive ρ-meson contributions. Once the mass parameters are fixed (M σ = 4M π , M ρ = 7M π in Ref. [88] and M σ = 500 MeV, M ρ = 770 MeV in Ref. [87]), the couplings α σ,ρ are tuned to reproduce the 1 S 0 phase shift or the scattering length and effective range [87], producing an overall good representation of the data. In this class of models it is relatively simple to compute the HOS T -matrix elements defined in Eq. We have also performed the calculation of T (q, p) and f S (q, p) with the AV18 potential [89], after subtracting the one-pion exchange. This is a representative of the class of high-quality potentials that reproduce the NN scattering data for |p| out to several hundred MeV.
Representative results from this analysis are summarized in Fig. 17. In the left panel we show f S (q, p) versus |q| at |p| = 1 MeV for the three potentials considered here, Kaplan-Steele, Reid, and AV18. In the right panel we report the ratio r(|k|), see Eq. (4.19), relevant for the calculation of A <,sing C . The difference in r(|k|) provides some indication of the level of model dependence in our matching calculation. While locally the differences can be large, the fall-off of the A <,sing C integrand multiplying r(|k|) dilutes the impact of any feature in r(|k|) above 500 MeV, so that ultimately the impact on the extracted LEC is relatively minor. We expect this conclusion to hold in general, beyond the limited set of potential models explored here.

B Left-right correlator
In the main text we defined the "left-left" correlator, which naturally appears when dealing with the weak currents. In the case of electromagnetic currents and for matching purposes, we also need the "left-right" correlator. We definê With these definitions, the time-ordered products at large k 2 can be written as, see Eq. (4.24), Similarly, we obtain: where the respective last step holds under symmetric integration. Finally, using Fierz identities we can express (Ô LL ) α α and (Ô LR ) α α in terms of the four-quark scalar operator basis used in Ref. [42], obtaining As discussed in Ref. [42], the four-quark local operators O 1,4.5 can be treated at low energy in chiral EFT. The discussion of O 1 is given in Section 4.2 and we focus here on O 4,5 . At LO only non-derivative four-nucleon and pion-pion vertices appear, The ππ couplings induce a pion-range transition operator, while the NN couplings give a shortrange transition operator. The leading chiral expression for the nn → pp matrix element ofŌ LR is given by where [42] g N N LR = The scaling of the pion couplings follows from naive dimensional analysis and is confirmed by a lattice QCD calculation [90]. On the other hand, the scaling of the NN couplings follows from the RGE in chiral EFT [42], which we report here for completeness in terms of the rescaled couplingsg N N 4,5 defined by g N N 4,5 = (m N C/(4π)) 2gN N 4,5 : The scaling for the nucleon couplings follows from the above RGE and the scaling m N C/(4π) ∼ 1/Q. As a consequence of the above discussion, the coupling g N N LR depends on both the shortdistance QCD renormalization scale µ and on the chiral renormalization scale ν χ .
Using the above relations we arrive at For matching purposes only A > C matters. In the numerical estimates we will retain only the term proportional to g N N LR and the singular part ofÂ C , proportional to g ππ LR . Further, for convenience in the matching analysis we express the couplings g N N LR and g ππ LR in terms of dimensionless quantities as follows In numerical estimates we will useḡ ππ LR = 8.23 (in the MS scheme at µ = 2 GeV, based on the lattice QCD results of Ref. [90]) and assume the rangeḡ N N LR ∈ [−10, +10].

C Estimating inelastic effects
As a representative of inelastic effects, we estimate the contribution from the diagram in Fig. 8. 17 The contribution to the generalized Compton amplitude in Eq. (2.9) is already a two-loop diagram, which with appropriate momentum routing can be written as (C.1) Upon changing variables to q and l = q − q one has Recalling that I C → 1/(8|l|) for p = p → 0, one obtains The integral is logarithmically divergent and we regulate the divergence with a cutoff scale ν ∼ Λ χ . We will mimic the effect of the associated (and unknown) counter term by varying the scale ν in order to obtain a rough estimate of the corresponding inelastic effect. The inelastic contribution to A < C is given by We reduce this integral by first integrating over k 0 via the residues at the pion and neutrino poles, and second performing the angular integral in d 3 l, obtaining The integral in d|l| can be done analytically, leading to In the M π → 0 limit the integrand becomes The logarithmic divergence signals that there must be, at this order, a counter term for the generalized Compton amplitude (the divergence arose at that level, before integrating over d 3 k). Naive dimensional analysis suggests that this counter term should be of the same order of the coefficient of log ν. To obtain a numerical estimate we vary the scale within ν ∼ (0.5-1.0) GeV. The new contribution to the integrand is clearly suppressed by two chiral orders compared to the leading terms, as long as |k| Λ χ , as one is comparing the LO 1/|k| to |k|/(4πF π ) 2 . However, this leads to a badly divergent integral in |k|. In order to estimate the contribution to the integral, and therefore to the LEC, we explore two options: first, one may use the chiral EFT form of the integrand and cut it off at Λ χ ∼ 1 GeV, or, second, one may replace 1 + 3g 2 A → g 2 V (k 2 ) + 3g 2 A (k 2 ) in Eq. (C.6) and integrate up to Λ. Numerically, the second option leads to |δC 1 | = 0.15, while the first option leads to |δC 1 | = 0.35 (the largest absolute value is obtained for ν = 0.5 GeV). Based on these observations, we will take into account the effects of inelastic channels by adding a ±0.5 uncertainty to our estimate ofC 1 (µ χ ).

E.1 Low-energy NN scattering in the presence of Coulomb interactions
We begin by recalling some results on the Coulomb-modified effective range expansion for pp scattering [97][98][99]. We specialize to the S-wave and indicate the pure Coulomb phase shift by σ and the full phase shift by σ + ν. Further, we denote the total scattering amplitude by T , the purely Coulomb component by T C , and the Coulomb-modified strong amplitude by T SC T = T C + T SC = 4π m p e 2i(σ(k)+ν(k)) − 1 2ik , . (E.1) The phase shift ν(k) obeys a modified effective range expansion [97][98][99] 1 k (cot(ν(k)) − i) The standard effective range expansion is recovered formally by setting α → 0, which also implies C 2 0 (η) → 1. For α = 0, the function h(η) goes smoothly to zero for k → 0, so that a C can be interpreted as a scattering length. However, the Sommerfeld factor C 2 0 (η)-the square of the Coulomb wave function at the origin-goes to zero for k → 0 due to Coulomb repulsion and this prevents the usual threshold analysis of the amplitude. It will prove useful to work with the momentum-dependent quantityã pp (k) = C 2 0 (η) 1 a C + αm p h(η) , which is a proxy for the very low-energy Coulomb-subtracted pp scattering amplitude.

E.2 Low-energy NN amplitudes in chiral EFT with isospin breaking
Using the techniques of Ref. [23], we extract the phase shift ν(k) from the asymptotic behavior of the solution to the Schrödinger equation with inclusion of Coulomb potential and the LO chiral potential. The 1 S 0 chiral potential includes the LO channel-dependent contact interaction (denoted by C nn,pp,np in the following), as well as the channel-dependent Yukawa potential with the electromagnetically induced pion mass splitting. By including these effects in the Schrödinger equation we effectively iterate the isospin-breaking corrections to all orders. This is mandatory for the Coulomb potential, while for the other terms a perturbative treatment should be numerically close to our analysis, see Ref. [43] for additional details.
The phase shift ν(k) is related to the complex coefficients y pp (k) and z pp (k) controlling the large-r behavior of the regular and irregular solutions of the Schrödinger equation with Coulomb and Yukawa potentials ψ reg (r) → y pp (k) e i(kr−η log(2kr)+σ(k)) kr + c.c. , ψ irr (r) → z pp (k) e i(kr−η log(2kr)+σ(k)) kr + c.c. (E.8) By matching the low-energy form of the amplitude obtained in EFT to the (modified) effective range form, one obtains a np = a π np + G(C np , y np , z np ) , a nn = a π nn + G(C nn , y nn , z nn ) , a pp (k) = a π pp + G(C pp , y pp , z pp ) , G(C ij , y ij , z ij ) = m ij 16πy * 2 (E.9) In the above relations a π ij is the Yukawa-induced scattering length (obtained in terms of the phase of y ij (k) at small k). m ij is the appropriate channel-dependent mass for the nucleon system (m nn = m n , m pp = m p , m np = 2m n m p /(m n + m p )). The quantities z ij depend on the regulator λ introduced when imposing boundary conditions for the irregular solution of the Schrödinger equation [23]. The results quoted below are obtained with 1/λ = 0.001 fm and we have checked stability in the range 1/λ = 0.001 → 0.05 fm. The dependence on λ is canceled by the matching factors ∆ ij (µ/λ), to connect to the MS scheme, in which the contact couplings C ij (µ) are defined. The explicit form of these scheme-changing factors is In the presence of isospin breaking, the short-range 1 S 0 NN couplings can be decomposed as follows where the CIB combination C 1 + C 2 arises form I = 2 interactions, while the CSB term C CSB stems from I = 1 interactions. The latter term can originate from strong isospin-breaking ∝ (m u − m d ) or electromagnetic interactions. For the CIB coupling, renormalization arguments (cancellation of divergences) enforces the LO scaling ∼ e 2 /Q 2 [40,43,100]. For the CSB coupling, the Coulomb potential also induces a UV divergence that requires it to scale as ∼ e 2 /Q 2 .

E.3 Fitting the couplings
Equation (E.9) can be used to extract the three contact couplings C np,nn,pp (or equivalently C, C 1 + C 2 , C CSB ). The output is very stable when using momenta near threshold, in the range k ∈ [1 MeV, 10 MeV]. For definiteness, we quote results for k = 1 MeV. Using the reference scale µ = M π ≡ (M π 0 + 2M π ± )/3, this analysis leads to The result indicates that both CIB and CSB effects provide corrections to the isosymmetric coupling C at the (10-15)% level.

E.4 Predicting the CIB combination of scattering lengths
Alternatively, we can use Eq. (E.9) to validate our calculation of (C 1 + C 2 )(µ) as follows: 1. We use one linear combination of Eq. (E.9) to extract the isospin-conserving coupling C.
2. We construct the CIB combinationsã CIB (k) and a CIB . The theoretical expressions constructed from the right-hand side of Eq. (E.9) depend in general on C, (C 1 + C 2 ), and C CSB . We will fix C to the fit value (see previous step) and consider the dependence of the theoretical expressionsã th CIB (k) and a th CIB on the I = 2 coupling (C 1 + C 2 )(µ). Using our range for (C 1 + C 2 )(µ) we can predictã CIB (k) and a CIB , and compare to the experimental value.
3. Working to first order in isospin-breaking quantities, the observableã th CIB (k) should be insensitive to the I = 1 coupling C CSB . However, in Eq. (E.9) we are effectively iterating the insertions of electromagnetic isospin-breaking sources to all orders (including C CSB ), since we are iterating the full couplings C ij . As indicated by C CSB /C ∼ 15%, second-order effects due to two insertions of C CSB might not be negligible. To quantify the impact of this on our CIB analysis, in Fig. 19 we show for k 0 = 1 MeV, µ = M π (left panel), and µ = 2M π (right panel) the two functions Not including (C 1 +C 2 ) would lead to a scale-dependent prediction forã CIB (k). Moreover, at µ = M π the theory would predictã CIB (k = 1 MeV) ∼ 30 fm, more than a factor of two larger than the experimental value. (ii) The two choices of treating C CSB lead to a shift ∆ã CIB ∼ 2.5 fm at µ = M π (growing to ∼ 4 fm at µ = 2M π ), subdominant but not negligible. (iii) One could use the intercept of the curves in Fig. 19 with the horizontal experimental constraint to perform an alternative extraction of (C 1 +C 2 ). Due to the missing CSB contribution, usingã th (1) CIB implies a value of (C 1 +C 2 )(M π ) ∼ 3.7, to be contrasted with the 5.1 from the previous analysis. (iv) Finally, we note that usingã th(2) CIB produces results that are more stable under variation of µ (see below). In summary, we have derived a prediction of the CIB scattering length that overshoots the experimental value by about 35% and has a comparable uncertainty, thus comparing well with data. This is a significant phenomenological success of our theoretical approach and supports the validity of our uncertainty estimates.  (2) CIB (upper curve) versus the dimensionless coupling C 1 +C 2 . The curves are evaluated at µ = M π (left panel) and µ = 2M π (right panel). The horizontal line represents the experimental value a CIB = 10.35 fm. For reference, our theoretical ranges for the CIB LEC are (C 1 +C 2 )(M π ) = 2.9(1.2) and (C 1 +C 2 )(2M π ) = 5.4(1.2).