Renormalizability of the leading order operator for neutrinoless double beta decay with the effects of finite nucleon size

Introduction.—Neutrinoless double beta (0$\nu\beta\beta$) decay, a hypothetical second-order weak transition that converts two neutrons in a nucleus into two protons through the emission of two electrons but without neutrinos Furry (1939), is an ideal venue for probing lepton number violation (LNV). Besides demonstrating that neutrinos are Majorana fermions Schechter and Valle (1982), its observation provides insightful perspectives to explain the neutrino mass generation Mohapatra and Senjanovic (1980), determines the absolute neutrino mass scale and hierarchy Avignone et al. (2008), and explores the mechanism of the matter-antimatter asymmetry in early universe Davidson et al. (2008). Currently, numerous underground experiments with high sensitivity Ni et al. (2019); Agostini et al. (2020); Giachero et al. (2021); Abe et al. (2023) all over the world are actively searching for 0$\nu\beta\beta$ decay signals for different even-even nuclei. In the future, more ton-scale facilities Abgrall et al. (2017); Albert et al. (2018); Armatol et al. (2022); Cao et al. (2024) will be deployed to improve the discovery potential of experiments and measure the half-life of different candidate isotopes (for a recent review, see Ref. Agostini et al. (2023)).

These experiments will help us make better understanding of this process. In general, whatever operator including the LNV physics may induce the nuclear 0$\nu\beta\beta$ decay. To reveal the underlying mechanism of this process, the precise calculation of nuclear matrix elements (NMEs) that incorporates various LNV contributions is essential. However, its estimation based on the different many-body methods suffers from a discrepancy of a factor of three to five thus far Engel and Menéndez (2017). One should note that the uncertainty of NMEs is rooted in the realization of nuclear structure and LNV mechanisms.

It has been shown that the chiral effective field theory (ChEFT) Weinberg (1990, 1991) offers a promising tool to quantify such uncertainties in nuclear calculations, for a comprehensive review see Ref. Hammer et al. (2020). Moreover, the electroweak currents can be introduced systematically in ChEFT Pastore et al. (2009); Baroni et al. (2016); Krebs et al. (2017, 2019). By using series of EFTs, the various LNV sources at energy well above the electroweak scale can be expressed as the 0$\nu\beta\beta$ operators at the nucleon level Prezeau et al. (2003); Cirigliano et al. (2017); Cirigliano et al. (2018b, c); Dekens et al. (2020); de Vries et al. (2022). With ChEFT, ab initio calculations under the standard light Majorana neutrino-exchange mechanism Weinberg (1979) for NMEs were performed Yao et al. (2020); Belley et al. (2021); Novario et al. (2021); Wirth et al. (2021); Belley et al. (2024).

Recently, within the standard mechanism of 0$\nu\beta\beta$ decay, the ChEFT description based on naive dimensional analysis (NDA) Cirigliano et al. (2018b) has been proven to violate the requirement of renormalizability Cirigliano et al. (2018a); Cirigliano et al. (2019). As a consequence, a contact operator that encodes the effect of hard-neutrino exchange is promoted to leading order (LO). The corresponding low-energy constant (LEC) is currently still undetermined by experiments or lattice QCD simulations, introducing additional uncertainties to the estimation of NMEs besides those from nuclear structure.

The prescription to estimate the size of the LEC was first presented by Cirigliano et al. in Cirigliano et al. (2021a, b). Although this approach by matching the $nn\to ppe^{-}e^{-}$ amplitude onto the Cottingham formula Cottingham (1963), is somewhat model-dependent, the uncertainty is acceptable at LO calculation for ChEFT. On the other hand, Yang et al. Yang and Zhao (2023) developed a relativistic method that can renormalize the LO decay amplitude by resorting to the light neutrino exchange rather than utilizing a contact term. The two distinct approaches show decent consistency in the prediction of the LO amplitude.

An essential ingredient used to evaluate the 0$\nu\beta\beta$ amplitude is the coupling of lepton-neutrino current to a nucleon weak current. This coupling can be decomposed into all structures permitted by Lorentz invariance, whose kinematics is factorized into the momentum-dependent functions, i.e., nucleon form factors (FFs), encoding the effects of finite nucleon size. As an application of ChEFT for the low-energy electroweak processes for nucleons, the nucleon FFs can be expressed in terms of powers of momentum (see Ref. Bernard et al. (1995) for a review), in which the LO term corresponds to the point-like assumption of the nucleon. It is important to emphasize that the prerequisite for this perturbative expansion is the observability of the external source particles. If one of the particles is virtual, the nucleon FFs will vary based on the momentum of the intermediate state. From this perspective, a non-perturbative approach to incorporate nucleon structure effects may provide novel insights into the renormalization of the 0$\nu\beta\beta$ amplitude. Furthermore, it is worth stressing that the aforementioned nucleon FFs are inappropriate for assessing the impact of nucleon structure on nuclear forces because nuclear forces are so strong that the internal structure of the nucleon changes during strong process.

In this study, we revisit the $nn\to ppe^{-}e^{-}$ transition mediated solely by the light Majorana neutrinos to examine the renormalizability of 0$\nu\beta\beta$ operators. The finite size of the nucleon is taken into account in a non-perturbative manner, similar to methodologies employed in nuclear many-body calculations Simkovic et al. (1999, 2009); Hyvärinen and Suhonen (2015); Yao et al. (2022). Using the momentum-dependent FFs, we obtain the finite and regulator-independent amplitude at LO in the absence of short-range contribution. It shows the model-free advantage of this method and the validity of the way to treat nucleon structure in conventional nuclear calculation.

Theoretical framework.—We focus on the scenario of light Majorana neutrino mediated $nn\to ppe^{-}e^{-}$ transition. The scattering of initial and final state particles is restricted to the ${}^{1}S_{0}$ channel because it has been proved in Refs. Cirigliano et al. (2018a); Cirigliano et al. (2019) that the conflict between NDA and renormalizability only appears in this channel. For this purpose, the heavy-baryon formulation Jenkins and Manohar (1991); Bernard et al. (1992) of ChEFT is employed and the Coulomb interaction is ignored.

We start with the chiral Lagrangian, the nucleon-nucleon potential derived from two-nucleon irreducible diagrams comprises the long-range pion exchange and the short-range four-nucleon contact terms. At LO, the two-nucleon potential is given by

where $g_{A}=1.27$ is the nucleon axial coupling, $f_{\pi}=92.2$ MeV the pion decay constant, $m_{\pi}=138$ MeV the pion mass. The initial and final relative momenta, denoting as $\boldsymbol{p}$ and $\boldsymbol{p}\,^{\prime}$ respectively, define the momentum transfer $\boldsymbol{q}=\boldsymbol{p}\,^{\prime}-\boldsymbol{p}$. In the ${}^{1}S_{0}$ channel, the combination of LECs $C_{S}$ and $C_{T}$, i.e. $\widetilde{C}=C_{S}-3C_{T}$, is needed.

In the context of the standard mechanism of 0$\nu\beta\beta$ decay, the electron-neutrino Majorana mass term, which is the dominant contribution of LNV at low energy, can be written as

where $C=i\gamma_{2}\gamma_{0}$ denotes the charge conjugation operator, and the effective mass $m_{\beta\beta}=\sum_{i}U^{2}_{ei}m_{i}$ is defined with masses of three generations of light neutrino and elements of the neutrino mixing matrix. The definition of the neutrino potential is similar to that of the two-nucleon potential Cirigliano et al. (2017). At LO in NDA, only following neutrino-exchange potential

contributes to the $nn\to ppe^{-}e^{-}$ amplitude. In light of the non-perturbative property of the nucleon structure in the 0$\nu\beta\beta$ process, we have introduced the momentum-dependent FFs $G_{V}\left(\boldsymbol{q}^{2}\right)$ and $G_{A}\left(\boldsymbol{q}^{2}\right)$, which are normalized as the vector and axial-vector coupling respectively, to take into account the finite-nucleon-size effects.

The commonly used dipole parametrization $G_{\alpha}\left(\boldsymbol{q}^{2}\right)=g_{\alpha}/(1+\boldsymbol{q}^{2}/\Lambda_{\alpha}^{2})^{2}$ ($\alpha=V,A$) is employed in this work. Here we adopt $\Lambda_{V}=0.84$ GeV and $\Lambda_{A}=1.05$ GeV respectively. Given that the prediction of dipole parametrization is inconsistent with data in the high-momentum region Jones et al. (2000); Gayou et al. (2001, 2002); Puckett et al. (2010); Meyer et al. (2016), we additionally utilize the parametrization scheme presented by Kelly Kelly (2004) to describe the structure of the nucleon.

The Kelly parametrization of the axial FF takes the form

where $x=\boldsymbol{q}^{2}/(4m_{N}^{2})$ with the nucleon mass $m_{N}=939$ MeV, the dimensionless model parameters $a_{1},b_{1},b_{2}$ and $b_{3}$ are determined by fitting of data from Ref. Hashamipour et al. (2019). On the other hand, in conjunction with the fits in Refs. Kelly (2004); Punjabi et al. (2015) for the Sachs FFs of proton and neutron, one can obtain the vector FF as a function of momentum via the relation

where $G_{Ep}$, $G_{Mp}$, $G_{En}$ and $G_{Mn}$ are the so-called electric and magnetic Sachs FFs, for proton and neutron, respectively.

The results of the dipole and fitted Kelly forms are shown in Fig.1.

The prediction of the Kelly parametrization for the nucleon FFs is in reasonable agreement with the data. One of the important points is that both types of parametrization schemes for vector and axial FFs all predict a general behavior:

By this scale rule of FFs, we will demonstrate the convergence of the 0$\nu\beta\beta$ amplitude in the ultraviolet (UV) region. Comparing the transition amplitudes based on both parametrizations, we will demonstrate the validity of assessing the effects of finite nucleon size with dipole parametrization on NMEs.

The LO 0$\nu\beta\beta$ amplitude based on the kinematics $n(\boldsymbol{p})n(-\boldsymbol{p})\to p(\boldsymbol{p}^{\prime})p(-\boldsymbol{p}^{\prime})e^{-}\left(\boldsymbol{p}_{e1}=0\right)e^{-}\left(\boldsymbol{p}_{e2}=0\right)$ can be schematically expressed as in Cirigliano et al. (2019):

where $G_{0}$ is the two-nucleon free propagator, $T_{s}$ the two-nucleon scattering amplitude obtained by iterating the LO strong potential $V_{NN}$ through $G_{0}$. We name the four terms in Eq.(7), from left to right, $A_{\mathcal{A}}$, $A_{\mathcal{B}}$, $\bar{A}_{\mathcal{B}}$, and $A_{\mathcal{C}}$.

We are interested in the renormalizability of the 0$\nu\beta\beta$ amplitude. From Eq.(7), the sub-amplitude $A_{\mathcal{A}}$ as a tree-level contribution must be convergent. Because of the finiteness of the strong amplitude $T_{s}$, the UV behavior of the remaining terms is dominated by the loop integrals with the insertion of neutrino potential. The nonrelativistic propagator $G_{0}$ behaves like $\mathcal{O}\left(\Lambda^{-2}\right)$ in the UV region, where $\Lambda$ denotes a UV energy scale much greater than the typical momentum in the process. For the neutrino potential in the ${}^{1}S_{0}$ wave, its UV scale rule is given by

where $g=0$ (4) corresponds to the case without (with) the effects of finite nucleon size. Therefore, in the nonrelativistic case, the superficial degree of divergence for the last three terms in Eq.(7) has the following form:

where $L$ is the number of loops, and $N=1$ is the number of insertions of the neutrino potential. Given that $L=1$ and $g\geq 0$ for $A_{\mathcal{B}}$ and $\bar{A}_{\mathcal{B}}$, it follows that the both sub-amplitudes could reach the convergence as $\mathcal{O}\left(\Lambda^{-1-2g}\right)$. However, $A_{\mathcal{C}}\sim\mathcal{O}\left(\Lambda^{-2g}\right)$ means that its convergence depends on the treatment of the nucleon structure. In the case without the effects of finite nucleon size, i.e. $g=0$, $A_{\mathcal{C}}$ diverges logarithmically such that a LEC $g_{\nu}^{NN}\sim\mathcal{O}\left(f_{\pi}^{-2}\right)$ is required to renormalize the LO amplitude $A_{\nu}$ Cirigliano et al. (2018a); Cirigliano et al. (2019). With finite nucleon size, however, $A_{\mathcal{C}}\sim\mathcal{O}\left(\Lambda^{-8}\right)$ is finite so that $A_{\nu}$ can be renormalized without inclusion of any short-range counter terms.

Results and discussions.—Following the usual numerical strategy, the renormalizability of $A_{\nu}$ is demonstrated by its independence of the cutoff value in momentum or coordinate space. In this work, we introduce the momentum cutoff $\Lambda$ via the way that regulates the UV part of the LO strong potential with a separable Gaussian regulator:

To determine the LEC $\widetilde{C}$ in the strong interaction, we fit $T_{s}$ to the ${}^{1}S_{0}$ wave neutron-proton scattering length $a_{np}=-23.74$ fm. The cutoff independence of the phase shift has been achieved at $\Lambda\to\infty$.

Choosing the kinematic points $|\boldsymbol{p}|=1$ MeV and $|\boldsymbol{p}^{\prime}|=38$ MeV, we depict in Fig.2 the cutoff variation of the LO amplitude $A_{\nu}$.

The ChEFT results without the contributions of nucleon structure display a logarithmic divergence, while the amplitude with the finite-nucleon-size effects shows a satisfactory convergence, as same as the expectation above. In addition, the amplitudes obtained from both the dipole and Kelly parametrizations are remarkably similar at any cutoff value, with a difference around only 1%. This indicates that the impact of different parametrizations of nucleon structure on 0$\nu\beta\beta$ decay is negligible, highlighting the robustness and reliability of the parametrization in calculating the 0$\nu\beta\beta$ amplitude. Therefore, the dipole parametrization is accurate enough to assess the contribution of the nucleon structure to NMEs.

To validate the present scheme, we calculate the amplitude for $|\boldsymbol{p}|=25$ MeV and $|\boldsymbol{p}^{\prime}|=30$ MeV and perform a comparison with the previous results Cirigliano et al. (2021a, b); Yang and Zhao (2023). As shown in Fig.3,

again, we capture the renormalized amplitude without introducing any undeterminate contact operators. The limit value of absolute amplitude $|A_{\nu}|=0.01956(8)$ $\text{MeV}^{-2}$ at this kinematic point, where the error stems from the dependence on parametrization. Although the approach used to renormalize $A_{\nu}$ in this study differs from that in Refs. Cirigliano et al. (2021a); Yang and Zhao (2023), our results are consistent with the estimations of Refs. Cirigliano et al. (2021a); Yang and Zhao (2023) within the uncertainty 10%, a relatively strict upper bound for ChEFT truncation error at LO.

The renormalized on-shell and off-shell amplitudes as functions of momentum are shown in Fig.4. At various kinematic points, the results accounting for the effects of nucleon structure are in quite good agreement with the predictions based on the renormalized scheme of Refs. Cirigliano et al. (2021a, b).

Furthermore, the present scheme can also recover the transition amplitude for higher partial-wave channels obtained in Ref Cirigliano et al. (2019).

The primary motivation for employing a non-perturbative approach to treat the finite-size effects of the nucleon arises from the inadequacy of perturbative expansion of the nucleon FFs caused by virtual neutrinos, rather than from the behavior of FFs decreasing as $\mathcal{O}(\Lambda^{-4})$ in the UV region. The numerical results indicate that the approach to including the suppressed effects of FFs is equivalent to the renormalized scheme of Ref. Cirigliano et al. (2021a) at LO. Therefore, to ensure the correctness of the result, we argue that the calculation of the 0$\nu\beta\beta$ amplitude should be conducted within a uniform renormalization scheme.

We note that conventional many-body calculations estimate the contributions of finite nucleon size to NMEs in the same manner as this work. Based on the above discussions, we conclude that these calculations introducing the nucleon structure effects in a non-perturbative way are consistent with the scheme that renormalizes the amplitude with short-range terms. In contrast, combining the nucleon structure with the short-range contribution in many-body calculation may lead to an overestimation of the amplitude, thereby amplifying the uncertainty of NMEs.

Moreover, since we evaluate the LO 0$\nu\beta\beta$ amplitude without relying on any unknown LECs, our approach provides a model-independent description for 0$\nu\beta\beta$ decay. This is crucial for improving the precision of NMEs calculated in many-body methods. Besides circumventing model dependence, the realistic size of the short-range operator as discussed in Refs. Cirigliano et al. (2018a); Cirigliano et al. (2019); Cirigliano et al. (2021a, b) can be addressed by matching the estimation of the decay rate based on the present scheme.

Our framework is more straightforward for advancing higher-order calculations, in contrast to the relativistic scheme Yang and Zhao (2023) within the framework of manifestly Lorentz-invariant ChEFT. The latter would introduce certain complexities when moving to higher orders, especially for the treatment of the renormalization Hammer et al. (2020); Epelbaum et al. (2020). Meanwhile, our method facilitates a more intuitive combination with existing many-body techniques compared to the relativistic scheme.

Summary.—In this study, we present a nonrelativistic approach to renormalize the LO 0$\nu\beta\beta$ amplitude for the standard mechanism using the effects of nucleon structure instead of the short-distance operators with undetermined LECs. Our calculations of the 0$\nu\beta\beta$ amplitude at various kinematic points demonstrate the efficacy of the proposed scheme, showing its comparability with established frameworks Cirigliano et al. (2021a, b). The most significant implication of this equivalence is the validation of conventional many-body calculations based on the light Majorana neutrino mechanism, confirming the absence of any oversight regarding LO LNV physics in these calculations. Finally, our scheme has the capacity to transcend the specific scenario of light Majorana neutrino exchanges, generalizing to the non-standard decay mechanisms and helping to reduce the uncertainties of corresponding NMEs in many-body calculations.