Neutrinoless double beta decay and QCD running at low energy scales

There is a common belief that the main uncertainties in the theoretical analysis of neutrinoless double beta ($0\nu\beta\beta$) decay originate from the nuclear matrix elements. Here, we uncover another previously overlooked source of potentially large uncertainties stemming from non-perturbative QCD effects. Recently perturbative QCD corrections have been calculated for all dimension 6 and 9 effective operators describing $0\nu\beta\beta$-decay and their importance for a reliable treatment of $0\nu\beta\beta$-decay has been demonstrated. However, these perturbative results are valid at energy scales above $\sim 1$ GeV, while the typical $0\nu\beta\beta$-scale is about $\sim 100$ MeV. In view of this fact we examine the possibility of extrapolating the perturbative results towards sub-GeV non-perturbative scales on the basis of the QCD coupling constant"freezing"behavior using Background Perturbation Theory. Our analysis suggests that such an infrared extrapolation does modify the perturbative results for both short-range and long-range mechanisms of $0\nu\beta\beta$-decay in general only moderately. We also discuss that the tensor$\otimes$tensor effective operator can not appear alone in the low-energy limit of any renormalizable high-scale model and then demonstrate that all five linearly independent combinations of the scalar and tensor operators, that can appear in renormalizable models, are infrared stable.


I. INTRODUCTION
From neutrino oscillation experiments it is nowadays well-known that at least two neutrinos have non-zero masses. Oscillation experiments, however, can not decide whether neutrinos are Dirac or Majorana particles. Lepton number violating (LNV) processes have to be studied to explore the neutrino nature [1][2][3]. Neutrinoless double beta decay (0νββ), for recent reviews see for instance [4][5][6], is by far the most powerful available probe of lepton number violation, and its non-observation allows to constrain LNV beyond standard model (BSM) physics. From several experimental searches for 0νββ [7][8][9][10] the best current bound on the 0νββ half-life, T 0ν 1/2 , comes from the KamLAND-Zen experiment [8]: From the theoretical perspective, there are two different kind of contributions to the 0νββ amplitude: the short-range mechanisms (SRM) [11], which are mediated by heavy particle exchange; and the long-range mechanisms (LRM) [12], in which a light neutrino is exchanged between two point-like vertices. It has been recently demonstrated that QCD corrections to 0νββ are important, especially in the SRM case [13] due the presence of the color-mismatch effect and the corresponding mixing of different operators, with numerically very different nuclear matrix elements (NME). This effect leads to differences in the limits on the Wilson Coefficients (WC), which amount in some cases up to 3 orders of magnitude. On the other hand, the LRM operate between two different and distant nucleons, so that no color-mismatch appears and only QCD vertex corrections have to be taken into account. Their effect on the extracted limits does not exceed 60% [14], less than the typical estimate of the uncertainties of the nuclear matrix elements (NMEs).
These QCD RGE results [13,14] are valid for energy scales above ∼ 1 GeV -the limit of perturbative QCD. On the other hand, the typical scale of 0νββ-decay is about ∼ 100 MeV. It is then natural to ask, if these perturbative results still can be considered a reasonable approximation or whether they could be drastically affected by some non-perturbative effects. In this paper, we try to address this issue within the semi-empirical approach of "freezing" the running of the strong coupling constant at low energies. Freezing is motivated by some approaches to non-perturbative QCD effects, in particular, by the Background Perturbation Theory (BPTh) [20]. We do not pretend that our work will determine the final conclusion about this question, but rather hope that our analysis will allow us at least to visualize this problem, overlooked in the literature, and stimulate efforts for theoretically grounded studies in this direction. The conventional approach would rely on the operator product expansion for observables, such as 0νββ half-life, and on a proper matching of the quark and nucleon level theories at a certain scale µ 0 . This is a low-energy scale, down to which perturbative QCD for the quark-level theory is valid. In a self-consistent approach both Wilson coefficients and the nucleon matrix elements of the effective operators depend on µ 0 in such a way that the effect in any observable should cancel, in theory. Unfortunately, the nucleon matrix elements of 0νββ-operators have not been studied yet from this perspective. The hadronization prescriptions existing in the literature (see, for instance, Refs. [21,22]) parameterize the matrix elements in terms of phenomenological nucleon form-factors loosing connection with QCD and, therefore, a dependence on the matching scale µ 0 remains. Presumably, lattice QCD is able to shed light on this issue. Some recent lattice QCD publications about matrix elements of 0νββoperators can be found in Refs. [23][24][25]. Recently an approach to hadronization in 0νββ-decay, which probably suits better for the matching with lattice QCD, has been developed in Refs. [26,27] on the basis of chiral effective field theory. However, the issue of the matching-scale dependence still remains unaddressed. This is the situation in which we adopted the idea of "freezing". In a sense it can be thought of as a rough modeling of the matching scale-dependence of the nucleon matrix elements alleviating the dependence of the 0νββ half-life on this scale.

II. QCD RUNNING COUPLING CONSTANT IN THE INFRARED LIMIT
The high-energy behavior of QCD is well understood. Thanks to asymptotic freedom, perturbation theory allows an accurate calculation of the running QCD coupling constant α s (µ). In the infrared region, however, non-perturbative effects become important and lead, in particular, to color confinement at some energy scale below ∼ 1 GeV [15][16][17]. At such low energies, perturbative QCD (pQCD) suffers a singularity, the so-called Landau pole, i.e. α s (µ) tends to infinity. This singularity prohibits even a formal extrapolation of the pQCD results to the region below ∼ 1 GeV. Unsurprisingly, in the literature there is a plethora of approaches aiming at the extrapolation of the pQCD results towards the infrared (IR) region taking into account at least some non-perturbative effects suggested by QCD. We will not discuss all of these attempts and instead refer to the recent review [17] on this subject.
However, different authors have discussed the possibility, for a list of references see for example [17], that the predictions of QCD can be formulated in terms of a modifiedα s (µ), which is finite in the IR regime. Basically, the strong increase of α s , predicted by pQCD, is thought to be regularized by non-perturbative effects summed up in an effectiveα s (µ), which stops growing at some infrared scale. This IR behavior, dubbed "freezing" of α s in the literature is compatible with the lattice QCD simulations [18,19].
To give an example let us consider the Background Perturbation Theory (BPTh) [20] applied to QCD. In this approach the gluon field A µ is separated into two parts A µ = a µ + B µ , the perturbative part a µ , and an effective non-perturbative background field B µ . The BPTh one-loop QCD running coupling constant has been calculated in Refs. [28,29]. It has been shown that the non-perturbative background field generates an effective mass M B for the gluon in such a way that the argument of α s (µ 2 ) is replaced by µ 2 → µ 2 + M 2 B . In the UV limit the effects of the background field become negligible, while at low energies it essentially leads to a cutoff in the running. Then, the modified one-loop QCD running coupling constant can be written as: Here, β 0 is the usual β-coefficient, while the mass parameter, M B , is non-universal and depends on the specific process in consideration [29,30]. The parameter M B can be related to the mass of the glueball M 2g (0 ++ ) [29], a bound-state of two gluons connected with the adjoint string In Refs. [31,32] α s has been calculated in the BPTh up to 3 loops and the value M B = 1.06 GeV has been extracted from the analysis of bottomium fine structure. We mention also [33], quoting M B = 0.95 GeV, and [28], which give M B = 2.15 GeV. Below we revisit the QCD corrections to the SRM and LRM of 0νββ-decay [13,14] applying Eq. (2) in order to extrapolate the pQCD results towards µ ∼ 100 MeV, the characteristic scale of 0νββ-decay. In order to estimate the uncertainty of our results, we will not use some specific value of M B and instead show our results as functions of the calculated value of the "frozen" value of α F S , corresponding toα s (µ 2 ) for µ ≤ 0.1 GeV. Although not numerically relevant, for completeness we mention that for the normalization of the expression in Eq. (2), we use the experimental value α s (µ = M Z ) from [34].

III. QCD RUNNING OF WILSON COEFFICIENTS
Here we apply the IR extrapolation based on "freezing", discussed above, to the pQCD results for 0νββdecay derived in refs [13,14,35].

A. pQCD Wilson Coefficients
First let us recall the pQCD results for the short-range and long-range mechanisms of 0νββ-decay derived in [13,35].
(A) Short-range mechanisms of 0νββ-decay are mediated by the exchange of heavy particles of mass M I between two nucleons of the decaying nucleus. At low energy scales µ < M I this mechanism is described the low-energy effective Lagrangian [11,13]: with the complete set of dimension-9 non-equivalent 0νββ-operators [13,35] O (9)XY 1 where X, Y = L, R and the leptonic currents are The QCD corrections induce operator mixing due to the color-mismatch effect, which is equivalent to mixing different Wilson coefficients C XY i (µ) in Eq. (3) in the amplitude. This effect has a dramatical impact on the numerical predictions of high-scale models for the 0νββ-decay half-life [13,35] due to the significant differences in the numerical values of NMEs of different operators (4)- (8). The QCD-corrected 0νββ-decay half-life formula, reads [13] T 0νββ [11,38], and summation over the different chiralities X, Y = L, R is implied. The parameters β (SRM ) i are defined in [13] and represent linear combinations of the NMEs of the operators (4)- (8).
(B) Long-Range Mechanisms describe contributions via light neutrino exchange between two quarks belonging to two different and distant nucleons of a decaying nucleus. One of the vertices connected by the neutrino propagator is the standard V −A vertex while the other is an LNV ∆L = 2 beyond the SM (BSM) vertex originating from some heavy particle exchange. At low energies the BSM vertices are given in terms of the effective operators [35] O with X = R, L. Then 0νββ-decay amplitude is given by the second-order of perturbation theory in the effective Lagrangian [12,35]: Here the first term is the SM low-energy 4-fermion effective interaction of the currents The QCD-corrected 0νββ-decay half-life formula for the LRM reads [14] T 0νββ where G 0i and (NME) i are the phase-space factors [36]. The coefficients β (LRM ) i (µ 0 , Λ) involve NMEs of the operators (11)-(13) and the pQCD RGE running parameters. Unlike the analogous β-coefficients in Eq. (10) they do not mix the NMEs of different operators. This is a drastic difference between the SRM and LRM caused by the absence of the color-mismatch in the latter case as demonstrated in Ref. [35].

B. "Freezing" Wilson Coefficients
The Wilson coefficients C i in Eqs. (3) and (14) are calculable in the framework of a particular underlying model above some high-energy scale Λ. These coefficients can be expressed in terms of the model couplings and heavy particle masses (for a comprehensive analysis see Ref. [35]). At the scale Λ the model is matched to the effective theory given by the Lagrangians (3) and (14). Then the Wilson coefficients C i should be RGE evolved, taking into account the QCD loop-corrections, down to the characteristic scale µ 0 of 0νββ, typically about 100 MeV. Refs. [13,14] stopped the RGE evolution at 1 GeV. This kind of truncation is a common practice in the literature applying pQCD to observable processes. We now will discuss the crucial question if the pQCD results are stable on the way from 1 GeV to 100 MeV.
The discussion in section (II) suggests a method for a rough assessment of the effect of the IR extrapolation by simple freezing the one-loop pQCD running coupling. This amounts to the replacement of α s (µ) byα s (µ) from Eq. (2) in the expressions for all the coefficients β(µ 0 , Λ) in Eqs. (10) , (16). Sinceα s (µ) is finite in the IR limit, one can then extend the running down to the required value of µ 0 100 MeV.
Our main results are then shown in Fig. 1 for the 9 SRM coefficients and in Fig. 3 Fig. 3, because the relative changes of C L i and C R i are the same. Note also, that the coefficients are calculated for the particular case of 136 Xe, but the plots for other isotopes are very similar (since we plot ∆(C AB i )).
As Fig. 1 shows, all 9 short-range coefficients change only moderately, when α F S is varied in the window (0.2-0.8). The two C 3 are the most stable coefficients, while the other C i change within (less than) a factor of roughly two. We would like to point out that it is usually argued in the 0νββ decay community, that nuclear matrix elements (NME) have uncertainties of a factor of roughly 2 as well, thus we call the change of the coefficients under variation of α F S "moderate", since it is of comparable size.
Nevertheless, we would like to mention that for values of α F S > ∼ 1 larger changes of the Wilson coefficients result. Thus, the stability of our extrapolation rests on whether or not the idea of a finite, frozen value of α S is indeed

correct.
We would also like to discuss that among the operators in Eq. (5), forming the low-energy operator basis, the tensor operators O are special in the sense that they can never appear alone, i.e. without O (9) 1 , in the low-energy limit of any renormalizable high-scale model (HSM). This is because renormalizable models do not contain fundamental tensor interactions. Then, at low energies the effective operators O can arise only from a Fierz transformation of the scalar operator O (9) 1 in the Lorentz and color indices. As was shown in [39] there are 5 possible linear combinations of these two operators, which can originate from renormalizable high-scale models (for details see Ref. [14]) 1 + c 2 O 1 + c2O that appear in the tree-level decomposition of the 0νββ decay operator, see [39]. For a discussion see text. We show the five combinations of coefficients c 1 , c 2 in Fig. 2 together with the corresponding limits on these combinations, again as a function of α F S . The plot demonstrates that these limits are stable (within again roughly a factor 2) with respect to variations of α F S for all 5 cases.  Fig. 1. Thus, one can conclude that both, running and freezing, do not play an important role for the long-range part of the 0νββ decay amplitude.

IV. DISCUSSION AND CONCLUSIONS
We have calculated QCD corrections to 0νββ for both cases, the SRM and LRM, with particular emphasis on the IR behaviour of the QCD running coupling. We discussed the regulated form of the QCD running coupling, using Background Perturbation Theory, which introduces the background mass parameter M B . In this treatment, the resultingα S (Q 2 ) "freezes" at values of Q 2 smaller than Q 2 ≤ 0.1 GeV. However, since M B in this setup is both, a model-dependent as well as a process-dependent parameter, the exact value of α F S is not predicted. We, therefore, showed our results as a function of α F S . We conclude that the Wilson coefficients depend only moderately on the exact value of α F S and it seems we can extract reliable limits on these coefficients from 0νββ. We noted that in renormalizable ultra-violet completions ("high scale models") the tensor operator is never generated alone, ie. it always appears in combination with O 2 that appear in the (tree-level) decomposition of the 0νββ decay operator are also stable in the above quoted sense.
On the other hand, as a word of caution we need to mention that the idea of a frozen α S at low energy is not at present universally accepted, see the detailed discussion in [17]. Our conclusions will remain valid if α F S can be shown to lie definitely below, say, α F S < ∼ (0.7 − 0.8). For values of α F S much larger than this, our simpleminded treatment can not be considered to be valid.
Finally, we would like to stress again that the present study does not pretend to present a well theoretically grounded analysis of the non-perturbative QCD effects in the transition region from the quark to the nucleon level description of 0νββ-decay. The problem of the matching, at some low-energy scale, of the perturbative QCD results with the effective low-energy theory in terms of the nucleon fields remains open and is waiting for its clarification. Future work in this direction would be very valuable.