Two- and three-loop anomalous dimensions of Weinberg's dimension-six CP-odd gluonic operator

We apply, for the first time, a fully automated extension of the $R^*$-operation capable of calculating higher-loop anomalous dimensions of n-point Green's functions of arbitrary, possibly non-renormalisable, local Quantum Field Theories. We focus on the specific case of the CP-violating Weinberg operator of the Standard Model Effective Field Theory whose anomalous dimension is so far known only at one loop. To avoid mixing with gauge-variant operators we employ the background field method. Based on both the 't Hooft-Veltman scheme and the Larin scheme we develop two different strategies to consistently handle the Levi-Civita symbol in the $R^*$-operation. With this setup, we calculate the two- and three-loop Weinberg anomalous dimension in Yang-Mills theory. We find sizeable two-loop corrections and large three-loop corrections, because of the appearance of a new quartic group invariant. We discuss phenomenological implications for measurements of the electric dipole moments of nucleons and future applications of the method.


INTRODUCTION
The absence of evidence for beyond-the-Standard Model (BSM) physics in high-energy proton-proton collisions at the LHC, and in large classes of low-energy precision measurements all indicate that the scale of BSM physics (Λ) is significantly higher than the electroweak scale (v): Λ ≫ v ≃ 246 GeV. In such a scenario, the effects of BSM physics at low energies E ≪ Λ, can be described in terms of effective operators consisting of SM fields that obey the Standard Model (SM) gauge and Lorentz symmetries [1][2][3]. The resulting framework is called the SM effective field theory (SMEFT). The SMEFT Lagrangian contains an infinite number of operators that can be ordered by their dimension. Effects of higher-dimensional operators on low-energy observables are suppressed by additional powers of E/Λ.
The connection between observables at (relatively) low energies and the SMEFT operators at the scale Λ, where the effective operators can be matched to specific UVcomplete BSM models, is determined by renormalisationgroup equations (RGEs). The RGEs depend on anomalous dimensions that can be calculated in perturbation theory in an expansion in small coupling constants. The complete one-loop anomalous dimension matrix of dimension-six SMEFT operators has been obtained [4][5][6]. Already for dimension-six, the number of operators is large and the general mixing structure of the RGEs is rather complex. It has been observed that the one-loop anomalous dimension matrix is almost holomorphic [7], but it is not clear whether this feature extends to higher order. Higher-order anomalous dimensions have been cal-culated for subsets of dimension-six operators [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23], but due to the hard technical nature of the calculations the complete matrix is not known. Higher-order anomalous dimensions can be used to (1) improve the precision of SMEFT contributions to LHC processes or low-energy precision observables, (2) to study the convergence of the perturbative expansions, and (3) to investigate the structure of the SMEFT mixing pattern.
In this letter we extend the R * -operation to the framework of the Standard Model Effective Field Theory (SMEFT), and develop an efficient and highly automated method to calculate higher-order QCD anomalous dimensions of SMEFT operators. The R * -operation provides a way to subtract UV and IR divergences from Euclidean Feynman diagrams, taking care of the combinatorics of overlapping divergences [24][25][26]. Recently, it has been extended to Feynman diagrams with arbitrary numerator structure [27]. So far the R * -operation has been used extensively in calculations of anomalous dimensions in QCD; see, e.g., [28][29][30][31]. However, the R * -method is not limited to pure QCD and can be applied to arbitrary local quantum field theories.
To avoid the complicated mixing structure of SMEFT operators, we focus on a specific SMEFT dimension-six operator. Once tested and developed, the method can be extended to a larger set of operators without too many additional complications. The operator we now consider is the CP-violating gluonic operator, often called "the Weinberg operator" [32], defined as in terms of the gluon field strength G a αβ , the Levi-Civita tensor (LCT) ǫ µναβ , and the gauge group structure constants f abc . The Weinberg operator is induced in various classes of BSM models with additional CPviolating phases [33][34][35], and leads to nonzero electric dipole moments (EDMs) of nucleons, nuclei, and diamagnetic atoms (such as 199 Hg [36]). Current experimental EDM limits [37] set strong constraints on BSM models that induce the Weinberg operator. The one-loop anomalous dimension of the Weinberg operator was obtained in the original paper by Weinberg [32], albeit with the wrong sign, finding sizeable QCD corrections from the evolution of C W from highto low-energy scales. The calculation was corrected in Refs. [38][39][40] that also calculated the mixing, proportional to the small quark masses, of the Weinberg operator into the quark chromo-electric dipole moment. The only other operator C W can mix with is the QCD theta term ∼ ǫ µναβ G a αβ G a µν , but this mixing is of little phenomenological use as the bare theta term is an unknown SM parameter. Furthermore, the renormalised theta term vanishes after a Peccei-Quinn mechanism [41]. We will not consider the mixing into the theta term in this letter.
The potential phenomenological implications of the unknown higher-loop anomalous dimensions of C W , in addition to the high complexity of the Feynman rules induced by the Weinberg operator, make determining the higherorder corrections a suitable real-world test case for the R * -method.

THE BACKGROUND FIELD METHOD
The renormalisation of Green's functions with a single insertion of the Weinberg operator O W requires in general a counterterm matrix Z ij . The matrix Z ij takes into account mixing with all the operators of equal or smaller mass dimension, which share the same quantum numbers as O W [42,43]: The operators that contribute to eq. (2) are divided into three classes [44,45]: • gauge invariant (GI) physical operators, • operators that vanish after applying the classical equations of motion (eom), • BRST-exact operators.
The eom and the BRST-exact operators are unphysical, since they have vanishing S-matrix elements. Nevertheless they have non-zero Green's functions and associated UV counterterms which mix with O W , as in eq. (2).
We focus on the calculation of the matrix Z ij in Yang-Mills theory without fermions (n f = 0), which probes all the gluonic diagrams. These diagrams are the most computationally demanding. The fermionic contributions mix with quark (chromo)-electric operators and a larger sector of four-quark operators [46] and will be the subject of a separate work. In pure Yang-Mills theory, there is one eom operator mixing with O W already at one loop: where D a µ is the covariant derivative. We use the background gauge field method [47][48][49][50] to simplify the mixing pattern of eq. (2). The main advantage of this method is that the 1PI Green's functions with external background gauge fields are GI and therefore the matrix Z ij involves only mixing with the GI operators. As a consequence, the counterterms proportional to O E or to BRST-exact operators vanish. Note that in general eom operators can still mix with O W if they are GI, but up to total derivatives the only independent operator we can write is which vanishes by the Bianchi identity.
In conclusion, the UV counterterms in eq. (2) must be proportional to the Weinberg operator itself, therefore the only relevant contribution is the diagonal one Z W W ≡ Z W . Beyond one loop this picture is more complicated because individual diagrams have subdivergences that do not come from diagrams with external background gauge fields and therefore give rise to gaugevariant contributions. Instead of renormalising these divergences by including the gauge-variant operators in the mixing matrix, we apply the R * -operation, which subtracts recursively all the UV subdivergences of the diagrams in a fully automated way. In this manner we isolate the local UV divergence of the correlators, which is again GI.
The renormalisation constant Z W can be extracted from the 1PI correlator Γ nb W of n background fields with a single insertion of O W -see figure 1 for examples of Feynman diagrams -by acting on it with the UVcounterterm operation Z: where Z b is the well-known wave function renormalisation of the background gauge field [49,50] and the Zoperation is defined to include counterterms for all the UV-subdivergences of Γ nb W . Here O nb W denotes the Feynman rule of the n background-field vertex generated by the Weinberg operator. The calculation of renormalisation constants can be simplified by nullifying the external momenta of Γ nb W after applying a Taylor expansion in the external momenta whose order equals the superficial degree of divergence of the correlator; for the case of Γ nb W this is ω(Γ nb W ) = 6 − n, e.g. ω(Γ 3b W ) = 3. After the external momenta are nullified, a convenient scale can be reintroduced into the correlator by introducing either arbitrarily chosen external momenta into each Feynman diagram or inserting a mass into a single propagator. This is known as the procedure of infrared rearrangement (IRR) [51]. Nullifying the external momenta introduces new IR-divergences, which are fully automatically subtracted by the local R * -operation.
In contrast to earlier works which made use of the R * -operation, we now act it on Feynman diagrams before contracting any of the Feynman rules. This makes the algorithm more efficient and allows us to directly compute the UV counterterm of a particular diagram in MS. As discussed in Ref. [27], this is not possible if the Feynman rules are contracted before the action of R * . The algorithm of [27] was therefore only capable of extracting the pole terms of self-energy diagrams from simpler self-energy diagrams. The more general method used in this work allows us to reduce the calculation of UV-counterterms of correlators with arbitrary numbers of external legs to the calculation of massless self-energy diagrams of one loop less. Therefore, this method has quite some advantages compared to the former, but its implementation introduces new complications. For instance, the Taylor expansion must be applied before the contraction of the Feynman rules, which leads to a new set of 'differentiated' Feynman rules (essentially new vertices and propagators). A comprehensive overview of the method will be given in an upcoming publication [52]. Employing this new formalism we write where the operationR * (see, e.g. [27] for a precise definition) acting on a Feynman diagram Γ subtracts from it counterterms for all UV subdivergences and all IR divergences. The operation K extracts the singleand multi-pole contributions in the dimensional regulator ǫ = (4 − D)/2 and T (ω) p1,...,pn denotes the Taylor expansion operator for the order ω-term in the expansion around the (external) momenta p 1 , . . . , p n .

CALCULATION AND RESULTS
The actual calculation of the Feynman diagrams contributing to Γ nb W , which we first generate using QGRAF [53], is done via two independent codes. The Levi-Cevita tensor appearing in the Weinberg operator is not strictly defined in D dimensions and one must fix a scheme when encountering it within dimensional regularization. In the first code, written in Maple, we use the Larin [54] scheme for the LCT appearing in the Feynman rules. The second code is written in form [55] and applies the 't Hooft-Veltman (HV) scheme [56]. The implementation of these schemes is further discussed in supplementary material.
Since the implementation of these two schemes results in rather different algorithms, obtaining a consistent result provides a powerful check. For the reduction to Master integrals both our implementations heavily rely on the forcer program [57].
To calculate the two-and three-loop anomalous dimensions of the Weinberg operator one can extract Z W from the correlator Γ nb W for n = 3-6. Ward identities ensure that the resulting anomalous dimension is independent of n. A smaller n implies a lower number of diagrams to compute. However, this is counterbalanced by the fact that the n-background field correlator must be differentiated 6 − n times with respect to external momenta in order for IRR to be applicable. Even though the Taylor expansion proliferates terms for the n = 3 case, it is nevertheless the least computationally demanding and involves 729 diagrams. At the two-loop level we obtain the result where C A is the adjoint Casimir, C A = N c for the gauge group SU(N c ), with N c the number of colors. We have checked our two-loop results in four ways. First, as discussed above we have applied two independent codes that obtain the same result. Second, in the form code we investigated gauge invariance by performing the computation with a single power of the gauge parameter ξ 1 and verified that it cancels. Third, we have extracted Z W from the n = 4 case. The evaluation of the associated 9303 diagrams leads to the same two-loop result for Z W . Finally, the 1/ǫ 2 poles can be determined from one-loop results, and we have verified that our results match the one-loop predictions.
Using the form code, which is optimized for large expressions, we evaluated the three-loop correction to Z W . This required the computation of 26443 diagrams involving up to O(10 9 ) terms in intermediate expressions. The total computation time came to 48 hours on a 24-core machine with 2.4 GHz Intel Xeon E5-2695v2 CPUs and 150 GB of memory, whereas the two-loop computation only took 20 minutes. We obtain the result Here we encounter the Riemann zeta value, ζ 3 ≃ 1.202, and the quartic Casimir d abcd A d abcd A /(N A C A ), see Ref. [58] for more details. For SU(N c ) this becomes . We checked this result by verifying that the 1/ǫ 2 and 1/ǫ 3 poles match one-and two-loop predictions. Furthermore, the fact that Z W is proportional to the three-gluon Feynman rule of the Weinberg operator is another non-trivial crosscheck.

DISCUSSION
The UV counterterm determines the anomalous dimension of the Weinberg operator up to three loops where we included the complete dependence on n f , the number of active quark flavors, at one loop. It is interesting that ζ 3 and the quartic group invariant enter at three loops, unlike the QCD beta function where they appear only at the four-loop order. However, to the best of our knowledge, there is nothing that forbids their appearance at lower order in the Weinberg operator. We notice that the coefficient of the quartic group invariant is directly proportional, up to a factor 2C A /9, to the same quartic group invariant appearing in the four-loop QCD beta function. While this could simply be a coincidence, it would be interesting to see if similar patterns appear at higher loop orders. As discussed above, the two-and three-loop results are evaluated at n f = 0. Setting n f = 0, C A = 3 we obtain the series 8πγ W α s (µ 2 )C A = 1 + 3.15657α s − 23.72872α 2 s . (10) The next-to-next-to-leading-order (NNLO) coefficient can be decomposed as 23.72872 = 5.46537 − 29.19409 where the underlined number stems from the contribution of the quartic group invariant, which is responsible for the large negative correction. The size of the coefficients increases drastically with the loop order, and undermines the convergence of the α s expansion, unless large cancellations occur in the anomalous dimension due to the diagrams involving fermion loops that are not computed in this work.
As an example, we calculate the evolution of C W , determined by from a specific high-energy scale µ H = 1 TeV, where we assume C W (µ H ) = 1, to various low-energy scales. We use α s (M Z ) = 0.118 and M Z = 91.2 GeV [59], and apply the QCD beta function at two [60][61][62] and three loops [63,64]. When evaluating the beta function we adjust the number of fermions n f at the top, bottom and charm thresholds The Wilson coefficient at the different energy scales is given in Tab. I, where we kept the LO n f dependence. Around and above the electroweak scale, α s is sufficiently small such that NLO and NNLO corrections are suppressed. For µ ≤ 100 GeV, NNLO corrections are as large, or larger, than NLO corrections. At lower energies higher orders become relevant, and in particular for µ = 1 GeV, the scale where the Weinberg operator is often matched to hadronic quantities, the NLO and NNLO correction are −21% and +33% of the LO result, respectively. The total result, however, is not far from the LO result due to cancellations between NLO and NNLO corrections.  The main phenomenological impact of a nonzero Weinberg operator is its contribution to the neutron EDM d n . The QCD matrix element connecting d n to C W is difficult to calculate, but future lattice-QCD calculations might be up to the task [65][66][67]. Two techniques have been used to estimate the matrix element. A QCD sumrule estimate [33] gives |d n | ≃ e 20 MeV C W (µ L ) with an undetermined overall sign. Another technique that is often applied is Naive Dimensional Analysis (NDA) as introduced in [68]. NDA predicts the relation [32,69] where Λ χ ≃ 1.2 GeV denotes the chiral-symmetrybreaking scale and µ match a matching scale at hadronic energies. The matching should occur where tree-and loop-level corrections are of the same order, and the scale where α s (µ match ) = 2π/3 was suggested as appropriate [32]. Using the LO or NLO anomalous dimension to evolve the Weinberg operator to the matching scale leads to such that the two-loop evolution suppresses the nEDM contribution from the Weinberg operator by a small amount. NNLO corrections completely change this picture because of the large three-loop anomalous dimension and we obtain We must conclude that higher-loop corrections destabilize the NDA predictions and that the matrix element suffers from a larger uncertainty than typically considered. Lattice-QCD calculations are direly needed.

CONCLUSION
In this letter we have reported a new method to calculate higher-order QCD anomalous dimensions of SMEFT operators in a highly automated manner. To develop and test the method we focused on one particular operator, the CP-violating gluonic Weinberg operator, whose anomalous dimension is hard to calculate, even at one loop [32,[38][39][40]. Due to its CP-violating nature this operator does not mix into lower-dimensional operators in the n f = 0 limit. By also applying the background field method, we avoid complications associated to operator mixing. We extracted the two-loop anomalous dimension by calculating 729 diagrams contributing to the three background-gluon vertex. We verified our result by extracting the same anomalous dimension of the 9303 diagrams contributing to the four background-gluon vertex, as predicted by gauge invariance. Finally, due to the automated nature of the framework we were able to immediately calculate the three-loop anomalous dimension by evaluating 26443 diagrams.
While we found a sizeable positive two-loop correction to the anomalous dimension, the three-loop correction was negative instead and sufficiently large to threaten the perturbative convergence, as the two-and three-loop evolution factors almost cancel. The lack of convergence motivates a four-loop calculation. We have discussed the impact of the anomalous dimension on EDM phenomenology, and concluded that the large anomalous dimensions destabilize often-used NDA predictions for the neutron EDM.
The R * -operation combined with the background field method provides a powerful framework for higher-order loop calculations, and our calculation can be extended into several directions without too much additional effort. In fact, our setup works up to five loops. Currently the only hurdle is computing time. Beyond five-loopsshould the need for such corrections ever arise -there exists neither a complete basis of master integrals nor a suitable reduction onto such a basis. However, the R *operation itself is valid to all loop orders.
So far, we have only performed higher-order loop calculations in pure Yang-Mills theory, i.e., without including virtual quark loops. It is crucial to develop the framework to include virtual quarks in order to get the full anomalous dimensions, and to be able to renormalise SMEFT operators containing quark fields. There are no inherent additional complications associated to the inclusion of quarks, but a consistent treatment of γ 5 in the R *method needs to be developed. This is left to future work. Here we have taken a first big step by presenting methods for the consistent use of LCTs beyond one loop in both the L and HV schemes. Once developed, we can not only include the n f -corrections to Eq. (9), but we can also calculate the (so far unknown) higher-loop mixing of the Weinberg operator into the quark electric and chromo-electric dipole moments, and CP-odd fourquark operators. It would also be interesting to compute four-loop corrections to the anomalous dimension, to see whether the trend of large coefficients continues and whether other surprising relations, as the one found for the quartic group invariants, to the QCD beta function appear.
While this work focused on the Weinberg operator, this is by no means an inherent limitation of the method. We envision calculations of higher-order anomalous dimensions of a much larger class of SMEFT operators. Isolating the footprints of SMEFT operators left behind at the LHC, or future colliders, from SM contributions is an active field of research [70,71]. Higher-order QCD corrections to SMEFT contributions can be sizeable, as shown here, and are important to disentangle SMEFT operators [72]. Higher-loop anomalous dimensions will further reduce theoretical uncertainties and make the theoretical framework of the SMEFT more robust.