Neutrinoless double-beta decay in a finite volume from relativistic effective field theory

The pionless EFT is based on the tenet that the few-nucleon processes at very low energies, i.e., $Q\ll m_{\pi}$, are not sensitive to the details associated with pions or other meson exchanges [24]. Then, all mesons can be integrated out and the effective Lagrangian contains nucleon and lepton degrees of freedom, organized according to the number of derivatives. The observables are expanded in $Q/m_{\pi}$, where $Q$ is the low-energy scale of the order of the binding momentum $\gamma=\sqrt{m_{N}B_{NN}}$ or of the inverse of the scattering length $a$. At LO, the effective Lagrangian reads

$$\begin{split}\mathcal{L}_{\Delta L=0}&=\overline{\Psi}({\rm i}\gamma^{\mu}\partial_{\mu}-m_{N})\Psi-\sum_{\alpha}\frac{C_{\alpha}}{2}(\overline{\Psi}\Gamma_{\alpha}\Psi)^{2}+\frac{1}{2}\overline{\Psi}(l_{\mu}\gamma^{\mu}+g_{A}l_{\mu}\gamma^{\mu}\gamma_{5})\Psi,\end{split}$$

(1)

with $\Psi$ the nucleon field , $m_{N}$ the nucleon mass, LECs $C_{\alpha}$, $\alpha=S,V,AV,T$, and $\Gamma_{\alpha}$ the corresponding Dirac gamma matrices. The nucleons are coupled to the electroweak current $l_{\mu}$ via both vector coupling $g_{V}=1$ and axial coupling $g_{A}=1.27$. Here, we neglect the dependence of $g_{A}$ on the unphysical pion mass. This dependence can be easily taken into account once its value is provided by the LQCD calculations at the corresponding pion mass since it only provides constant factors on the neutrino potential. The electroweak current reads $l_{\mu}=-2\sqrt{2}G_{F}V_{ud}\tau^{+}\overline{e}\gamma_{\mu}v_{eL}+{\rm h.c.}$, with the Fermi constant $G_{F}$ and the $V_{ud}$ element of the Cabibbo-Kobayashi-Maskawa (CKM) matrix [78, 79].

After expanding the nucleon field in terms of the free Dirac spinors $u(\bm{p},s)$ with positive energies in momentum space and keeping only the leading term of $u(\bm{p},s)$ expanding in powers of small momenta $\bm{p}$, the LO strong potential takes the following form,

$$\begin{split}V_{S}(\bm{p}^{\prime},\bm{p})&=\frac{m_{N}}{\omega_{p^{\prime}}}\frac{m_{N}}{\omega_{p}}[C_{S}+C_{V}-(C_{AV}-2C_{T})\bm{\sigma}_{1}\cdot\bm{\sigma}_{2}],\end{split}$$

(2)

where $\bm{p}^{\prime}$ and $\bm{p}$ are the nucleon’s final and initial momenta in the center-of-mass frame, respectively, and $\omega_{p}=(m_{N}^{2}+p^{2})^{1/2}$. At the leading order of pionless EFT, the interaction only contributes to the $s$ wave. Here, we consider the ${}^{1}S_{0}$ channel relevant for the $nn\rightarrow ppee$ process, in which $\bm{\sigma}_{1}\cdot\bm{\sigma}_{2}=-3$, and the LEC $C=C_{S}+C_{V}+3(C_{AV}-2C_{T})$ determines the interaction strength in this channel. As in the nonrelativistic case, the LEC $C$ scales as $\mathcal{O}(4\pi/(m_{N}Q))$ so that the LO amplitudes consist of a resummation of the LO interaction $V_{S}$ and the low-energy pole in the two-nucleon ${}^{1}S_{0}$ amplitude can be reproduced. In particular, for the two-nucleon scattering process, the amplitude reads

$${\rm i}\mathcal{M}_{S}(E)=-{\rm i}V_{S}(p,p)+\int\frac{{\rm d}^{3}k}{(2\pi)^{3}}(-{\rm i})V_{S}(p,k)\frac{{\rm i}}{E-2\omega_{k}+{\rm i}0^{+}}{\rm i}\mathcal{M}_{S}(E).$$

(3)

with $p=\sqrt{\frac{1}{4}E^{2}-m_{N}^{2}}$. After regularizing the contact term with a separable regulator, $V_{S}(p^{\prime},p)\rightarrow V_{S}(p^{\prime},p)f_{\Lambda}(p^{\prime})f_{\Lambda}(p)$, with momentum cutoff $\Lambda$, the scattering amplitude takes the form

$$\mathcal{M}_{S}(E)=-\frac{1}{C_{\Lambda}^{-1}-I_{\Lambda}(E)},$$

(4)

where the “bubble integral” reads

$$I_{\Lambda}(E)=\int\frac{{\rm d}^{3}k}{(2\pi)^{3}}\frac{m_{N}^{2}}{\omega_{k}^{2}}\frac{1}{E-2\omega_{k}+{\rm i}0^{+}}[f_{\Lambda}(k)]^{2}.$$

(5)

The divergence of the “bubble integral” is absorbed into the cutoff dependence of the LEC $C_{\Lambda}$ and the resulting scattering amplitude is cutoff-independent as $\Lambda\rightarrow\infty$. For the two-nucleon bound state, the energy eigenvalue $E$ and the momentum-space wave function $\phi_{E}(\bm{p})$ are obtained by solving the following eigen equation

$$(E-2\omega_{p})\phi(\bm{p})=\int\frac{{\rm d}^{3}k}{(2\pi)^{3}}V_{S}(p,k)\phi(\bm{k}).$$

(6)

Beyond LO, the strong interactions arise from the Lorentz-invariant contact Lagrangian with an increasing number of derivatives (See Ref. [80] for the expressions of the contact Lagrangian up to fourth-order derivatives). Here, we consider the ordering of the contact terms in the relativistic pionless EFT to be the same as its nonrelativistic counterpart, because they both should reproduce the effective range expansion order by order, $p\cot\delta(p)=-\frac{1}{a}+\frac{1}{2}rp^{2}\ldots$, with the ellipse denoting high-order momentum dependences. The subleading contribution is given by the effective range $r\sim O(m_{\pi}^{-1})$, and it is of the order $O(Q/m_{\pi})$ relative to the LO contribution. Once properly renormalized, the difference between the LO relativistic and nonrelativistic pionless EFTs is that the former includes higher-order terms dictated by the Lorentz invariance. For example, the LO relativistic EFT yields a small effective range of $r\sim O(m_{N}^{-1})$, while the nonrelativistic theory yields exactly zero effective range $r=0$.

In this work, we consider the standard mechanism of $0\nu\beta\beta$ decay, in which the lepton number violation at low energies is dominated by a Majorana mass term

$$\mathcal{L}_{\Delta L=2}=-\frac{m_{\beta\beta}}{2}\nu^{T}_{eL}C\nu_{eL},$$

(7)

where $C={\rm i}\gamma_{2}\gamma_{0}$ denotes the charge conjugation matrix and $m_{\beta\beta}$ the effective neutrino mass. Following naive dimension analysis, the LO neutrino potential is contributed only by the long-range exchange of a potential neutrino,

$$V_{\nu{\rm L}}(\bm{p}^{\prime},\bm{p})=(1+3g_{A}^{2})\frac{m_{N}}{\omega_{p^{\prime}}}\frac{m_{N}}{\omega_{p}}\frac{1}{|\bm{p}^{\prime}-\bm{p}|^{2}},$$

(8)

and its contribution is $O(Q^{-2})$. Similar to the strong potential, there are $m_{N}/\omega_{p}$ factors coming from the expansion of the nucleon field without performing the nonrelativistic reduction. However, in pionless EFT, it is known that contact operators in weak processes are typically enhanced when $S$ waves are involved [81, 24]. This leads to the existence of a contact term in the LO neutrino potential–despite being expected two orders down in the naive dimension analysis

$$V_{\nu{\rm CT}}(\bm{p}^{\prime},\bm{p})=-2\frac{m_{N}}{\omega_{p^{\prime}}}\frac{m_{N}}{\omega_{p}}g_{\nu}^{NN}.$$

(9)

with $g_{\nu}^{NN}$ the corresponding LEC. However, different from the nonrelativistic case, $g_{\nu}^{NN}$ does not serve as a counterterm that absorbs the ultraviolet divergence, since the LO long-range amplitude itself is ultraviolet finite [33]. As a result, the contact-term contribution to the LO amplitude is finite in relativistic pionless EFT and can be estimated by integrating out the pion contributions in the relativistic chiral EFT, as we will demonstrate below.

The LO long-range $nn\rightarrow ppee$ amplitude can be schematically written as

$$\begin{split}\mathcal{M}_{0\nu,{\rm L}}&=-m_{\beta\beta}(1+3g_{A}^{2})\left(V_{\nu{\rm L}}-\mathcal{M}_{S}I^{\infty}-\overline{I}^{\infty}\mathcal{M}_{S}+\mathcal{M}_{S}J^{\infty}\mathcal{M}_{S}\right),\end{split}$$

(10)

with $\mathcal{M}_{S}$ is the LO two-nucleon scattering amplitude [Eq. (4)], and $I^{\infty}$ and $J^{\infty}$ are the one- and two-loop integrals with a neutrino exchange. Here, the four terms correspond to the four diagrams depicted in Fig. 1. The first three terms are the contributions in which a neutrino propagates between two external nucleons. The last term, where a neutrino propagates between two nucleons dressed by strong interactions in both the initial and final states, is the subject of matching to LQCD matrix elements [59],

$$\mathcal{M}_{0\nu,{\rm L}}^{({\rm int})}(E_{f},E_{i})=-m_{\beta\beta}(1+3g_{A}^{2})\mathcal{M}_{S}(E_{f})J^{\infty}(E_{f},E_{i})\mathcal{M}_{S}(E_{i}),$$

(11)

where $J^{\infty}$ takes the form,

$$J^{\infty}(E_{f},E_{i})=\int\frac{{\rm d}^{3}k_{1}}{(2\pi)^{3}}\frac{{\rm d}^{3}k_{2}}{(2\pi)^{3}}\frac{m_{N}^{2}}{\omega_{k_{1}}^{2}}\frac{1}{E_{f}-2\omega_{k_{1}}+{\rm i}0^{+}}\frac{m_{N}^{2}}{\omega_{k_{2}}^{2}}\frac{1}{E_{i}-2\omega_{k_{2}}+{\rm i}0^{+}}\frac{1}{|\bm{k}_{1}-\bm{k}_{2}|^{2}}.$$

(12)

This two-loop integral is ultraviolet finite in the present relativistic formulation, but divergent in the nonrelativistic case. This is because, in the nonrelativistic approach, the $1/m_{N}$ expansion is carried out for the integrand making its ultraviolet behavior more singular and resulting in a logarithmic divergence.

Now, we discuss the contact-term contribution to the LO $nn\rightarrow ppee$ amplitude. It is well known that the coupling of the axial current to pions gives rise to the Goldberg-Treiman relation between the pseudoscalar and axial contribution to the weak form factor. In the relativistic chiral EFT, this effect is included in the LO long-range neutrino potential,

$$\delta V_{\nu\rm L}(\bm{p}^{\prime},\bm{p})=-g_{A}^{2}\frac{m_{N}}{\omega_{p^{\prime}}}\frac{m_{N}}{\omega_{p}}\frac{|\bm{p}^{\prime}-\bm{p}|^{2}+2m_{\pi}^{2}}{\left(|\bm{p}^{\prime}-\bm{p}|^{2}+m_{\pi}^{2}\right)^{2}},$$

(13)

while in the relativistic pionless EFT, the contributions from $\delta V_{\nu\rm L}$ are integrated out and will manifest in the form of a contact term. Therefore, the LO contact term $g_{\nu}^{NN}$ can be estimated by the contribution of $\delta V_{\nu\rm L}$. This is achieved by inserting $\delta V_{\nu\rm L}$ into the four diagrams in Fig. 1 and considering $m_{\pi}$ as a hard scale. The first tree-level diagram is just $\delta V_{\nu\rm L}\sim O(m_{\pi}^{-2})$, i.e., next-to-next-to-leading order (NNLO). The contribution of the second diagram can be written as $-g_{A}^{2}\delta I^{\infty}_{\pi}(E_{f},E_{i})\mathcal{M}_{S}(E_{i})$ with $\delta I^{\infty}_{\pi}(E_{f},E_{i})$ the one loop integral with an insertion of $\delta V_{\nu\rm L}$,

$$\delta I^{\infty}(E_{f},E_{i})=\int\frac{{\rm d}^{3}k}{(2\pi)^{3}}\frac{m_{N}^{2}}{\omega_{k}^{2}}\frac{1}{E_{f}-2\omega_{k}+{\rm i}0^{+}}\frac{|\bm{k}-\bm{p}_{i}|^{2}+2m_{\pi}^{2}}{\left(|\bm{k}-\bm{p}_{i}|^{2}+m_{\pi}^{2}\right)^{2}}.$$

(14)

Expanding this integral around the threshold $E_{f}=E_{i}=0$, we have $\delta I^{\infty}(E_{f},E_{i})=\delta I^{\infty}(0,0)+O(Q^{2}/m_{\pi}^{2})$. Then, $I_{\pi}^{\infty}(0,0)$ is just a function of the hard scales $m_{\pi}$ and $M$, and dimension analysis determines $\delta I^{\infty}(E_{f},E_{i})\sim O\left(m_{N}/(4\pi m_{\pi})\right)$. Since the scaling of $\mathcal{M}_{S}$ is $O(4\pi/(m_{N}Q))$, the contribution of second diagram is subleading, $O((m_{\pi}Q)^{-1})$. Following the above analysis, the contribution of the third diagram is also subleading, $O((m_{\pi}Q)^{-1})$, but the contribution of the fourth diagram is instead LO, $O(Q^{-2})$. The latter takes the form

$$\delta\mathcal{M}_{0\nu}^{({\rm int})}(E_{f},E_{i})=m_{\beta\beta}g_{A}^{2}\mathcal{M}_{S}(E_{f})\delta J^{\infty}(E_{f},E_{i})\mathcal{M}_{S}(E_{i})$$

(15)

with

$$\delta J^{\infty}(E_{f},E_{i})=\int\frac{{\rm d}^{3}k_{1}}{(2\pi)^{3}}\frac{{\rm d}^{3}k_{2}}{(2\pi)^{3}}\frac{m_{N}^{2}}{\omega_{k_{1}}^{2}}\frac{1}{E_{f}-2\omega_{k_{1}}+{\rm i}0^{+}}\frac{m_{N}^{2}}{\omega_{k_{2}}^{2}}\frac{1}{E_{i}-2\omega_{k_{2}}+{\rm i}0^{+}}\frac{|\bm{k}_{1}-\bm{k}_{2}|^{2}+2m_{\pi}^{2}}{(|\bm{k}_{1}-\bm{k}_{2}|^{2}+m_{\pi}^{2})^{2}}.$$

(16)

This two-loop integral has a mass dimension of two and thus scales as $O\left(m_{N}^{2}/(4\pi)^{2}\right)$. As a result, $\delta\mathcal{M}_{0\nu}^{({\rm int})}(E_{f},E_{i})\sim O(Q^{-2})$ and needs to accounted for by a LO contact term. The contact term contributes to the amplitude as

$$\mathcal{M}_{0\nu,{\rm CT}}^{({\rm int})}(E_{f},E_{i})=2m_{\beta\beta}\tilde{g}_{\nu}^{NN}\left(\frac{m_{N}}{4\pi}\right)^{2}\mathcal{M}_{S}(E_{f})\mathcal{M}_{S}(E_{i})$$

(17)

after defining the dimensionless LEC $\tilde{g}_{\nu}^{NN}$,

$$\tilde{g}_{\nu}^{NN}=\left(\frac{4\pi}{m_{N}C}\right)^{2}g_{\nu}^{NN}.$$

(18)

By matching Eqs. (15) and (17) at the threshold $E_{f}=E_{i}=0$, the dimensionless LEC can be determined,

$$\tilde{g}_{\nu}^{NN}=\frac{(4\pi)^{2}}{2m_{N}^{2}}g_{A}^{2}\delta J^{\infty}(0,0).$$

(19)

In addition to the contribution originating from the coupling of pions, there could be other unknown short-range contributions to the LO LEC $\tilde{g}_{\nu}^{NN}$ in the relativistic pionless EFT. By comparing the pionless and chiral EFT amplitudes, the unknown short-range contributions in the pionless EFT corresponds to the $\tilde{g}_{\nu}^{NN}$ contact-term contribution in the chiral EFT. It then follows that, based on Weinberg power counting, $\tilde{g}_{\nu}^{NN}$ contribution is expect to be suppressed by two orders in the chiral expansion, i.e., $O(Q^{2}/\Lambda_{\chi}^{2})$ with $\Lambda_{\chi}\sim m_{N}$ the break down scale of chiral EFT. Based on this estimation, the uncertainty of the estimation by Eq. (19) can be considered subleading, since $Q^{2}/m_{N}^{2}\lesssim Q/m_{\pi}$ when $m_{\pi}\leq 806$ MeV (see Table 1).

Finally, adding up the long-range and contact-term contributions, the amplitude is given by

$$\mathcal{M}_{0\nu}^{({\rm int})}(E_{f},E_{i})=-m_{\beta\beta}\mathcal{M_{S}}(E_{f})\left[(1+3g_{A}^{2})J^{\infty}(E_{f},E_{i})-2\tilde{g}_{\nu}^{NN}\left(\frac{m_{N}}{4\pi}\right)^{2}\right]\mathcal{M_{S}}(E_{i}).$$

(20)

In this work, we implement the relativistic pionless EFT in a finite-volume (FV) cubic box with spatial extent $L$ and the periodic boundary conditions. Considering the $nn\rightarrow ppee$ process with the kinematics that the two electrons in the final state are at rest, the Euclidean four-point function accessible from LQCD can be analytically continued to the Minkowski spacetime [59],

$$\mathcal{T}_{L}^{(M)}(E_{f},E_{i})=\int{\rm d}z_{0}\int_{L}{\rm d}^{3}z[\langle E_{f},L|T[\mathcal{J}(z_{0},\bm{z})S_{\nu}(z_{0},\bm{z})\mathcal{J}(0)]|E_{i},L\rangle]_{L},$$

(21)

where the subscript $L$ on the spatial integral indicates that the integral is performed over the FV cubic box and $|E,L\rangle$ is the normalized FV $s$-wave two-nucleon state with the center-of-mass energy $E$. Here, $\mathcal{J}$ denotes the hadronic part of the weak current, and the neutrino propagator $S_{\nu}(z_{0},\bm{z})$ in a finite volume is given by the Fourier transformation

$$S_{\nu}(z_{0},\bm{z})=\frac{1}{L^{3}}\sum_{\begin{subarray}{c}\bm{q}\in\frac{2\pi}{L}\mathbb{Z}^{3}\\ \bm{q}\neq\bm{0}\end{subarray}}\int\frac{{\rm d}q_{0}}{2\pi}{\rm e}^{{\rm i}\bm{q}\cdot\bm{z}-{\rm i}q_{0}z_{0}}\frac{-{\rm i}m_{\beta\beta}}{q_{0}^{2}-\bm{q}^{2}+{\rm i}0^{+}},$$

(22)

ignoring the small nonzero neutrino mass. Because the space is limited to a box with the periodic boundary conditions, the momentum modes are discretized, only taking the values with $2\pi/L$ times three-dimensional Cartesian vectors with integer components. The infrared divergence is regulated by removing the zero-momentum mode of neutrinos.

The energy eigenvalues of the two-nucleon states are also discretized in a finite volume. Their discrete values $E_{n}$ are directly related to the two-nucleon scattering amplitudes in the infinite volume, by the Lüscher quantization condition $\mathcal{F}^{-1}(E_{n})=0$ [82, 83]. For the ${}^{1}S_{0}$ channel considered in this work, it is

$$\mathcal{F}^{-1}(E)=\frac{4\pi}{m_{N}}\left(-\frac{1}{\pi L}\mathcal{Z}_{00}\left[1,\left(\frac{pL}{2\pi}\right)^{2}\right]+{\rm i}p\right)^{-1}+\mathcal{M}_{S}(E),$$

(23)

where $\mathcal{Z}_{00}$ is the zeta function defined in Ref. [82] and $\mathcal{M}_{S}(E)$ is the scattering amplitude defined in Eq. (4).

For the case in which the initial and final states are “scattering” states, the Minkowski matrix element $\mathcal{T}_{L}^{(M)}$ is calculated as following [59],

$$L^{6}|\mathcal{T}_{L}^{(M)}(E_{f},E_{i})|^{2}=|\mathcal{R}(E_{f})||\mathcal{M}^{({\rm int},L)}_{0\nu}(E_{f},E_{i})|^{2}|\mathcal{R}(E_{i})|,$$

(24)

where two FV quantities $\mathcal{M}^{({\rm int},L)}_{0\nu}$ and $\mathcal{R}$ are involved. The former one corresponds to the amplitude (Eq. 20) in a finite volume [59, 60],

$$\mathcal{M}_{0\nu}^{({\rm int},L)}(E_{f},E_{i})=-m_{\beta\beta}\mathcal{M_{S}}(E_{f})\left[(1+3g_{A}^{2})J^{L}(E_{f},E_{i})-2\tilde{g}_{\nu}^{NN}\left(\frac{m_{N}}{4\pi}\right)^{2}\right]\mathcal{M_{S}}(E_{i}).$$

(25)

The function $J^{L}$ resembles the function $J^{\infty}$ in the infinite volume [Eq. (12)], with the momentum integrals replaced by the sum of discrete momentum modes,

$$J^{L}(E_{f},E_{i})=\frac{1}{L^{6}}\sum_{\begin{subarray}{c}\bm{k_{1}},\bm{k_{2}}\in\frac{2\pi}{L}\mathbb{Z}^{3}\\ \bm{k_{1}}\neq\bm{k_{2}}\end{subarray}}\frac{m_{N}^{2}}{\omega_{k_{1}}^{2}}\frac{1}{E_{f}-2\omega_{k_{1}}}\frac{m_{N}^{2}}{\omega_{k_{2}}^{2}}\frac{1}{E_{i}-2\omega_{k_{2}}}\frac{1}{|\bm{k}_{1}-\bm{k}_{2}|^{2}}.$$

(26)

Here, the imaginary part of the propagator is dropped since now the denominator takes nonzero discrete values. The FV quantity $\mathcal{R}(E)$ is the generalized Lellouch-Lüscher residue matrix [61],

$$\mathcal{R}(E_{n})=\lim_{E\rightarrow E_{n}}(E-E_{n})\mathcal{F}(E)=\left(\left.\frac{\rm d\mathcal{F}^{-1}}{{\rm d}E}\right|_{E=E_{n}}\right)^{-1},$$

(27)

which is the residue of the FV function $\mathcal{F}$ (23) at FV energies $E_{n}$.

For the case in which the initial and final states are bound states, $E=2M-B$ with $B>0$, the Minkowski matrix element $\mathcal{T}_{L}^{(M)}$ is calculated by

$$\begin{split}\mathcal{T}_{L}^{(M)}(E_{f},E_{i})&=m_{\beta\beta}\frac{1}{L^{6}}\sum_{\begin{subarray}{c}\bm{k_{1}},\bm{k_{2}}\in\frac{2\pi}{L}\mathbb{Z}^{3}\\ \bm{k_{1}}\neq\bm{k_{2}}\end{subarray}}\phi^{*}_{E_{f},L}(\bm{k}_{1})V_{\nu{\rm L}}(\bm{k}_{1},\bm{k}_{2})\phi_{E_{i},L}(\bm{k}_{2})\\ &-2m_{\beta\beta}\tilde{g}_{\nu}^{NN}\left(\frac{MB}{4\pi}\right)^{2}|\phi(\bm{0})|^{2}.\end{split}$$

(28)

with $\phi_{E,L}(\bm{k})$ the normalized momentum-space wave function of the FV state $|E,L\rangle$ in Eq. (21). On the right hand side, the first and the second terms are respectively the expectations of the long-range neutrino potential $V_{\nu\rm L}$ in Eq. (8) and the contact term in Eq. (9) regulated with the same separable regulator as the one for the strong interaction. The wave function $\phi_{E,L}$ is solved by Eq. (6) with discrete momentum modes. Although there is no two-nucleon bound state in the ${}^{1}S_{0}$ channel at the physical pion mass, the above matrix element is relevant for the study at the unphysical pion masses. At the unphysical pion masses, two nucleons might exhibit a ${}^{1}S_{0}$ bound state predicted by the LQCD calculations [73, 74, 75, 76]. Note that there is an ongoing discussion on whether such a bound state exists at the unphysical pion masses, as several newer works [84, 85, 77, 86] do not identify such a bound state.

The momentum integrals (5), (12), (6) associated with the calculations in the infinite volume are all calculated numerically using Gaussian quadrature. For the calculations of the scattering amplitudes, the real and imaginary parts of the two-nucleon free propagator are calculated separately,

$$\frac{1}{E-2\omega_{k}+{\rm i}0^{+}}=P\left(\frac{1}{E-2\omega_{k}}\right)-{\rm i}\pi\delta(E-2\omega_{k}).$$

(29)

Here, $P$ denotes principle-value integral and it is eliminated by a standard subtraction technique [87]. The eigen equation (6) for the bound states is solved by matrix diagonalization on the Gaussian grids.

However, special care has to be taken for the infrared singularity of the neutrino potential. For the calculation of the scattering amplitude $\mathcal{M}_{0\nu}^{({\rm int})}$, inserting the separation (29) into the expression of the two-loop integral $J^{\infty}$ (12), we have

$$\begin{split}{\rm Re}J^{\infty}(E_{f},E_{i})&=\int_{0}^{\infty}\frac{k_{1}^{2}{\rm d}k_{1}}{2\pi^{2}}\left[\frac{m_{N}^{2}}{\omega_{k_{1}}^{2}}P\left(\frac{1}{E_{f}-2\omega_{k_{2}}}\right)\int_{0}^{\infty}\frac{k_{2}^{2}{\rm d}k_{2}}{2\pi^{2}}\right.\\ &\left.\frac{m_{N}^{2}}{\omega_{k_{2}}^{2}}P\left(\frac{1}{E_{i}-2\omega_{k_{2}}}\right)\frac{1}{4k_{1}k_{2}}\ln\frac{(k_{1}+k_{2})^{2}}{(k_{1}-k_{2})^{2}}\right]-\frac{m_{N}^{2}}{32\pi^{2}}\frac{m_{N}}{\omega_{p_{f}}}\frac{m_{N}}{\omega_{p_{i}}}\ln\frac{p_{f}+p_{i}}{|p_{f}-p_{i}|}.\end{split}$$

(30)

On the right-hand side, there are logarithmic divergences in the two terms when $E_{f}=E_{i}$. We introduce a subtraction technique making use of the analytic expression of the following two-loop integral,

$$\begin{split}I^{\infty}_{\Lambda}&=\int_{0}^{\Lambda}\frac{k_{1}^{2}{\rm d}k_{1}}{2\pi^{2}}P\left(\frac{1}{E_{f}-k_{1}^{2}/m_{N}}\right)\int_{0}^{\infty}\frac{k_{2}^{2}{\rm d}k_{2}}{2\pi^{2}}P\left(\frac{1}{E_{i}-k_{2}^{2}/m_{N}}\right)\frac{1}{4k_{1}k_{2}}\ln\frac{(k_{1}+k_{2})^{2}}{(k_{1}-k_{2})^{2}}\\ &=\frac{m_{N}^{2}}{32\pi^{2}}\ln\frac{\Lambda^{2}-p_{i}^{2}}{|p_{f}^{2}-p_{i}^{2}|}.\end{split}$$

(31)

Denoting the integral term in Eq. (30) as $I^{\infty}$, then we have

$${\rm Re}J^{\infty}(E_{f},E_{i})=\left(I^{\infty}-\frac{m_{N}^{2}}{\omega_{p_{f}}\omega_{p_{i}}}I^{\infty}_{\Lambda}\right)+\frac{m_{N}^{2}}{\omega_{p_{f}}\omega_{p_{i}}}\frac{m_{N}^{2}}{32\pi^{2}}\ln\frac{\Lambda^{2}-p_{i}^{2}}{(p_{f}+p_{i}^{2})}$$

(32)

Now, the two terms are both infrared convergent when $E_{f}=E_{i}$ and we have confirmed the numerical stability using the above expression.

For the calculations of the matrix element $\mathcal{T}_{L}^{(M)}$ between bound states, the infrared singularity of the neutrino potential is treated with the Lande subtraction [88, 89, 90, 91],

$$\int_{0}^{\infty}\frac{{\rm d}k}{2\pi^{2}}\frac{1}{4kp}\ln\frac{(k+p)^{2}}{(k-p)^{2}}k^{2}f(k)=\int_{0}^{\infty}\frac{{\rm d}k}{2\pi^{2}}\frac{1}{4kp}\ln\frac{(k+p)^{2}}{(k-p)^{2}}[k^{2}f(k)-p^{2}f(p)]+\frac{1}{8}pf(p)$$

(33)

with $f(p)$ an arbitrary smooth function.

For the calculations of the matrix elements $\mathcal{T}_{L}^{(M)}$ between the “scattering” states, the FV quantity $J^{L}$ (26), which involves summation over the discrete three-momenta $\bm{k}_{1},\bm{k}_{2}=\bm{n}2\pi/L$ with $\bm{n}\in\mathbb{Z}^{3}$, is calculated using the method of tail-singularity separation (TSS) described in Ref. [92]. In this method, the summation is split into two pieces. One piece contains the singular contributions around $\bm{k}_{1}=\bm{k}_{2}$, but it is exponentially decaying when $|\bm{k}_{1}|,|\bm{k}_{2}|\rightarrow\infty$. The other piece contains a power-law decaying tail at $|\bm{k}_{1}|,|\bm{k}_{2}|\rightarrow\infty$, but it is sufficiently smooth so that it can be approximated by its integral counterpart. Based on this method, we calculate $J^{L}$ as the following

$$\begin{split}J^{L}(E,E)&=\frac{m_{N}^{2}}{4(2\pi)^{6}}\left\{\sum_{\begin{subarray}{c}\bm{n}_{1}\in\mathbb{Z}^{3}\end{subarray}}\frac{\tilde{m}}{\tilde{m}^{2}+n_{1}^{2}}\frac{1}{\tilde{\omega}_{n_{1}}-\tilde{\omega}_{p}}\left[\mathcal{X}_{\rm sum}(\bm{n}_{1},\tilde{p}^{2})+{\rm e}^{-\alpha(n_{1}^{2}-\tilde{p}^{2})}\mathcal{X}_{\rm int}(n_{1}^{2},\tilde{p}^{2})\right]\right.\\ &\left.+4\pi\int_{0}^{\infty}n_{1}^{2}{\rm d}n_{1}\frac{\tilde{m}}{\tilde{m}^{2}+n_{1}^{2}}\frac{1-{\rm e}^{-\alpha(n_{1}^{2}-\tilde{p}^{2})}}{\tilde{\omega}_{n_{1}}-\tilde{\omega}_{p}}\mathcal{X}_{\rm int}(n_{1}^{2},\tilde{p}^{2})\right\}+O({\rm e}^{-\pi^{2}/\alpha}),\end{split}$$

(34)

where $\tilde{m}=2\pi m_{N}/L$, $\tilde{p}=2\pi p/L$, $\tilde{\omega}=(\tilde{m}^{2}+\tilde{p}^{2})^{1/2}$, and

$$\begin{split}\mathcal{X}_{\rm sum}(\bm{n}_{1},\tilde{p}^{2})&=\sum_{\begin{subarray}{c}\bm{n}_{2}\in\mathbb{Z}^{3}\\ \bm{n}_{2}\neq\bm{n}_{1}\end{subarray}}\frac{\tilde{m}}{\tilde{m}^{2}+n_{2}^{2}}\frac{1}{\tilde{\omega}_{n_{2}}-\tilde{\omega}_{p}}\frac{1}{|\bm{n}_{1}-\bm{n}_{2}|^{2}}\left[1-(1-{\rm e}^{-\alpha(n_{2}^{2}-\tilde{p}^{2})})(1-{\rm e}^{-\alpha|\bm{n}_{1}-\bm{n}_{2}|^{2}})\right]\\ \mathcal{X}_{\rm int}(n_{1}^{2},\tilde{p^{2}})&=\int{\rm d}^{3}n_{2}\frac{\tilde{m}}{\tilde{m}^{2}+n_{2}^{2}}\frac{1-{\rm e}^{-\alpha(n_{2}^{2}-\tilde{p}^{2})}}{\tilde{\omega}_{n_{2}}-\tilde{\omega}_{p}}\frac{1-{\rm e}^{-\alpha|\bm{n}_{1}-\bm{n}_{2}|^{2}}}{|\bm{n}_{1}-\bm{n}_{2}|^{2}}-\alpha\frac{\tilde{m}}{\tilde{m}^{2}+n_{1}^{2}}\frac{1-{\rm e}^{-\alpha(n_{1}^{2}-\tilde{p}^{2})}}{\tilde{\omega}_{n_{1}}-\tilde{\omega}_{p}}.\end{split}$$

(35)

In the expression of $\mathcal{X}_{\rm int}$, the second term removes the value at the pole $\bm{n}_{2}=\bm{n}_{1}$ when replacing $\sum_{\bm{n}_{2}}\rightarrow\int{\rm d}^{3}n_{2}$. We use $\alpha=0.01$ and truncate the integer Cartesian coordinates at $|n_{x}|,|n_{y}|,|n_{z}|\leq 32$ in the present calculations. Under this condition, we evaluate the geometric constants $\mathcal{X}_{2}$ with a single sum and $\mathcal{R}_{24}$ with double sums, as defined in Eq. (A1) of Ref. [93], by using the TSS method. We obtain $\mathcal{X}_{2}=91.18$ and $\mathcal{R}_{24}=170.9$, which agrees with the corresponding results of Ref. [93] up to four significant figures.

For the calculations of the matrix element $\mathcal{T}_{L}^{(M)}$ between the bound states, they can be straightforwardly calculated using Eq. (28), once the bound-state wave function $\phi_{E,L}$ is solved. When solving the bound-state wave function, we truncate the integer sum at $|n_{x}|,|n_{y}|,|n_{z}|\leq\Lambda L/(2\pi)$. The number of momentum modes could reach several thousand, making direct diagonalization intractable. Therefore, we use the imaginary-time propagation starting from an initial wave function $\phi_{i}$ to solve for the bound-state wave function

$$\begin{split}\phi_{E,L}=\lim_{N_{\tau}\rightarrow\infty}({\rm e}^{-H\Delta\tau})^{N_{\tau}}\phi_{i},\end{split}$$

(36)

where $\Delta\tau$ is a small imaginary-time step, and ${\rm e}^{-H\Delta\tau}$ is expanded up to $O(\Delta\tau^{2})$. This is particularly efficient because of the separable form of the potential,

$$H\phi(\bm{p})=2\omega_{p}\phi(\bm{p})+C_{\Lambda}f_{\Lambda}(p^{2})\frac{1}{L^{3}}\sum_{\bm{k}\in\frac{2\pi}{L}\mathbb{Z}^{3}}f_{\Lambda}(k^{2})\phi(\bm{k}),$$

(37)

so the numerical complexity scales linearly with the number of momentum modes, instead of cubicly when using direct diagonalization. We take the initial wave function to be $\phi_{i}(\bm{p})={\rm e}^{-p^{2}/m_{N}^{2}}$, but using a different form should not affect the final results once the imaginary-time projection converges.

At the leading order of relativistic pionless EFT, there are two LECs $C_{\Lambda}$ and $\tilde{g}_{\nu}^{NN}$ that need to be determined for predicting the scattering amplitudes and matrix elements for the $nn\rightarrow ppee$ process. For $C_{\Lambda}$, it is the strength of the short-range strong potential and can be fixed by one low-energy observable in the ${}^{1}S_{0}$ channel. The observables used to determine $C_{\Lambda}$ are shown in Table 1. At the physical pion mass, there is no ${}^{1}S_{0}$ bound state, so we fix it using the experimental scattering length $a=-23.74$ fm. At the unphysical pion masses, several LQCD calculations yielded deeply bound ${}^{1}S_{0}$ two-nucleon state [75, 76, 73, 74]. However, many other LQCD studies [94, 95, 85, 77, 86] did not obtain such bound states, raising concerns on whether or not the previous works correctly determined the two-nucleon spectrum. There are several explanations for this issue and it is still not completely conclusive whether ${}^{1}S_{0}$ two-nucleon state is bound or unbound at unphysical large pion masses [96]. Nevertheless, the results from the present work could be easily adjusted to updated LQCD values of two-nucleon binding energies or scattering lengths. For example, in Table 1, we considered both the older (without asterisk) [74] and the latest (with asterisk) [77] LQCD values of $B_{nn}$ at $m_{\pi}=806$ MeV. For $m_{\pi}=300$, 400, 510, and 806 MeV, the LEC $C_{\Lambda}$ is fixed using the two-nucleon binding energy $B_{nn}$ in the infinite volume, extrapolated from the FV energies provided by the LQCD calculations. For $m_{\pi}=806^{*}$ MeV, the LQCD result of $B_{nn}$ is only available at a single finite volume $L=4.6$ fm and, thus, the fitting of the LEC $C_{\Lambda}$ is performed at the same finite volume.