Evaluation of neutrinoless double beta decay: QCD running to sub-GeV scales

We evaluate QCD effects in the neutrinoless double beta ($0\nu\beta\beta$) decay, originating from new physics short-range mechanism in the form of five dimension-9 operators. For this, we employ the one-loop and two-loop renormalization group equations (RGEs) for the corresponding Wilson coefficients, performing the RGE-evolution from the new physics scales (estimated as $\Lambda ~ \sim 10^2$ GeV) to the typical spacelike $0\nu\beta\beta$-scale $Q \sim 0.1$ GeV. Since the latter scale is clearly nonperturbative, we apply various infrared-safe (IR-safe) variants of QCD where the running coupling has no Landau singularities at low spacelike $Q$. We point out that the correct treatment of the IR-safe analogs of the (noninteger) powers of the couplings is important. It turns out that in most cases of the considered operators the resulting QCD effects can be significant in this process, i.e., can be stronger than the effects of the present uncertainties in the nuclear matrix elements.


I. INTRODUCTION
One of the basic questions of high energy physics is whether the neutrinos, and/or their more exotic fermionic relatives if they exist, are Majorana or Dirac particles. The question of the existence of Majorana neutrinos is closely related with the question of whether the lepton number violating (LNV) processes exist. At present, the most powerful probe of LNV processes is the neutrinoless double beta (0νββ) decay (cf. [1, 2] for recent reviews), i.e., the process where two d quarks of a nucleus transform into two u quarks with the simultaneous production of two low-energy electrons. Such processes have not (yet) been observed, and one of the best lower bounds on the half-life for 0νββ is from the KamLAND-Zen experiment [3] for the decay of 136 Xe decay process can be regarded as a low-energy spacelike process. This means that the half-life T 0ν 1/2 ≡ D(Q 2 ) can be regarded as a spacelike observable with positive Q 2 (≡ −q 2 ) ∼ µ 2 f ∼ 10 −2 GeV 2 . Furthermore, if the mass of the exchanged neutrino N is heavy (M N 0.1 GeV), the process can be regarded as an effective pointlike process dd → uuee. It can be called a short-range process, due to the high masses of the exchanged particles involved. On the other hand, such short-range (pointlike) process dd → uuee can originate also from an exotic physics [4] which can be described effectively in terms of dimension-9 operators where the scale of the new LNV-physics is expected to be Λ LNV 10 2 GeV. There are five classes of such effective (pointlike) operators (see the next Section). Since there is a very large difference between the new physics scale Λ LNV and the 0νββ decay scale µ f , the effects of the QCD corrections on the corresponding Wilson coefficients (which  [49]), at effective Fermi motion scales Q 2 f = µ 2 f ∼ 0.01 GeV 2 .

II. EFFECTIVE LAGRANGIAN IN 0νββ DECAY
The effective Lagrangian within the Operator Product Expansion (OPE) formalism for the dimension-9 operators, which originate from short-range new physics and contribute to 0νββ decay, have the generic structure [48] L 0,νββ where G F = 1.166 × 10 −5 GeV −2 is the Fermi constant, m p is the proton mass. The expansion (3) contains five types of dimension-9 operators; the indices X, Y = L, R indicate the chirality. The dimension-9 operators O XY i can be shown in the compact notation [49,50] O XY where j =ē (1 ± γ 5 ) e c , j µ =ēγ µ γ 5 e c are the lepton currents. In Eqs. (4b) and (4d) we use the convention σ µν = (i/2)[γ µ , γ ν ]. In general, these operators mix under renormalization through QCD when we express them in terms of a color singlet structure. In this procedure, the following property is used: (λ a ) αβ (λ a ) ηξ = − 2 N δ αβ δ ηξ + 2δ αξ δ ηβ , where λ a = 2t a are the Gell-Mann matrices. This leads to the original operator [the first term on the RHS of Eq.(5)] plus a color mismatch part [the second term on the RHS of Eq. (5)]. Note that O XY 2 = 0 for X = Y . The effective Lagrangian (3) at high physics scales µ = Λ LNV (∼ 10 2 -10 3 GeV) represents the new short-range physics. When these contributions are evolved to lower scales µ, the QCD effects are the dominant contributions to the RGE evolution. The effective Lagrangian (3) must be evaluated down to a spacelike scale µ 2 = Q 2 (≡ −q 2 ) that enters in the 0νββ process (and in the corresponding NMEs), typically of the order of the Fermi motion scale µ ∼ µ f ∼ 0.1 GeV. In practice, when we use the Lagrangian (3) in the perturbation theory within pQCD, it is applicable only down to µ ∼ 1 GeV in the best scenario. This RGE-running for 0νββ decay was performed in pQCD, at one-loop level of anomalous dimensions, in [51] for the set of operators O 1 -O 3 , and in [49] for the set O 1 -O 5 . The restriction µ 1 GeV is due to unphysical singularities, known as Landau singularities, in the pQCD running coupling at µ 2 ∼ Λ 2 QCD ≈ 10 −1 GeV 2 for n f = 3 active flavors. In the vicinity of these singularities our physical predictions are jeopardized. The experience shows that if we do not include some nonperturbative effects, the applicability of this series extends only down to µ ≈ 1 − 2 GeV.
Based on the Lagrangian (3), we can calculate the amplitude and then the 0νββ half-life as [52] T 0νββ Here, Q 2 f ∼ 0.01 GeV 2 is the squared energy of the (spacelike) process of 0νββ decay, G j are the phase space factors (G 1 = G 2 = G 3 , and G 4 = G 5 ) [53], and M j is the Nuclear Matrix element (NME) of the operator O j Eqs. (4a)-(4e) at an effective squared energy Q 2 f , M j = A fin |O j (Q 2 f )|A in . These constant parameters depend on the isotopes we are considering. For the considered isotope 136 Xe, the values of NMEs are given in Table I. The Wilson coefficient C j (Q 2 f ) depends on the typical scale of the 0νββ process, and as we mentioned above, this spacelike scale is quite low and some extension of the usual QCD should be taken into account. In Ref. [50] the authors considered a freezing of the QCD running coupling based on the inclusion of an effective glueball mass M , where M 2 ∈ (0.4, 5.0) GeV 2 . This inclusion was made by the shift Q 2 → Q 2 + M 2 in the one-loop pQCD coupling, cf. Eq. (B28). In Ref. [54] the authors cut the effective theory of QCD at a reasonable scale (see the argumentation above) Q = 2 GeV, and below it down to µ = 0.1 GeV they considered a new effective theory without quarks called Chiral Perturbation Theory (ChiPT) and many-body methods (cf. also [55,56]). Finally, in Ref. [57] the pion mechanism is considered, where the hadronization of quarks and gluons is produced within the effective vertices given by operators (4a)-(4e).
In the present work, we propose an alternative method to deal with this low-energy problem. We propose to extend the applicability of QCD through the Dispersion Relations, which are integrals in the complex Q 2 -plane, which allow us to avoid the appearance of the Landau singularities in a natural way. For details on the construction of such models, we refer to Appendix B.

III. RENORMALIZATION GROUP EQUATIONS WITHIN QCD
The renormalized effective operators (4a)-(4e) are scale independent. Then the Renormalization Group Equation (RGE) will define the anomalous dimension matrixγ in the form where the renormalization constant matrices Z of the effective operators imply that we will have in general some mixing between them. The scale Q 2 is considered to be spacelike, i.e., Q 2 ≡ −q 2 is regarded to be nonnegative. Now, the RGE for the Wilson coefficient follows from the fact that the Lagrangian in (3) is independent of the (spacelike) renormalization scale µ 2 ≡ Q 2 . As a consequence, we obtain the RGE in the matrix form The anomalous dimension matrixγ(Q 2 ) is extracted from the renormalization of the composite operators (4a)-(4e). The corresponding available anomalous dimension factors and matrices are collected in Appendix A: for the operators O 1 -O 3 from Ref. [58] (at the one-loop and two-loop level); for the operators O 4 -O 5 from Refs. [54,59] (at the one-loop level). If we rewrite Eq. (8) in terms of the pQCD running coupling a(Q 2 ) ≡ α s (Q 2 )/π, the RGE can be solved at the twoloop level explicitly, and it is given in the form (for the case of no mixing, i.e., ignoring the problems of diagonalization) where a ≡ a(Q 2 ) and a 0 ≡ a(Q 2 0 ); for β(a) which appears in the renormalization group equation (RGE) Eq. (B1) for the running coupling a(Q 2 ), we took the two-loop truncated form β(a) = −β 0 a 2 (1 + c 1 a). The constants ν and k (1) appearing in Eqs. (9) are where γ (j) (j = 0, 1) are the one-loop and two-loop coefficients, respectively, in the anomalous dimension matrixγ For more details and for different cases of mixing, we refer to Appendix C. We note that in the expansion in Eq. (9b) the terms O(a ν+2 ) are not known if the three-loop anomalous dimension coefficientγ (2) is not known.
In the case of mixing, the analogous formulas for pQCD are obtained in Appendices C 1 and C 2 for the nondegenerate (ν 1 − ν 2 = 1) and degenerate case (ν 1 − ν 2 = 1): cf. Eqs. (C3), (C7), and (C13)-(C14) for the nondegenerate case, and additionally Eq. (C28) for the degenerate case. According to our knowledge, the solution of the two-loop RGE for Wilson coefficients in the degenerate case [which occurs in the n f = 3 regime for the (31) XY mixing of operators O XY 3 and O XY 1 (X = Y )] has not been addressed in the literature hitherto. Within the evolution procedure, the heavy quark thresholds should be taken into account. For this purpose, the evolution matrix U (Q 2 f , Λ 2 LNV ), which connects the "bare" C ≡ C(Λ 2 LNV ) at high momenta with the physical C(Q 2 f ) at Fermi-motion monenta can be written in the following way: where the first equality is given for matching scale of the order of W-boson mass M W = 80.379 GeV [60], and the second equality for large scales, where the theories beyond the standard model play a crucial role. In Eqs. (13), the heavy quark thresholds are at Q t = κm t = 163κ GeV; Q b = κm b = 4.20κ GeV; and Q c = κm c = 1.27κ GeV [60], where we will choose κ = 2 (in general, κ ∼ 1). Note that the variation of the threshold parameter κ is numerically not important in comparison with variation of other parameters. We will use Eq. (13a), i.e., we will take Λ LNV = M W throughout. 4 In the n f = 3 regime, we will use AQCD, because the realistic Fermi motion scale Q 2 f ≈ 0.01 GeV 2 in this regime is quite low and the deviation of the AQCD couplings from the underlying pQCD couplings is significant. In the regimes n f ≥ 4 we use the underlying pQCD in 3δ AQCD because there the A(Q 2 ) coupling practically coincides with the underlying pQCD coupling a(Q 2 ) [Eq. (B5) has N = 5]. In two other cases (2δ AQCD, FAPT) we use at n f ≥ 4 the corresponding AQCD couplings out of convenience (because those coupling are available for all n f ). In the one-loop massive QCD (Massive Perturbation Theory: MPT) we used for n f ≥ 4 the underlying pQCD, for simplicity. 5

IV. EVALUATION OF RGE WITH IR-SAFE COUPLINGS
As mentioned in the Introduction, in AQCD the coupling a(Q 2 ) = α s (Q 2 )/π gets replaced by a coupling A(Q 2 ) where the latter reflects correctly the holomorphic (analytic) behavior of the spacelike QCD physical quantities D(Q 2 ). This means that A(Q 2 ), in contrast to a(Q 2 ), has no Landau singularities in the complex Q 2 -plane, or equivalently, A(Q 2 ) is a holomorphic function for Q 2 ∈ C\(−∞, −M 2 thr ] where M 2 thr is a positive threshold scale, M thr ∼ 0.1 GeV 2 . 6 Here we refer to Appendix B for various AQCD variants. Usually they are constructed with the dispersion relation approach, i.e., starting with a specific form of the discontinuity function ρ A (σ) = ImA(Q 2 = σ exp(−iπ)) for positive σ, and the holomorphic coupling A(Q 2 ) is a dispersion integral involving ρ A (σ), cf. Eq. (B3). Due to the asymptotic freedom, ρ A (σ) at large σ > 1 GeV 2 (practically) coincides with the discontinuity function ρ a (σ) of the underlying pQCD coupling a(Q 2 ) (the latter is defined in a specific chosen renormalization scheme). At low positive σ 1 GeV 2 , the discontinuity function ρ A (σ) is in principle unknown and can be parametrized with Dirac-delta functions, cf. Eqs. (B7) for n = 2, 3 (2δ and 3δ AQCD), Eq. (B29) for one-loop "massive" coupling (MPT). In (Fractional) Analytic Perturbation Theory [(F)APT], the discontinuity function ρ A (σ) is considered to coincide with its pQCD analog ρ a (σ) for all σ values (all the way down to σ = 0), cf. Eq. (B6).
Various AQCD variants [nδ AQCD (n = 2, 3), FAPT, and massive one-loop AQCD (MPT)] are summarized in Appendix B. In the following we will argue that in AQCD the result for the Wilson coefficient C(Q 2 ) is really obtained from the pQCD result (9b) by the replacements (14). We will show this in the case of no mixing, while the extension to the case of mixing of operators is given in Appendix C.
The renormalization group equation (RGE) for a Wilson coefficient C(Q 2 ) as a function of the effective spacelike scale Q 2 in pQCD has the form 7 where γ (n) are the (n + 1-loop) coefficients of the anomalous dimension, and the pQCD expansion of C(Q 2 ) in terms of a(Q 2 ) has the form [cf. Eq. (9b)] where C is a Q 2 -independent quantity. Using this expansion in the RGE (15), and the RGE running of the pQCD coupling a(Q 2 ) according to Eq. (B1), it is straghtforward to see that the index ν and the expansion coefficients k (j) are [cf. Eq. (10) for the two-loop case] etc. The RGE in AQCD is obtained by making analytic the LHS and the RHS of the RGE (15) where the expansion of C(Q 2 ) has the form (16). This is performed with the replacements a ν → A ν as explained in Appendix B One may wonder whether in AQCD the index ν and coefficients k (j) (j = 1, 2, . . .) are the same as in pQCD Eq. (17); they turn out to be the same. Namely, the LHS of the above RGE (i.e., the first line), when using the AQCD relations (B21)-(B22), can be shown be equal to When we equate this expression with the RHS [i.e., the last line in Eq. (18), we obtain for ν and k (j) (j = 1, 2) the same expressions Eqs. (17) as obtained by the pQCD approach.
The conclusion of this exercise is that the solution of the RGE for Wilson coefficients C(Q 2 ) in AQCD is the same as in pQCD, with the replacements a(Q 2 ) ν+m → A ν+m (Q 2 ) in the pQCD expansion (16). Therefore, the relation (9b) in AQCD obtains the form where the above expression U (Q 2 ; Q 2 0 ) (A) is the RGE-evolution matrix in AQCD for the Wilson coefficient from an effective (higher) scale Q 2 0 to an effective (lower) scale Q 2 . In the case of mixing, the analogous formulas for AQCD are obtained in Appendices C 1 and C 2 for the nondegenerate (ν 1 − ν 2 = 1) and degenerate case (ν 1 − ν 2 = 1): cf. Eqs. (C3), (C7), and (C15)-(C16) for the nondegenerate case, and additionally Eqs. (C31)-(C32) for the degenerate case.
In the general approach, applied in nδ AQCD (n = 2, 3) and in one-loop "massive" AQCD (MPT), where the general power analogs A ν are constructed via the generalized logarithmic-derivative analogs A ν+m , Eqs. (B17), it is important to apply the truncations in the evaluation of A ν in Eq. (B17b) consistent with the loop-level in the expression for the Wilson coefficients. When the anomalous dimension γ(a) is known only at one-loop level, then we have and the expression in Eq. (B17b) has only one term On the other hand, when the anomalous dimension is known at the two-loop level, Eq. (20), then the expression in Eq. (B17b) has two terms 8 In practice, this implies that the (two-loop) expression (20) obtains the form where we consistently ignore the terms ∼ A ν+2 . It turns out that with such evaluation we get, even at low |Q 2 | < 1 GeV 2 , reasonable convergence behavor for the Wilson coefficient C(Q 2 ) when going from the one-loop to the two-loop case, see Sec. V. This is probably related with the fact that in AQCD the numerical hierarchy . . is valid in general not just for high |Q 2 | but even for low |Q 2 | 1 GeV 2 .
There is no such hierarchy in pQCD, because of the Landau singularities at or close to |Q 2 | 1 GeV 2 . In the case of mixing, analogous approach is applied for ν = ν 1 and ν = ν 2 , and we refer to Appendix C for more details.
In the case of FAPT, although the use of Eqs. (22)-(23) is an entirely acceptable option, we will follow the more special FAPT-type approach as described in Eqs. (B26)-(B27). In the case of FAPT, this is equivalent to the evaluation of A As mentioned at the end of Sec. III, AQCD will be applied here always in the n f = 3 (low-Q 2 ) regime. In the regimes n f ≥ 4 in general the underlying pQCD approach will be applied; in FAPT and 2δ AQCD, the AQCD approach will be applied also in the regimes n f ≥ 4 for the aformentioned reasons of conveniency.

A. Evolution matrix elements for Wilson coefficients
For evaluation of QCD correction to the 0νββ-decay, we need the physical observable, i.e., the half-life quantity based on OPE. The first question is how the evolution factors or matrices U (Q 2 f ; Λ 2 LNV ) [cf. Eqs. (12) and (13a)] behave when the Fermi motion scale Q 2 f varies downwards towards the realistic values Q 2 f ∼ 0.01 GeV 2 . We will apply a variety of AQCD frameworks: 3δ AQCD [33,61] which has the zero limit in deep IR regime, A(0) = 0; 2δ AQCD [16,17,61] and the one-loop "massive" AQCD (MPT) Eq. (B28), all these having finite positive IR limit A(0) > 0; and FAPT in the MS scheme, cf. Eq. (B27), which gives a nonholomorphic A(Q 2 ) in the point Q 2 = 0 [but has also The MPT coupling (B28) used is taken in two variants: (I) The first one is with M = 1.5 GeV (for n f = 3) and with the scale Λ 3 (i.e., at n f = 3) fixed in such a way that the underlying one-loop pQCD coupling achieves at Q 2 = M 2 Z the value a(M 2 Z ; n f = 5) = 0.1181/π, resulting in Λ 3 = 0.1588 GeV. 9 . This variant will be denoted as MPT(1.5). (II) The second variant is with M = 0.3 GeV (for n f = 3) and Λ 3 = 0.234 GeV, which is suggested by the works of Refs. [62]. 10 This variant will be denoted MPT(0.3). Most of the variants of MPT used in the literature have 0.3 GeV ≤ M ≤ 1.5 GeV (cf. also [50]).
In all other cases, the couplings are normalized in such a way that, at the high scale Q 2 = M 2 Z (and n f = 5) their underlying pQCD coupling (when tranformed to the MS scheme, if needed) achieves the value α s (M 2 Z ; MS) = 0.1181 which is the central value of the present world average [60].
In Fig. 2 we present, for illustration, the four elements of the evolution matrix U (Q 2 ; Λ 2 LNV ) (12) as a function of the Fermi motion scale Q 2 (0 < Q 2 < 5 GeV 2 ), for the case of the mixing of the operators O XX . At relatively large (unrealistic) values Q 2 ≈ 5 GeV 2 , the matrix elements approximately coincide, as they should (asymptotic freedom). However, at more realistic Fermi motion scale values Q 2 < 1 GeV 2 , the 3δ AQCD and MPT(1.5) predictions in general differ significantly, especially for the two-loop (i.e., NLO) anomalous dimension case. The other interesting feature is that the predictions with two-loop (NLO) and one-loop (LO) anomalous dimension do not differ much for 3δ AQCD model, and even less so for MPT(1.5) model. 11 This apparent convergence suggests that the IR-safe versions of QCD (AQCD), once specified, will in general give us definite quantitative predictions for the evolution matrices U (Q 2 ; Λ 2 LNV ) even for very low (realistic) Fermi motion scales Q 2 ∼ 0.01 GeV 2 , starting at least at the two-loop level of the anomalous dimension; and at the one-loop level the predictions can be taken at least as first qualitatively correct estimates. On the other hand, the specific details of the applied AQCD in the deep IR regime can affect quite significantly the values of the evolution matrices (for Q 2 < 1 GeV 2 ); for example, if the AQCD has zero value of the coupling A(Q 2 ) at Q 2 → 0 (such as 3δ AQCD) or a finite nonzero value [such as MPT(1.5), and even more so MPT(0.3) and 2δ AQCD].
In Table II we present the values of the elements of the evolution matrix U (Q 2 f ; Λ 2 LNV ) for Q 2 f = 0.01 GeV 2 , for the various operators (4) or operator mixings, for the described AQCD frameworks. We included in the Table also the values of the coupling πA(Q 2 ) at low scales Q 2 = 0.01 GeV 2 and Q 2 = 0. The results for (one-loop) pQCD are not included, because the pQCD coupling has Landau singularity at Q 2 ≈ 0.025 GeV 2 which is larger than the Fermi motion scale Q 2 f = 0.01 GeV 2 . In Table II we can see that the "strength" of AQCD in the deep infrared regime 9 As explained in the previous Section, the quark thresholds are taken at Q 2 = (κmq) 2 with κ = 2. The threshold condition for the one-loop pQCD coupling is simply continuity. 10 This coupling is parametrized so that it describes in the infrared an effective charge appearing in the DGLAP equation for the parton distribution functions in the pion, and in the ultraviolet it behaves as a (one-loop) pQCD coupling. The underlying pQCD coupling at Q 2 = M 2 Z is then a(M 2 Z ; n f = 5) = 0.1264/π. 11 The values of U (12)22 may at first sight suggest otherwise, but in this case we should keep in mind that all these values are all not far from zero.
, for 3δ AQCD and MPT(1.5), at one-loop (leading order) and two-loop level (next-to-leading order) of the anomalous dimension.
(i.e., the values of coupling A at very low Q 2 ) significantly affect the values of the evolution matrix elements U ij . For example, the results for U ij in 2δ AQCD and in MPT(0.3) are similar, and appear to be influenced largely by the high values of their coupling A(Q 2 ) in the deep IR regime. On the other hand, however, we see that the results for U ij in 3δ AQCD and MPT(1.5) do differ significantly (but not drastically) although the values of A(Q 2 ) in the deep IR regime in both of these frameworks are low. 12 This probably has to do with the fact that 3δ AQCD has significantly more complicated behavior of A(Q 2 ) in the deep IR than MPT(1.5) has. 13 The case of FAPT appears to be intermediate between 3δ and MPT(1.5) on one hand and 2δ and MPT(0.3) on the other hand.
Another interesting aspect which can be inferred from Table II is that all AQCD frameworks give a reasonable convergence of the results when going from the one-loop to the two-loop anomalous dimension case, despite the very low (nonperturbative) Fermi motion scale Q 2 f = 0.01 GeV 2 ; see also Figs. 2-4. This is connected with the holomorphic nature of A(Q 2 ).
In Fig. 5 we present the running couplings πA(Q 2 ) for low positive Q 2 for the mentioned AQCD frameworks. In addition, the large-volume lattice results [63] are included (rescaled to the usual Λ MS -scaling convention) on which the low-Q 2 behaviour of the 3δ AQCD is motivated (see Appendix B for more details). Included is also the one-loop pQCD running which is the underlying coupling for MPT(1.5). All curves are for N f = 3. As mentioned, the high-energy reference strength for the underlying couplings is in all cases except MPT(0.3): α s (M 2 Z ; MS; N f = 5) = 0.1181. The large-volume lattice results are not reliable for Q 2 1 GeV 2 . The couplings have the same scaling convention (Λ MS ), but are in general in different renormalization schemes: lattice and 3δ AQCD in the Lambert MiniMOM (LMM) 12 At first sight, one notable exception are the values of U XY (31)12 which are 0.272 and −0.007, respectively, representing a large relative difference. However, since the reference values of the elements of U matrices are ∼ 1, we see that both values ( 0.272 and −0.007) can be considered close to zero and thus similar. 13 MPT(1.5), due to the high value M = 1.5 GeV, has the coupling A(Q 2 ) almost "frozen" in a wide IR region 0 ≤ Q 2 1 GeV 2 , while 3δ AQCD achieves a maximum at relatively low Q 2 ≈ 0.135 GeV 2 , cf. Fig. 11(a) in Appendix B. scheme; 2δ AQCD in the c 2 = −4.9 Lambert scheme; (F)APT in the MS scheme; cf. Appendix B for more details.
The results of Table II, in conjunction with the curves of Fig. 5, reflect the fact that the conclusions on the QCD effects in 0νββ decay from short-range physics significantly depend on the specific behavior of the AQCD coupling in the deep IR regime. The variation of these QCD efects, when different IR-safe AQCD frameworks are used, can at the moment be regarded as an estimate of the uncertainty of these effects. On the other hand, one could adopt the view that those AQCD frameworks which incorporate more physically motivated information in the IR regime than the others should be preferred in these considerations. 3δ AQCD would then definitely fit into this category, because it takes into account in a natural way the additional information provided by the large-volume lattice calculations    [ [63][64][65][66] in the deep-IR regime (see also Appendix B).
In Table III Table  II). We point out that the naive powers A ν do not treat the nonperturbative contributions correctly, in contrast to the power analogs A ν , as argued in Appendix B. The nMPT(1.5) is close to the approach taken in Ref. [50] where one-loop anomalous dimensions were used. The nMPT(1.5) and nMPT(0.3) are difficult to compare with any of the AQCD frameworks. Further, the values of U ij in nMPT(1.5) almost do not vary when Q 2 f increases from 0.01 GeV 2   Table II. MPT(1. 5) nMPT(1.5)

B. Bounds on Wilson coefficients
The upper bounds on the Wilson coefficients can now be obtained by requiring that the expression on the RHS of Eq. (6) is larger than the lower bound on the half-life T 0ν 1/2 ( 136 Xe) Eq. (1). The RHS of Eq. (6) involves the NMEs (Table I), the space factors G 1 = 2.92 × 10 −14 yr −1 and G 4 = 1.57 × 10 −14 yr −1 [53], and the Wilson coefficients at the Fermi motion scale We recall that we use throughout this work Λ 2 LNV = M 2 W . When using the expansion (25) on the RHS of Eq. (6), the following expression for the half-life in terms of the "bare" Wilson coefficients C k ≡ C k (Λ 2 LNV ) is obtained (cf. also [49]) where We used here the simplified notation ] is zero at one-loop and nonzero at two-loop level of anomalous dimension, cf. Table II. We mention that factor 2 appears in the terms with C XY j (j = 1, 3); this is so because the operators O XY j (j = 1, 3) for XY = LR and RL, Eqs. (4), are symmetric under the interchange of L and R, and hence: [49]). Further, we recall that the values of NMEs are given in Table I.
In the two-loop running of the evolution factors or matrices U ≡ U (Q 2 f ; Λ 2 LNV ), we used the two-loop anomalous dimensionγ of Ref. [58] in the naive dimensional regularization MS (NDR-MS) scheme. On the other hand, the NMEs are often evaluated in a different, Regularization-Independent (RI, also named MOM) scheme. Since the values of NMEs M j have large uncertainties (by about a factor of 2), we used the NDR-MS expressions for the anomalous dimensions, which have the attractive feature of being independent of the gauge-fixing parameter (in contrast to the case of IR scheme).
We can now obtain the upper bounds on the values of the "bare" (new physics) Wilson coefficients |C j | (≡ |C j (Λ 2 LNV )| by assuming that only one operator contributes dominantly to the half-life. This then gives us the upper bounds for various values of the Fermi motion scale Q 2 f = 1.0, 0.1 and 0.01 GeV 2 as given in Tables IV and V, for the cases of one-loop and two-loop anomalous dimension matrices, respectively, for various AQCD frameworks. In these Tables we included, for comparison, the results of pure pQCD approach (only for Q 2 f = 1 GeV 2 and 0.1 GeV 2 ), 15 and for the "bare" case when there are no QCD effects (the evolutions factors or matrices are unity).
In Figs. 6-9 we present the upper bounds as a function of Q 2 f in an extended interval 0.01 GeV 2 ≤ Q 2 f < 5.0 GeV 2 . At (artificially) large values of the (Fermi motion) scales Q 2 ≈ 5 GeV 2 , we can see in these Figures that the upper bounds for various AQCD variants approximately coincide, as it should be due to the asymtotic freedom.  15 We recall that the values Q 2 f < 1 GeV 2 in pQCD become very unreliable or impossible to obtain, due to the Landau singularities.        reduced by this factor of 2. Therefore, in order to discern whether the (A)QCD effects for the upper bound of a Wilson coefficient C j (M 2 W ) are more important than the uncertainties of the values of NMEs, we will consider that where is the corresponding upper bound value when there are no QCD effects (U is unity then). The first inequality in Eq. (28) gives a more stringent upper bound on |C j | (by at least a factor of 2) than when QCD is ignored, and the second one gives a less stringent upper bound (by at least a factor of 2).
With these naive criteriums, we can infer from Table V  . It is interesting that most (but not all) of these qualitative conclusions for the mentioned Wilson coefficients are also valid when regarding the upper bounds obtained with the use of the one-loop anomalous dimensions, cf. Table IV; C XX 1 is here a notable exception. For the coefficients C 4 and C 5 , only the one-loop anomalous dimensions are available. If we regard the upper bounds for them obtained in this way as indicative, then we conclude that in all AQCD variants and for almost all these coefficients the QCD effects are important, where for C 4 the upper bounds become less restrictive and for C 5 more restrictive, cf. Table IV. Some of the upper bounds obtained in the described analysis are quite high, at least in certain specific ranges of values of (Fermi motion scale) Q 2 , as seen in some of Figs. 6-9, and in some cases in Tables IV-V. This is so because in such cases in the corresponding coefficients β XY j Eq. (27), which consist mostly of a sum of two terms of the type  Table II (cf. also Table III).

VI. NONPERTURBATIVE CONTRIBUTIONS IN THE SUB-GEV REGIME
In this work, we used for NMEs the results of the works [2, 48] (cf. also [49]), which were derived by a different method than those of Ref. [54]. On the other hand, in Ref. [54], N N , ππ and πN interactions are used. In addition, in [54] the effective Standard Model with QCD is matched to Chiral Perturbation Theory (ChiPT) already at scales Q = Λ χ ∼ 1 GeV, and ChiPT couplings (g N N j , g ππ j and g πN j ) are then evolved to the Fermi motion scales Q = Q f ∼ 0.1 GeV by nonperturbative renormalization, where they are used in the calculation of NMEs 16 by many-body methods. In a simplified notation, NMEs of Ref. [54] [without the Wilson coefficients at the scale Λ 2 χ , C XY k (Λ 2 χ )], here denoted as M (χ) , can then be represented schematically as M (χ) = g αβ j (Q 2 f )M red , where g αβ j (Q 2 f ) are the RGE-evolved ChiPT couplings (α, β = N, π) and M red are "reduced" NMEs. As argued in [54], the N N ChiPT couplings g N N j with j = 2, 3, 4, 5 have significant effect of the running, i.e., where we took into account that m π ∼ Q f ∼ 0.1 GeV, and g ππ j ∼ 1 GeV 2 as determined by lattice calculations [68]. On the other hand, the estimate g N N j (Λ 2 χ ) ∼ 1 is obtained by the naive dimensional analysis (NDA). 17 The enhancement Eq. (29) comes from contributions of g ππ j couplings, i.e., from pion exchanges. These important effects increase those NMEs which, in the expression for the inverse half-life (T 0νββ 1/2 ) −1 , appear (in our notation) at the Wilson coefficients of the sectors (12) XX and (31) XY (X = Y ). In our approach, where we evolve the Wilson coefficients down to the scales Q f ∼ 0.1 GeV, the ChiPT effects Eq. (29) of the approach of [54] can be reformulated in the (12) XX sector by adding to the evolution matrixÛ (1) (Q 2 ) A , cf. Eq. (C15b), a higher-twist (D = 2) OPE term whereμ 2 is a D = 2 matrix condensate. Such a nonperturbative term ∼μ 2 /Q 2 would appear also in the brackets on the right-hand side of RGE (15). The condensateμ 2 is expected to be the expectation value of certain operators, with quark structure, in the nucleon. The ascertain the explicit structure of this operator would require a special treatment which goes beyond the scope of this work. At Q 2 = Q 2 f ∼ 0.01 GeV 2 , the contribution D = 2 is expected to dominate [cf. Eq. (C15a)] 18 This implies that in this formulation, the Wilson coefficients (C15a) get amplified by a factor of ∼ 10 2 in the running from Q = Λ χ down to Q = Q f . Stated otherwise, the effect of RGE-running (enhancement) of the ChiPT coefficients g N N j in the decay amplitude in the formulation of Ref. [54] is reflected in our formulation by the enhancement of the corresponding Wilson coefficients Eq. (31). Schematically, (T 0νββ It turns out that a rough estimate can be made for the values of D = 2 matrix condensate of the expression (30). Namely, for Q 2 > 1 GeV 2 the leading-twist (D = 0) contribution is assumed to be dominant, and at Q 2 ∼ 1 GeV 2 the two contributions (D = 0, 2) are assumed to start competing (i.e., become comparable). This implieŝ This then gives, for both 3δ AQCD and 2δ AQCD, the values: (μ 2 ) 11 ≈ (μ 2 ) 21 ≈ 2.4-2.8 GeV 2 ; (μ 2 ) 22 ≈ 0.4 GeV 2 ; and (μ 2 ) 12 ∼ 10 −3 GeV 2 . These arguments were made for the sector (12) XX . However, very similar arguments can be made for the sector (31) XY (X = Y ).
As mentioned, the sectors (3) XX and (45) XY are apparently not affected by the possible enhancements mentioned above, i.e., D = 0 contribution appears to be the dominant one in these sectors. This suggests that our results in these sectors, including the upper bounds on the Wilson coefficients C k ≡ C k (Λ 2 LNV ), are more complete in the sense that the values of NMEs in these sectors would eventually become comparable in the various approaches [2, [54][55][56]. 19 Therefore, the results in these sectors may allow us to decide in the future which of the AQCD frameworks is better.
More detailed analysis of the consequences of different methods of the calculation of NMEs for our results goes beyond the scope of the present work.
We wish to add a further comment on the mentioned problem of the large theoretical uncertainties of the NMEs M j = A fin |O j (Q 2 f )|A in . No specific regulator or scheme is involved in the approaches which calculate NMEs. On the other hand, the anomalous dimension matrices we used in RGEs for the Wilson coefficients were in the aforementioned NDR MS scheme [58]. Therefore, this aspect can be regarded as part of the uncertainty of NMEs used in this work.
The present work is presented in such a way that the possible future improved results for NMEs M j = A fin |O j (Q 2 f )|A in at Q f ∼ 0.1 GeV (and preferably in the NDR MS scheme) can be readily used to generate the updated upper bounds on the Wilson coefficients, for various AQCD frameworks, at least in the mentioned sectors 18 It is expected that D > 0 terms are IR-regularized. A plausible scenario is thatμ 2 /Q 2 gets replaced byμ 2 /(Q 2 + M 2 2 ) where M 2 ∼ 0.1 GeV (∼ mπ);μ 4 /Q 4 gets replaced byμ 4 /(Q 2 + M 2 4 ) 2 where M 4 ∼ 0.01 GeV; etc. In such case, D = 2, 4, . . . contributions are all suppressed at Q ∼ 1 GeV; when Q decreases, D = 2 contribution freezes at Q ∼ 0.1 GeV to its maximal value and D = 4 term is still suppressed; D = 4 contribution freezes to its maximal value at Q ∼ 0.01 GeV; etc. In such a scenario, D = 2 term is dominant for Q ∼ 0.1 GeV, if we assume that the values of the dimensionless quantitiesμ 2 /M 2 2 ,μ 4 /M 2 4 , etc., are mutually comparable. 19 The contributions from g ππ j and g πN j are not present in Ref.
[2], but are present in Refs. [54][55][56]. They give contributions to NMEs comparable with those of the ChiPT-enhanced g N N j (Q 2 f ). However, g ππ j and g πN j are not ChiPT-enhanced, and we do not include them in the present discussion.
(3) XX and (45) XY where D = 2 contributions do not appear. Namely, the values of the evolution matrix elements, at Q = Q f = 0.1 GeV, are given in Table II, which means that they can be reused with the future new values of NMEs in Eqs. (27) to obtain β XY j 's, and the new upper bounds for |C j (Λ 2 LNV )| can be obtained by the simple manipulation described in this Section.
We recall that we used in our calculations the values of NMEs of Refs. [2,48], which were obtained in the factorization approximation and using many-body calculations with interactions between nucleons. These quantities correspond to the reduced NMEs M red (at Q ∼ 0.1 GeV) mentioned in this Section. 20 On the other hand, in the discussion in this Section we referred to the NMEs of Ref. [54] which have, due to the intermediate use of ChiPT (χ), the structure M = g αβ j (Q 2 f )M red , and suggested that at least some of the effects of the running of the ChiPT couplings g αβ j (Q 2 ) from Q ∼ 1 GeV down to Q ∼ 0.1 GeV may be contained in the D = 2 termμ 2 /Q 2 of the corresponding Wilson coefficients, Eq. (30). This then raises yet another question, namely how can M and M red (both at Q ∼ 0.1 GeV) be related in principle in our approach which has no ChiPT involved. One possibility is that in our approach the use of M red is made for NMEs (i.e., M → M red , as we did it in this work) and that all other effects of physics from the regime Q > 0.1 GeV are contained in the (OPE) expansions (30); in particular, that the nonperturbative physics from the regime 0.1 GeV < Q < 1 GeV is contained in some sectors dominantly (though not exclusively) in the condensate matrix 21μ 2 . This is, of course, only a conjecture, as these difficult aspects go beyond the scope of the present work.

VII. CONCLUSIONS
In this work we investigated possible QCD effects in 0νββ decays dd → uuee within the scenarios of new LNV physics which are parametrized as short-range dimension-9 operators O j , cf. Eqs. (3)-(4). These QCD effects are reflected in the running of the Wilson coefficients C j of such operators, from the new physics scales Λ 2 LNV (taken here as M 2 W ∼ 10 4 GeV 2 ) to the typical 0νββ-decay (sub-GeV) spacelike scales Q 2 f ∼ 0.01 GeV 2 . For some of these operators their anomalous dimension factors or matrices, which govern the RGE-evolution of the corresponding Wilson coefficients, are known up to two-loops (for O 1 -O 3 , cf. Ref. [58]); for other operators they are known only up to one-loop (for O 4 -O 5 , cf. Ref. [49,54,59]). The pure pQCD treatment of these RGEs is applicable only down to the (spacelike) scales Q 2 ∼ 1 GeV 2 , because below such scales the pQCD coupling a(Q 2 ) (≡ α s (Q 2 )/π) is significantly influenced by the artificial Landau singularities which are situated at 0 < Q 2 < Λ 2 Lan. ∼ 0.1 GeV 2 . In order to achieve the running of the Wilson coefficients C j (Q 2 ) down to Fermi motion scales Q 2 ∼ 0.01 GeV 2 , we employed various variants of QCD where the running coupling A(Q 2 ) [the analog of the pQCD coupling a(Q 2 )] has no such Landau singularities, i.e., various frameworks of AQCD: 3δ, 2δ, MPT(M ) and FAPT. We point out that in such evaluations, in order to evaluate correctly the low-momentum nonperturbative effects, it was important not to treat the analogs of the powers a(Q 2 ) ν as naive powers A(Q 2 ) ν , but rather as A ν (Q 2 ) ( = A(Q 2 ) ν ) which are linear combinations of the (generalized) logarithmic derivatives A ν+m (Q 2 ) (m = 0, 1, . . .). 22 The mentioned evolution of the Wilson coefficients, down to the Fermi motion scales, allowed us then in Sec. V to evaluate the 0νββ half-life of 136 Xe in terms of these coefficients C j (Q 2 f ) and of the corresponding nuclear matrix elements (NMEs) of the operators O j . Comparison of this expression with the presently available lower bound on the mentioned half-life then allowed us to extract the upper bounds for the Wilson coefficients C j (Λ 2 LNV ) at the new physics scale.
Our main conclusions are the following. The values of the evolution factors or matrices U (Q 2 f , Λ 2 LNV ) of Wilson coefficients, when the two-loop anomalous dimensions were used, were in all AQCD frameworks not far from (and often close to) the values obtained for U (Q 2 , Λ 2 LNV ) when one-loop anomalous dimensions were used. This conclusion holds even when the values of the Fermi motion scales are realistic, i.e., very low, Q 2 f ∼ 0.01 GeV 2 . As a consequence, similar conclusion can be made for the extracted values of the upper bounds of |C j (Λ 2 LNV )|. Further, as could be expected, the numerical results for different AQCD frameworks depend largely on the behavior of the coupling A(Q 2 ) in the IR regime Q 2 0.1 GeV 2 . Therefore, for example, the results of 2δ and MPT(0.3) AQCD variants were mutually comparable. The results of 3δ AQCD are not easily comparable with those of other AQCD frameworks, principally because the coupling A(Q 2 ) in 3δ AQCD goes to zero in the deep IR-regime (as suggested by large-volume lattice results). Yet another conclusion of this work is that the described QCD effects are important (more than the present uncertainty of the NMEs) in most of the cases of the considered Wilson coefficients: these effects affect in 20 They do include, though, specific factors which can be interpreted as g N N j ∼ 1, cf. Eq. (F.2) of Ref. [54]. 21 The sectors (12) XX and (31) XY (X = Y ) are expected to give two different condensate matricesμ 2 . 22 Appendix B is a summary of various AQCD frameworks and of the evaluation of Aν , all this information being available in the literature. such cases the upper bounds for |C j (Q 2 )| (when Q 2 = 0.01 GeV 2 ) by more than a factor of two.
Finally, we point out that at present there is a major uncertainty in the NMEs associated with Wilson coefficients of the operator sectors O XX . Namely, according to Ref. [54] such NMEs may involve an enhancement effect of ∼ 10 2 in comparison to other approaches [2, 48,53,55,56]. This effect may require, in our approach in the mentioned operator sectors, introduction of additional, D = 2 contributions ∼ 1/Q 2 in the evolution matrices.

Acknowledgments
This work was supported in part by the Chilean FONDECYT Regular Grants No. 1200189 (C.A.) and No. 1180344 (G.C.). The work of L.G. was supported by CONICYT Chile Grant No. 21160645 and DGIIP of UTFSM.
We are grateful to M. González and S.G. Kovalenko for helpful discussions.

Appendix A: Anomalous dimension at LO and NLO
In this Section we write down the anomalous dimension in the full one-loop approximation and in the currently known two-loop one. The results will be expressed in terms of the number of colors N , and the number of active flavors n f . First, we write down the result for the mixing of operators (4a) and (4b), as obtained in Ref. [58] γ (0),XX (12) We note that the off-diagonal elements here in γ (0) (12) and γ (1) (12) have the opposite sign to those of [58]; this is so because in Eq. (4b) we use the convention σ µν = (i/2)[γ µ , γ ν ], while in Ref. [58] the convention σ µν = (1/2)[γ µ , γ ν ] is used.
Then, we write down the result for the mixing of the operators (4c) and (4a), as obtained in Ref. [58] γ (0),XY (31) The final result that is known up to two-loop approximation corresponds to the operator (4c) [58] γ (0),XX (3) For the operators (4d) and (4e), only the one-loop anomalous dimension is known [59] γ (0),XX This same result can be deduced also from the results of [54], when taking into account the following relations between the operators O µ 6 -O µ 9 of [54] For the operators O µ 6 -O µ 9 of [54] the above relations turn out to be the same, but with L ↔ R on the RHS. These relations can be obtained by the use of Fierz transformations and of the color rearrangement identity involving the SU (N ) color generators

Appendix B: IR-safe couplings
This Appendix is a compendium and synthesis of several results obtained and described in various works on IR-safe holomorphic couplings, among them Refs. [9,14,17,33,36,47] and [69] (App. B there).
The pQCD running coupling a(Q 2 ) ≡ α s (Q 2 )/π is defined as a function of the squared momentum Q 2 ≡ −q 2 in the generalized spacelike region, where q 2 = (q 0 ) 2 − q 2 and q represents a typical momentum of a considered process. When q 2 < 0 (Q 2 > 0), the momentum q is considered to be spacelike in the restricted sense (e.g., appearing in deep inelastic scattering and other t-channel quantities, and in current correlators). When q 2 = s > 0 (Q 2 = −s < 0), the momentum is usually called timelike (e.g., appearing in the s-channel type decay widths and cross sections). The generalized spacelike (Euclidean) region of Q 2 is considered to be the entire complex plane with the exception of the timelike semiaxis: Q 2 ∈ C\(−∞, 0], and it is in this region that the running coupling a(Q 2 ) is considered. The running coupling in this region is a solution of the (perturbative) RGE = −β 0 a(Q 2 ) 2 1 + c 1 a(Q 2 ) + c 2 a(Q 2 ) 2 + . . . .
Here, the first two β-coefficients, β 0 = (1/4)(11 − 2N f /3) and β 1 = (1/16)(102 − 38N f /3), are scheme independent in mass independent schemes. The coefficients c j = β j /β 0 (j ≥ 2) characterize the pQCD renormalization scheme [70]. This means that the form of the function β(a; c 2 , c 3 , . . .) represents a definition of the renormalization scheme. Here, the momentum scale Λ QCD will not be regarded as a scheme parameter, but as the momentum (re)scaling definition. This scaling change can be described equivalently as a change of the renormalization scale. Throughout this work, the MS scaling definition (Λ 2 QCD = Λ 2 ) is adopted.
The pQCD coupling a(Q 2 ), which is the solution of the perturbative RGE in a given or chosen renormalization scheme, usually has singularities on the positive axis in the Q 2 -complex plane, 0 ≤ Q 2 Λ 2 QCD (∼ 0.01-1 GeV 2 ); this is in addition to the expected singularities on the negative Q 2 -axis. However, the spacelike QCD observables D(Q 2 ) (such as current correlators, t-channel process quantities, and nucleon structure functions and their sum rules) are are holomorphic (analytic) functions in the Q 2 -complex plane with the exception of a part of the negative semiaxis, Q 2 ∈ C\(−∞, −M 2 thr ] (with a threshold mass M thr ∼ 0.1 GeV). This follows from the general principles of Quantum Field Theories [5,6]. Stated otherwise, spacelike QCD observables D(Q 2 ) are holomorphic functions in the entire (generalized) spacelike region.
These properties of D(Q 2 ) are not reflected qualitatively in the pQCD coupling a(Q 2 ), because the latter has the mentioned Landau singularities (cut and branching points) on the positive axis. This behavior of a(Q 2 ) is unfortunate in the following sense: if in the evaluation of D(Q 2 ) at low |Q 2 | 1 GeV 2 we use the coupling a(Q 2 ) [or a(µ 2 ) with µ 2 = κQ 2 ∼ Q 2 ], we obtain either useless or unreliable results. These Landau singularities (also called Landau ghosts) are usually a cut on an interval 0 ≤ Q 2 ≤ Q 2 br , and the point Q 2 = Q 2 br (∼ 0.01-1 GeV 2 ) is called the Landau branching point. When we apply the Cauchy theorem to the integrand a(Q 2 )/(Q 2 − Q 2 ) in the Q 2 -complex plane along the path in Fig. 10(a), this gives the following representation of the pQCD coupling a(Q 2 ) in the form of a dispersion integral: Here, ρ a (σ) = Ima(Q 2 = −σ − i ) is the discontinuity (or spectral) function of the coupling a along its cut.
On the other hand, the holomorphic (in the spacelike region) coupling A(Q 2 ) [the analog of a(Q 2 )] has its cut only along the negative semiaxis −∞ < Q 2 < −M 2 thr , and hence the form of its dispersion integral is [cf. Fig. 10(b)] where ρ A (σ) = ImA(Q 2 = −σ − i ) is the discontinuity of A(Q 2 ) along its cut in the complex Q 2 -plane. This coupling has the cut threshold σ min (≡ M 2 thr ) ≥ 0. In contrast to a(Q 2 ), the couplings A(Q 2 ) represent qualitatively correctly the holomorphic properties of the QCD spacelike observables D(Q 2 ), and can thus be regarded as better suited for the evaluation of such quantities. However, A(Q 2 ) have to fulfill various physically-motivated requirements: (a) at high |Q 2 | > 1 GeV 2 they must reproduce the perturbative QCD; (b) at intermediate |Q 2 | ∼ 1 GeV 2 they must reproduce the corresponding QCD phenomenology, especially the physics of the τ lepton semihadronic decays which is well measured; (c) and at very low |Q 2 | < 1 GeV 2 we may require that they have the behavior as suggested by large-volume lattice results for the (Landau gauge) gluon and ghost dressing functions, if the running coupling there is defined in a natural way as a product of these dressing functions.
The high-momentum condition (a) can be also formulated in the following way: in a chosen renormalization scheme (i.e., for a chosen set of values of the scheme c j coefficients, j ≥ 2), the discontinuity function ρ A (σ) coincides at large σ with the pQCD discontinuity function of the underlying pQCD coupling a where M 2 0 can be called the pQCD-onset scale. Then at large |Q 2 | > 1 GeV 2 the requirement that the two running couplings practically coincide can be written as where Λ 2 ∼ 0.1 GeV 2 and index N must be relatively large, e.g. N = 5. The simplest holomorphic coupling (APT) [7] was constructed from the underlying pQCD coupling by equating ρ A = ρ a for all σ ≥ 0 (and necessarily setting equal to zero the Landau cut discontinuities ρ A (σ) at σ < 0) On the other hand, we constructed two types of couplings which fulfill the condition (a) [i.e., Eq. (B5) with N = 5] and (b) [16,17,33], one type of coupling fulfilling also the deep infrared condition (c) [33]. The discontinuity functions for these two types of couplings are parametrized in the unknown low-σ region (σ < M 2 0 ) by a combination of Dirac-delta functions ρ (nδ) where we expect 0 < M 2 1 < . . . < M 2 n < M 2 0 , and M 2 0 ∼ 1 GeV 2 is the pQCD-onset scale. The corresponding coupling is The Dirac delta functions in the spectral function give a nonperturbative contribution ∆A IR (Q 2 ), in the form of a linear combination of simple fractions ∼ 1/(Q 2 + M 2 j ), and this can be rewritten as a near diagonal Padé approximant ∆A IR (Q 2 ) = [n/n − 1](Q 2 ). Padé approximants [n/n − 1](Q 2 ) are known to approximate the holomorphic functions in the Q 2 -complex plane increasingly well when n increases [71].
On the other hand, the fifth requirement comes from physics at moderate momenta |Q 2 | ∼ m 2 τ (∼ 1 GeV 2 ): the requirement that the physics of the semihadronic τ lepton decays be reproduced, i.e., that the calculated (massless and strangeless) τ decay ratio r (D=0) τ gives the correct well-measured value. This is sufficient for 2δ AQCD [16,17] which has five parameters.
In 3δ AQCD [33], which has seven parameters (n = 3), two additional requirements were imposed, namely that A (3δ) (Q 2 ) ∼ Q 2 when Q 2 goes to zero, and that A (3δ) (Q 2 ) has for positive Q 2 a local maximum at Q 2 ≈ 0.135 GeV 2 , in the Lambert MiniMOM (LMM) scheme. These two requirements are suggested by the large volume lattice calculations [63] for N f = 0 24 of the Landau gauge dressing functions Z gl (Q 2 ) and Z gh (Q 2 ) of the gluon and ghost propagators in the MiniMOM (MM) scheme [72], 25 where the lattice coupling A latt. was defined naturally as: There is yet another, the (2n + 2)'th "hidden" parameter, involved in A (nδ) (Q 2 ): it is the strength of the underlying pQCD coupling a(Q 2 ) (at N f = 3); it can be characterized by the value of α s (M 2 Z ; MS) (at N f = 5) which we take in this work in general as α s (M 2 Z ; MS) = 0.1181 [60]. When a specific coupling A(Q 2 ) has been constructed (a → A), the couplings A n (Q 2 ) which are the analogs of the powers a(Q 2 ) n of the underlying pQCD coupling (a n → A n ), can be obtained in general holomorphic AQCD following the steps presented in Ref. [14] for integer n, and in Ref. [47] for general (noninteger) n. In this construction of A n (Q 2 ) from A(Q 2 ) for integer n, the logarithmic derivatives of A(Q 2 ) play a central role.
Here, first the construction given in Ref. [14] for integer n will be outlined. From the linearity of the "analytization" a(Q 2 ) → A(Q 2 ) follows that the logarithmic derivatives a n+1 (Q 2 ) of a(Q 2 ) 23 The schemes are perturbatively defined, by the (perturbative) β(a) function of the underlying pQCD coupling a. 24 Similar results were obtained also by another group [64], for N f = 0. Further, similar results, but in general with lower statistics, were obtained for N f = 2 [65] and N f = 4 [66]. 25 Lambert MiniMOM (LMM) scheme is the MiniMOM (MM) scheme (where the lattice calculations are performed), but with momenta rescaled to the usual MS scaling: are replaced (i.e., "analytized") in AQCD by the analogous logarithmic derivatives A n+1 (Q 2 ) of A(Q 2 ) where the expression (B10b) is obtained by using the definition (B10a) and the dispersion integral (B3). This construction allows us to evaluate such (truncated) AQCD series whose perturbation series starts with an integer power of a(Q 2 ), e.g., with a(Q 2 ) 1 . Namely, the (leading-twist part of the) spacelike observable D(Q 2 ) has in such a case the power expansion where κ ≡ µ 2 /Q 2 is the renormalization scale parameter (0 < κ ∼ 1). This series can be reorganized in a straightforward way as a series in the logarithmic derivatives (B9) instead 26 The resulting (truncated) series is then evaluated with the A-coupling A weak renormalization scale dependence (κ-dependence) appears here due to the truncation effect. The analog of this expression in pQCD is the series Eq. (B12) truncated at d N −1 (κ) a N (κQ 2 ). The truncated series (B13) differs from the full sum D(Q 2 ) formally by a term ∼ A N +1 (∼ a N +1 ∼ a N +1 ); this is suppressed in comparison to ∼ A N , because AQCD frameworks in general fulfill the hierarchy |A(Q 2 )| > | A 2 (Q 2 )| > | A 3 (Q 2 )| > . . ., for all (non-timelike) scales Q 2 (cf. also Figs. 11), which appears as a consequence of the holomorphic behavior of A(Q 2 ). The truncated series (B13) can be reorganized (rewritten) explicitly in terms of the coefficients d n (κ) of the original perturbation (power) series (B11) where the power analog A n+1 (the A-coupling analog of the power a n+1 ) is a specific linear combination of the logarithmic derivatives A n+m in complete analogy with the pQCD relations Here, the sums are truncated consistently at A N ; we note that A N = A N in this case. We point out that the truncated series (B14) is equal to (B13), and its pQCD analog are the series (B12) and (B11), both truncated at n + 1 = N . Since A(Q 2 ) has in general some nonperturbative contributions in comparison to its underlying pQCD coupling a(Q 2 ), we have A n (Q 2 ) = A(Q 2 ) n (n ≥ 2), and this holds even if the truncation index N in the relations (B15) is very high. On the other hand, at high |Q 2 | > 1 GeV 2 we have in general A n (Q 2 ) ≈ A n (Q 2 ) ≈ A(Q 2 ) n ≈ a(Q 2 ) n , because (2δ and 3δ) AQCD in the high-momentum regime practically coincides with the underlying pQCD, due to the relation (B5) with N = 5 there. If in the series (B14) the naive powers A(Q 2 ) n were used instead of A n (Q 2 ), this would give to the series spurious uncontrollable nonperturbative contributions at low |Q 2 | 1 GeV 2 [73]. Hence it is important to employ the series in logarithmic derivatives instead, i.e., Eq. (B13) [⇔ Eq. (B14)]. In Figs. 11(a),(b) we present the couplings A(Q 2 ), A 2 (Q 2 ), as a function of Q 2 > 0, for the considered 2δ and 3δ AQCD, respectively, and we included the corresponding underlying pQCD coupling a(Q 2 ) and the MS couplinḡ a(Q 2 ) (all are for N f = 3). The "strength" reference value for the pQCD couplings is α s (M 2 Z ; MS; N f = 5) = 0.1181. The coupling A 2 (Q 2 ) is generated by Eq. (B15) with three terms (i.e., N = 4 and n = 1). The naive power A(Q 2 ) 2 is also presented in these Figures; clearly: A 2 (Q 2 ) ≈ A(Q 2 ) 2 at low Q 2 . The pQCD coupling a(Q 2 ) in the LMM scheme has the branching point at Q 2 br ≈ 1.29 GeV 2 (not a pole). In the Lambert scheme with c 2 = −4.9, where 2δ AQCD was constructed, a(Q 2 ) has Q 2 br ≈ 0.066 GeV 2 (a pole), and in the MS scheme Q 2 br ≈ 0.36 GeV 2 (a pole). All these curves were obtained by using the programs [61], written in Mathematica, for the evaluation of the couplings. Until now we have described the case when ν = n + 1 in A n+1 and A n+1 is an integer. However, in many cases in physics, the physical (spacelike) quantities F(Q 2 ), such as here considered Wilson coefficients, have perturbation expansion in powers of a ν = a ν0+n where ν 0 (> −1) is not integer (and n = 0, 1, 2, . . .) In such cases, the results for integer ν = 1 + n can be analytically continued to ν = ν 0 + n [47], i.e., we obtain (a ν (Q 2 )) an. = A ν (Q 2 ), (B17a) where ν = ν 0 + n (n = 0, 1, . . . , N − 1; −1 < ν 0 ); and ρ a (σ) (1 ) is the discontinuity of the one-loop pQCD coupling and the explicit expressions for the coefficients k m (ν) appearing in the relation (B17b) are given in Ref. [47]. The unsubtracted part of the dispersive integral in Eq. (B17c) was obtained by simple continuation of the expression (B10b) to noninteger values (n + 1 → ν). The full dispersive integral in Eq. (B17c) converges in an extended regime of indices ν, namely ν > −1 (not just for: ν > 0). This is so because the basic (unsubtracted) dispersion integral was modified by subtracting and adding the one-loop (F)APT expression A (FAPT,1 ) ν = A (FAPT,1 ) ν which is known explicitly [10] (when ν > 0, this subtraction and addition are not needed) Here, the scale Λ 2 ∼ 0.1 GeV 2 is arbitrary and it appears also in the (one-loop) pQCD discontinuity function ρ a (σ) (1 ) .
The expressions A ν (Q 2 ), which are extensions of the logarithmic derivatives (B10) to noninteger n + 1 → ν, were shown to satisfy the recursive relations (B20) Furthermore, using the explicit expressions for the coefficients k m (ν) (m = 1, 2, 3, 4) obtained in Ref. [47], we can check that the following RGE-type relations hold for A ν : This turns out to be in complete analogy with the RGE in pQCD for the power a(Q 2 ) ν representing thus a cross-check of consistency of our construction of the power analogs A ν (Q 2 ) [AQCD analogs of the powers a(Q 2 ) ν+1 ].
The series (B16) in (IR-safe) AQCD, and its truncated version F [N ] , are then obtained by the simple replacements (B17a) AQCD + O( A ν0+N ). In this context, we point out that the AQCD frameworks in general fulfill the hierarchies | A ν0 (Q 2 )| > | A ν0+1 (Q 2 )| > | A ν0+2 (Q 2 )| > . . ., for all (non-timelike) scales Q 2 , a property which appears to be a consequence of the holomorphic behavior of these quantities (and of A(Q 2 )). The coefficient f n is a linear combination of the coefficients f n , f n−1 , . . . due to the relations (B17b).
Sometimes, as in the degenerate case Appendix C 2, in the perturbation expansion of physical observables we have the mixed powers a ν ln k a (where k = 1, 2, . . .), and they get analytized by the analogous approach [35] a(Q 2 ) ν ln k a(Q 2 ) an.
It is important to point out that the construction of the analytic analogs A ν (Q 2 ) of the powers a(Q 2 ) ν [cf. Eqs. (B10) and (B15) for integer n, and Eqs. (B17) for general n = ν − 1] is an operation which is linear in the (holomorphic) coupling A(Q 2 ), in contrast to the naive construction (A(Q 2 )) ν . This means that, when A → λA. we have: ρ A → λρ A , A ν → λ A ν and A ν → λA ν . Furthermore, in the case of integer n it is clear from the definition (B9) of the pQCD quantity a n+1 (Q 2 ) that its analytic version should be A n+1 (Q 2 ) of Eq. (B10), because the transition from pQCD to AQCD produces only the changes a(Q 2 ) → A(Q 2 ) and a(Q 2 + ∆Q 2 ) → A(Q 2 + ∆Q 2 ). More explicitly and for higher n analogously. One of the consequences of this construction is that A ν (Q 2 ) = (A(Q 2 )) ν . The construction of the AQCD analogs A ν (Q 2 ) of powers a(Q 2 ) ν described here can be applied in any AQCD. On the other hand, the case of APT Eq. (B6), where the discontinuity function ρ A (σ) is in its entirety (i.e., for all σ > 0) the pQCD discontinuity function ρ a (σ), exceptionally allows for a more direct evaluation of A ν (Q 2 ), namely as The extension of the convergence of this integral to the regime −1 < ν can be achieved by subtracting the one-loop (F)APT expression (B19) in the form of dispersive integral and adding it in its explicit form (B19) It turns out that this gives the same result as the aforedescribed general method of construction of A ν (Q 2 ) when applied to the APT case ρ A (σ) ≡ ρ a (σ) if the truncation index in the sum on the RHS of Eq. (B17b) is sufficiently high. We will apply the expression (B27) in the case of FAPT, using the four-loop MS pQCD coupling as the underlying coupling in the form given in Ref. [17] [Eq. (6) there]. In this context, we point out that the approach (B27) can be applied only in the case of the specific AQCD, namely FAPT (i.e., in the case where ρ A (σ) ≡ ρ a (σ) for all σ > 0), while the approach (B17) can be applied in any AQCD. Yet another, rather popular, AQCD coupling, i.e., coupling without Landau singularities, is the "massive" one-loop coupling (MPT) where M 2 ∼ 1 GeV 2 and Λ 2 ∼ 0.1 GeV 2 . The corresponding discontinuity function is In FAPT and in MPT, the deviation from the underlying pQCD at high |Q 2 | remains strong because it has N = 1 in the relation (B5). On the other hand, 2δ and 3δ AQCD described before have N = 5, i.e., they practically coincide with the underlying pQCD at high |Q 2 | > 1 GeV 2 .
Here we summarized the evaluation of the spacelike physical quantities D(Q 2 ). The timelike physical quantities can, in principle, be expressed as contour integrals involving the corresponding spacelike quantities, and can thus also be evaluated in AQCD (for example, cf. [33]).
(C1) can be rewritten as where a ≡ a(Q 2 ). LetV (0) be the "rotation" matrix which diagonalizes the one-loop matrixγ (0)T whereν is, by this definition, a diagonal matrixν When defining the RGE (C2) can be rewritten in the form where the matrixk (1) incorporates the two-loop effectŝ We recall that c 1 = β 1 /β 0 is the (universal) two-loop beta coefficient, cf. Eqs. (B1). Since the first matrix on the RHS of Eq. (C7) is in general nondiagonal, an additional, two-loop, "rotation" is needed to obtain fully decoupled system. This is achieved by a matrixĴ (1) which acts in the following way: such that the RGE for C (1) (a) is a decoupled system i.e., the total expression in brackets on the RHS of Eq. (C9) is a diagonal matrixk (1) D . This can be achieved by the following matrixĴ (1) :Ĵ (C10b) The decoupled system of RGEs (C9) can then be integrated, resulting in where C is a two-component (column) vector independent of Q 2 scale, and a(Q 2 )ν is a diagonal matrix according to Eq. (C4) a(Q 2 ) ν = exp[ν ln a(Q 2 )] = a(Q 2 ) ν1 0 0 a(Q 2 ) ν2 . (C12) Using the relation (C8), the solution for the original vector C(Q 2 ) of Wilson coefficients is where the matrixÛ (1) (a) isÛ (1) (a) = aν + k (1) D + where a ≡ a(Q 2 ), and the other parameters are Q 2 -independent. According to conclusions presented in Appendix B, in AQCD the same relations are valid, but under the consistent replacements a(Q 2 ) ν+m → A ν+m (Q 2 ) The evaluation of A νj and A νj +1 (j = 1, 2) in terms of A νj and A νj +1 is performed along the same lines as explained in Sec. IV Eqs. (23)- (24), but now separately for ν 1 and ν 2 : A νj = A νj + k 1 (ν j ) A νj +1 , and A νj +1 = A νj +1 (j = 1, 2). We recall that in Eq. (C15a) the vector C is Q 2 -independent. This allows us, equally as in pQCD [Eq. (C13)], to rewrite the solution in terms of the initial condition valuesÛ (1) (Q 2 0 ) (A) The matrixÛ (Q 2 ; Q 2 0 ) (A) is the two-loop (RGE-)evolution matrix for the Wilson coefficients C from a (higher) scale Q 2 0 to a (lower) scale Q 2 , in the case of (2 × 2) mixing, in AQCD with IR-safe and holomorphic coupling A(Q 2 ).

RGE for Wilson coefficients with mixing -degenerate case
In some exceptional cases, the eigenvalues of the matrixν [Eq. (C4)] can satisfy the relation (X = Y ) in Eqs. (4) for n f = 3, where the anomalous dimension matrix is known at the two-loop level. We recall that AQCD should be applied in the n f = 3 regime.
In such a case, the two-loop matrixĴ (1) Eq. (C10a), which is needed for the decoupling of the two RGEs, does not exist because one term there has zero in the denominator. In such a case, we have to proceed in a modified way. At the two-loop level, the matrixĴ (1) now has the limited form Here as earlier, a ≡ a(Q 2 ), and C 2 is a Q 2 -independent constant. The terms O(a ν2+2 ) are not specified in Eq. (C20b) because they are affected by the (unknown) three-loop contributions. Inserting the solution (C20) into the first RGE (C19a) gives us a nonhomogeneous differential equation for C where a ≡ a(Q 2 ) and a 0 ≡ a(Q 2 0 ). The first term on the RHS of Eq. (C22) represents a solution to the homogeneous version of Eq. (C21), and the second term a particular solution to the full (nonhomogeneous) Eq. (C21);Ũ 1 (a) is the evolution functionŨ 1 (a) = a ν1 +k The solution (C22) implies that the expression is a Q 2 -independent constant (C 1 ). This, and the relation (C20a), imply that the solution for C (1) (a) can be written in the form Using this solution, we can "rotate" back to the original basis of the Wilson coefficients using the relation (C8) and the explicit form (C18) ofĴ (1) in the considered degenerate case. In analogy with the algebra performed in the previous Subsection C 1, we obtain now C(a) = V (0)Û (1) (a) C, where C T = (C 1 , C 2 ) is the vector with the two Q 2 -independent constants, and the matrixÛ (1) (a) is now (the considered degenerate case ν 1 − ν 2 = 1) (1) 11 a ν1+1 ,k 12 a ν1 ln a +k (1) 21 a ν1+1 , 1 2k (1) 21k (1) 12 a ν1+1 ln a + (a ν2 +k where terms of higher order, which are affected by (unknown) three-loop contributions, were neglected. We recall that a ≡ a(Q 2 ). It can be shown that the result of the nondegenerate case considered in the previous Appendix C 1, Eq. (C14b), is in the case of ν 1 − ν 2 = 1 − ε (with → 0) the limiting case of the above result Eq. (C28), as it should be. Namely, when ν 1 − ν 2 = 1 − , we have a 1−ν1+ν2 1 − ν 1 + ν 2 = a = 1 + ln a + (−k As in the previous Subsection C 1, the transition to the AQCD is obtained by the replacements a ν+m → A ν+m and by a ν ln a [≡ (d/dν)a ν ] → (d/dν)A ν in Eqs. (C27)-(C28) where the matrixÛ (1) (Q 2 ) (A) in AQCD iŝ (1) 21 A ν1+1 (Q 2 ), 1 2k (1) 21k (1) 12 d dν A ν (Q 2 )| ν=ν1+1 + (A ν2 (Q 2 ) +k (1) 22 A ν2+1 (Q 2 ))   .
As in the previous Appendix C 1, the relation (C30) can be written in the form where the matrixÛ is the evolution matrix for the Wilson coefficients from the (upper) scale Q 2 0 to the (lower) scale Q 2 . Analogously as in the nondegenerate case in Appendix C 1 [and in Sec. IV Eqs. (23)- (24) in the case of no mixing], the evaluation of A νj and A νj +1 (j = 1, 2) in Eq. (C31) is performed in terms of A νj and A νj +1 as follows: A νj = A νj + k 1 (ν j ) A νj +1 , and A νj +1 = A νj +1 (j = 1, 2).