Non-perturbative renormalization scheme for the CP-odd three-gluon operator

We define a regularization-independent momentum-subtraction scheme for the $CP$-odd three-gluon operator at dimension six. This operator appears in effective field theories for heavy physics beyond the Standard Model, describing the indirect effect of new sources of $CP$-violation at low energies. In a hadronic context, it induces permanent electric dipole moments. The hadronic matrix elements of the three-gluon operator are non-perturbative objects that should ideally be evaluated with lattice QCD. We define a non-perturbative renormalization scheme that can be implemented on the lattice and we compute the scheme transformation to $\overline{\text{MS}}$ at one loop. Our calculation can be used as an interface to future lattice-QCD calculations of the matrix elements of the three-gluon operator, in order to obtain theoretically robust constraints on physics beyond the Standard Model from measurements of the neutron electric dipole moment.

While the observation of a non-zero EDM in any of these experiments will be a clear indication of new physics, connecting a nuclear or atomic EDM with the fundamental, high-energy mechanism of CP violation requires gaining control over hadronic and nuclear uncertainties. At the quark level, flavor-diagonal CP violation can be model-independently described by extending the SM Lagrangian with gauge-invariant higher-dimensional operators to an effective field theory (SMEFT) [28,29], which parametrizes the indirect effects of physics at scales Λ v, where v = 246 GeV is the Higgs vacuum expectation value. Heavy SM degrees of freedom can then be integrated out by matching the SMEFT Lagrangian onto an SU (3) c ×U (1) eminvariant low-energy effective field theory (LEFT) [30,31]. The complete one-loop matching was carried out recently in [32]. At the hadronic scale, the LEFT Lagrangian includes the dimension-four QCDθ term, the dimension-five electric and chromo-electric dipole moments (CEDMs) of the u, d, and s quarks (which arise from dimension-six SMEFT operators at the electroweak scale), and, at dimension six, the CP -odd three-gluon operator and several four-quark operators. The quark-level operators induce CP -violating hadronic interactions, such as the EDMs of the neutron and proton and CP -violating pion-nucleon and nucleonnucleon couplings, which then feed into the calculations of the EDM and Schiff moments of light and heavy nuclei. As QCD is nonperturbative at the hadronic scale, the hadronic matrix elements required to match the quark-level and hadronic/nuclear EFTs need to be evaluated via nonperturbative methods. In particular, lattice QCD (LQCD) has emerged as a powerful tool to compute hadronic matrix elements, in which all sources of systematic uncertainty can be quantified, controlled, and improved. This has led to the first LQCD calculations of the nucleon EDM from the u-and d-quark EDMs [33], with few-percent uncertainties, and to the first estimates of the nucleon EDM induced by the QCDθ term [34][35][36][37][38][39][40][41] and by the quark CEDM [42,43]. These estimates, though preliminary and still affected by large uncertainties, promise to deliver controlled EDM calculations for the next generation of experiments.
An important issue to be addressed in the interpretation of LQCD results is the mixing between different CP -violating operators. The lattice spacing in LQCD effectively works as a gauge-invariant cut-off, causing mixing between operators of different dimension, in addition to the familiar logarithmic mixing of dimensional regularization. To unambigously identify the effects of various LEFT operators, it is then necessary to define a renormalization scheme. This scheme needs to be interfaced with the MS scheme, which is employed for the calculation of the operator mixing and running between the weak and the hadronic scale in the LEFT [31] and above the weak scale in the SMEFT [44][45][46]. Such a matching calculation has been performed in [47] for the dimension-five quark CEDM operator. In this paper we focus on the CP -odd three-gluon operator [28,48] where G µν is the gluon field strength, the trace is taken in color space, and G µν = µναβ G αβ /2 is the dual field-strength tensor. The three-gluon operator induces a non-zero gluon chromo-electric dipole moment (gCEDM), and we will thus also denote O (6) 1 by gCEDM. The perturbative renormalization of the gCEDM has been calculated at one loop in [49][50][51]. The anomalous dimension has recently been calculated to two and three loops in [52]. A first exploratory study of the renormalization of O (6) 1 using the gradient flow has been presented in [53]. In this work we define the three-gluon operator in a regularization-independent (RI) momentum-subtraction (MOM) scheme, construct the complete basis of gauge-invariant and nuisance operators needed to carry out the nonperturbative renormalization, and calculate the matching between the MS scheme and the momentum-subtraction scheme at one loop.
The paper is organized as follows. In Sect. 2, we discuss mixing and different types of operators that need to be taken into account in our scheme. In Sect. 3, we discuss the construction of the operators and present the resulting basis of operators that mix with the gCEDM. In Sect. 4, we define a regularization-independent renormalization scheme and in Sect. 5 we present the results for the matching between the RI and the MS scheme at one loop, before we conclude in Sect. 6. More details on the operator basis construction are provided in the appendices.

Operator mixing
In this paper we define the CP -odd three-gluon operator (1.1) in a regularization-independent momentum-subtraction scheme [54]. In this scheme, the renormalization conditions are imposed on quark, gluon, and photon Green's functions, computed in a fixed gauge, with off-shell external states of large space-like virtualities. The renormalized operators O RI j thus defined are independent of the ultraviolet regulator, and, since the renormalization conditions can be implemented both on the lattice and in perturbation theory, one can convert them into the MS scheme with matching coefficients C ij computed in perturbation theory. The implementation of the RI-MOM schemes requires working off-shell in a fixed gauge. In this case a given gaugeinvariant operator mixes with two classes of operators of the same or lower dimension [55][56][57][58][59]: I. gauge-invariant and ghost-free operators that do not vanish by equations of motion (EOM) and have the same properties as O 1 under Lorentz, chiral, and discrete symmetries (C, P , and CP ), II. "nuisance" operators, which we denote by N . These operators are allowed by the solution of the Ward identities associated with BRST invariance. They do not need to be gauge invariant, they vanish by the EOM and can be constructed as off-shell BRST variations of operators with ghost number −1, with otherwise the same properties as O 1 .
We will discuss the construction of operators in class I and II in Sect. 3, where we will further divide the operators in class II into IIa. gauge-invariant operators that vanish by EOM, IIb. gauge-variant operators.
The mixing with gauge-variant operators (class IIb) can be avoided by working in backgroundfield gauge [60].
In lattice calculations, the traditional RI-MOM scheme [54] suffers from unwanted infrared effects, which can be suppressed by choosing subtraction points with non-exceptional kinematics as in the RI-SMOM prescription [61]. In our scheme, we will impose the renormalization conditions at non-exceptional but asymmetric kinematic points (dubbed RI-SMOM scheme in [47]). As the scheme involves momentum insertion into the operator, we also need to take into account mixing with operators that are total derivatives.
We now establish the conventions used throughout the paper. If we consider only singleoperator insertions, the relation between bare and renormalized operators in any scheme is linear: where the superscript (0) denotes bare operators. By general consideration, it can be proved that the renormalization matrix has triangular structure [55][56][57][58][59] 3) The matching coefficients for the translation between the MS and RI-MOM schemes are therefore given by In this paper we consider the matching at one-loop, where (2.4) simplifies to Since the matrix elements of nuisance operators vanish between physical states [58,59], when computing hadronic matrix elements we can neglect nuisance operators. In particular, the contribution to the neutron EDM will be extracted from where the operator O MS i arises from the insertion of the effective Lagrangian, which carries no external momentum. Hence, in (2.6) the summed index j only runs over gauge-invariant, physical operators that are not total derivatives. Note, however, that in order to determine the factors ∆ ij in (2.6), either perturbatively or nonperturbatively, and in particular the renormalized operators O RI j , one is forced to also determine the mixing with nuisance operators.
We need to calculate the mixing matrices in the two schemes, which can be obtained by considering the insertions of bare operators into amputated n-point Green's functions: where ψ denotes a generic field and Z ψ its field-renormalization factor.

Construction of the operator basis
The dimension-four QCD Lagrangian is given by where the quark field includes the three light quarks, q := (u, d, s) T . We define quark-mass and charge matrices as In the LEFT, the Lagrangian (3.1) is supplemented with QED as well as a tower of effective operators [30]. Here, we will be interested in the extraction of the neutron EDM from the matrix element (2.6) with the insertion of the dimension-six three-gluon operator (1.1). We consider the matrix element at O(e) and we will neglect higher-order QED corrections. This allows us to disregard the photon kinetic term and instead treat the photon field as an external source [62,63]. In this case, the covariant derivative is given by where t a = λ a /2 and λ a are the Gell-Mann matrices in color space, and v µ , a µ , l µ , and r µ are traceless, Hermitian 3 × 3 matrices in flavor space, fulfilling and taking the physical values In the following, we will construct the basis of operators that are needed to renormalize the CP -odd three-gluon operator. The symmetries of the Lagrangian strongly constrain the possible mixings. Neglecting the QCDθ term, the leading-order Lagrangian is P -and CPeven, implying that we only need to consider CP -odd operators as possible counterterms to the three-gluon operator. In addition, in the limit M → 0 and eQ → 0, the Lagrangian has an SU (3) L × SU (3) R chiral symmetry, i.e., it is invariant under the transformation where q L,R = P L,R q and the chiral projectors are P L = (1 − γ 5 )/2, P R = (1 + γ 5 )/2. While chiral symmetry is broken by quark masses and charges, one can formally recover chiral invariance by assigning spurion transformation properties to the mass matrix and external fields. Since the three-gluon operator is chirally invariant, it can mix only with operators that are chirally invariant in the spurion sense. Chiral symmetry applies to the continuum theory. If the lattice regularization breaks chiral symmetry, additional spurions are present in the effective Lagrangian, which can induce more mixings of the three-gluon operator. We will restrict our analysis to the case where chiral symmetry is preserved by the lattice regulator. In Sect. 3.1, we briefly describe the construction of the relevant set of gauge-invariant class-I operators, while in Sect. 3.2 we explain how we construct the class-II nuisance operators. More details on the construction of the operator basis are provided in App. A and B. In Sect. 3.3, we present the complete basis of operators that are needed to renormalize the gCEDM. Based on general considerations, we discuss the structure of the mixing matrix in Sect. 3.4.

Gauge-invariant operators
We construct the basis of operators up to dimension six that renormalize the CP -odd threegluon operator. The dynamical degrees of freedom that we need to consider are the gluon field and quark fields.
In order to implement the constraints of chiral symmetry, we rewrite the non-gauge part of the leading-order Lagrangian as where The mass matrix is promoted to a spurion field and the transformations formally make the leading-order Lagrangian invariant under chiral transformations. The fieldstrength tensors associated with the external fields are with the physical values Since we are interested only in effects of O(e), we will consistently restrict the basis to operators that are at most linear in the external photon field. Gauge-invariant operators will be expressed in terms of covariant derivatives, the external field-strength tensors, the gluon field-strength tensor (3.12) and the dual field strengths with 0123 = +1. The commutator of the covariant derivative is related to the field-strength tensors by (3.14) The covariant derivative in the adjoint representation is defined by i.e., the covariant derivatives of the field-strength tensors are We use the same symbol D µ for the covariant derivative in different representations. Note that the covariant derivative fulfills the Jacobi identity as well as the Leibniz rule where each D µ denotes the proper covariant derivative belonging to the representation of the object that it acts upon. The Jacobi identity and Leibniz rule imply the Bianchi identity All three brackets have to vanish separately. By contracting the Bianchi identity with the Levi-Civita tensor, one obtains These identities play an important role in identifying the minimal set of operators that mix with the three-gluon operator. The gauge-invariant operators mixing with the gCEDM are obtained by constructing an exhaustive list of operators that are Lorentz scalars, chirally invariant, P -odd, and CP -odd, using as building blocks the fields and spurions and subsequently removing all redundancies. Details on this construction are given in App. A. The counting of operators can be automatized using Hilbert series techniques [64][65][66][67][68]. We use these methods as a cross-check to count the number of operators that are invariant under the Lorentz group (which is isomorphic to From the complete list including total derivatives and EOM operators, we select the operators that are P -odd and CP -odd.

Nuisance operators
Because of gauge fixing and the peculiar nature of BRST symmetry, by which BRST variations of elementary fields are composite operators, gauge-invariant operators can mix with noninvariant operators [55][56][57][58][59]. The form of the operators is dictated by the Ward-Slavnov-Taylor identities associated with BRST symmetry, and we follow here the construction of [58]. The construction relies on the fact that the nuisance operators can be written as BRST variations of "seed operators" with ghost number −1. The seed operators need not be gauge invariant. Their building blocks consist of the dynamical fields, the ghost fields, the spurions, as well as external sources for the fields, which are set to zero after applying the BRST variation. Details on the derivation are given in App. B. The construction provides us with a list of nuisance operators of both classes IIa and IIb. The class-IIa operators (gauge-invariant operators that vanish by the EOM) are linear combinations of gauge-invariant operators constructed in Sect. 3.1. 2 They can be presented in a compact form by introducing the fields The complete list of operators is provided in Sect. 3.3.

Operator basis
The final matrix element (2.6) that is needed to extract the neutron EDM contains an external photon state. The photon is allowed to couple either to the electromagnetic current or directly to an effective operator. As we are working at leading order in the QED coupling, we disregard operators containing more than one photon field. In a cut-off scheme, the gCEDM mixes with CP -odd operators of dimension six or lower. There are no CP -odd, chirally invariant operators with dimension smaller than four. 2 These operators correspond to the EOM redundancies that are usually removed from the set of operators in the construction of EFT Lagrangians through field redefinitions.
In the following, we present the complete basis of operators that renormalize the gCEDM operator at leading order in the QED coupling. In order to make the operators manifestly Hermitian in D dimensions, we introduce the following symbols: 3 Dimension four At dimension four, we have two physical and one nuisance operator: (3.24) The operator O 1 is the QCDθ term, the operator O 2 is a total derivative and contributes due to momentum insertion. The nuisance operator belongs to class IIa.
Dimension five At dimension five, there is a single chiral invariant operator: For a diagonal mass matrix, the following relation holds [69,70]: Dimension six At dimension six, we find the following operator basis: (3.27) where N f = 3 is the number of quark flavors. The basis in (3.27) contains the gCEDM and three additional purely gluonic operators, O 4 , a mass correction to the QCDθ term, and O (6) 5,9 , which are total derivatives. O (6) 2 and O (6) 3 are the quark CEDM and EDM, respectively. Due to the different chiral properties, the gCEDM can mix into them only via insertions of M and MQ. O (6) 6,7,10 are derivatives of the axial current, with the appropriate number of mass insertions required for chiral invariance. An important result of our construction is that there are no SU (3) chirally invariant, CP -odd four-quark operators.
In addition to the physical operators, we find a set of 20 nuisance operators at dimension six. There are 10 gauge-invariant operators that vanish by EOM (class IIa): Finally, we find another 10 gauge-variant operators (class IIb): 1 is the CP -odd three-gluon operator.

Mixing structure
In Table 1, we give the structure of the mixing matrix at leading order in the QED coupling.
The structure is determined according to the following rules.
1. Dimensional argument: operators only mix into operators of the same or lower mass dimension.
2. Operators containing the mass matrix only mix into operators with at least the same power of the mass matrix.
5. At leading order in the QED coupling, photon operators only mix into photon operators.
Due to the choice of the operator basis, the second rule involves a subtlety: because of (3.22), operators proportional to the mass matrix can be obtained from linear combinations of EOM operators with class-I operators without mass matrices. E.g., the relations imply that the qCEDM operator O 1 is a total divergence, it is the divergence of a gauge-variant current. In background-field gauge, the axial current cannot mix into this gauge-variant current and the same is applies for their divergences. Since O (4) 1,2 are gauge invariant, the same conclusion holds in any gauge. The argument of course applies as well to the dimension-six operators involving G G.

Renormalization scheme
In order to calculate the matrix element of the MS three-gluon operator in (2.6), we need the first row of the conversion matrix between MS and RI-SMOM schemes, C 1j , as well as a definition of the renormalized physical RI-SMOM operators-nuisance operators as well as total-derivative operators (apart from the topologicalθ term) do not contribute to the physical matrix element. In this section, we formulate renormalization conditions in an RI-SMOM scheme, which can be implemented in lattice QCD. They define the renormalized physical operators in terms of bare operators, hence the renormalization conditions need to determine the entries of the mixing matrix, 4 (4.1) Table 1 shows the structure of the full renormalization matrix of all the operators that potentially mix with the three-gluon operator at leading order in the QED coupling.
Only those entries of the inverse mixing matrix (Z RI ) −1 ij are needed, where both i and j run over physical operators. However, in order to determine these entries, one still has to impose n conditions if the operator O i mixes with n operators O j , j = 1, . . . , n. Let us denote the renormalization condition for insertions of the operator O i into m-point Green's functions by where | S k denotes the evaluation at a certain kinematic point and appropriate contractions in Lorentz and Dirac space defined by condition k, with k = 1, . . . , n. The desired renormalization factors are then obtained by inversion of an n × n matrix A In the MS scheme, a relation similar to (4.2) holds, with the renormalization matrix (Z MS ) −1 ij chosen to cancel only the dimensionally regulated poles. In Sect. 4.1, we compute the counterterm vertex rules for the insertions of all operators of the basis, which allows us to determine the number of independent renormalization constraints that can be obtained from a particular Green's function. The explicit renormalization conditions R k O (6) 1 on the Green's functions with insertions of the three-gluon operator are formulated in Sect. 4.3. To carry out the full renormalization program, we also need to impose renormalization conditions R k [O i ] on the Green's functions with insertions of the other physical operators O i that mix with the three-gluon operator. They can be chosen as a subset of the conditions used for the gCEDM. Alternative conditions for these operators could be obtained in a straightforward way from [47].

Counterterm vertex rules
The renormalization conditions for the RI-SMOM scheme need to render all renormalized operators finite and determine the finite contributions to the mixing matrix. We will formulate the conditions as the requirement that at certain kinematic points the renormalized amputated Green's functions agree with their tree-level expressions. In order to determine the number c, γ Figure 1. Gluon two-point, gluon three-point, and ghost-gluon three-point functions with momentum insertion into the operator. Figure 2. Quark two-point and quark-gluon three-, four-, and five-point functions with momentum insertion into the operator. of independent conditions that can be obtained from each Green's function, in the following we calculate the n-point vertex rules for all the operators of the basis (3.24), (3.25), (3.27), (3.28), and (3.29). We insert momentum q into the operator. The convention for signs and factors of i is given by We only list the contact terms. We define kinematics and indices for all the necessary Green's functions as in Figs. 1, 2, and 3. Lorentz indices are denoted by Greek letters α, β, . . . , color indices by a, b, . . . , and quark-flavor indices by i, j. All the gluon and photon momenta are incoming, while for the quark and ghost lines the momentum flow is in the direction of the fermion-and ghost-number flow.
Gluon two-point function The gluon two-point function is given by 5 + 2n (6) 8 + 2n It provides a constraint on the dimension-four operator coefficient c 1 and three independent constraints on dimension-six operator coefficients.
Quark two-point function The quark two-point function is given by It provides the two missing constraints on dimension-four coefficients, a constraint on the dimension-five coefficient, as well as 10 constraints on dimension-six coefficients.
Quark-gluon three-point function The quark-gluon three-point function is given by 9 + 2n (6) 16 + 2n 17 + 2n 6 + 3n . (4.9) It provides 10 linearly independent constraints on dimension-six coefficients. Two of them are linearly dependent with the constraints from previously listed n-point functions.
Quark-gluon five-point function The quark-gluon five-point function is given by The condition that can be obtained from it is linearly dependent with the previously listed ones.
Quark-photon three-point function The quark-photon three-point function is given by 3 + 4n 6 + 3n 2 − 2n (4.12) Out of 7 conditions on the dimension-six coefficients, 3 are linearly independent of the previously listed ones.
Quark-gluon-photon five-point function The quark-gluon-photon five-point function is given by . (4.14) It only gives a condition that is linearly dependent with previously listed ones.

Projection of scalar structures
In order to renormalize the CP -odd three-gluon operator, we need to impose 3+1+30 linearly independent renormalization conditions on Green's functions, corresponding to the counterterms from dimension-four, -five, and -six operators. The number of available structures in the Green's functions is larger-we choose to use the lower n-point functions as far as possible, which leads to the set of structures listed in Table 2. We use structures from two-, three-, and four-point functions with additional momentum insertion into the operator. The fivepoint functions are not needed. The photonic Green's functionqqA is required to provide 3 conditions that fix the photonic counterterms O 3 , O 8 , and N 2 . In the following, we define 3 + 1 + 30 projections in Lorentz, Dirac, color, and flavor space out of the Green's functions. The explicit renormalization conditions will be formulated in Sect. 4.3 by requiring these projections to agree with their tree-level expressions. We remark that all the Lorentz contractions in the projections are performed in D = 4 dimensions, see Sect. 5.2.
The quark-mass and charge matrices take the values However, after taking possible derivatives with respect to the quark masses, all conditions will be understood in the chiral limit, M → 0. This leads to a mass-independent renormalization scheme.
Gluon two-point function g 2 We evaluate the gluon two-point function with momentum insertion into the operator O i as a function of the three Lorentz invariants p 2 a , p 2 b , and q 2 . The conditions are imposed on one Lorentz contraction of the amputated two-point function and on its partial derivatives with respect to q 2 , p 2 a , and the s-quark mass at the symmetric point in the chiral limit defined by i.e., all invariants take large space-like values. Denoting by λ(a, b, c) = a 2 + b 2 + c 2 − 2(ab + bc + ca) the Källén triangle function, the projections are In R 4 , the derivative is taken with respect to the renormalized MS s-quark mass m MS s (µ = Λ). On the lattice, this can be implemented as a derivative with respect to the bare mass times the appropriate renormalization factor connecting the bare lattice mass to the MS mass.
Gluon three-point function g 3 The three-point function with momentum insertion effectively has four-point kinematics and depends on six Lorentz invariants, e.g., p 2 a , p 2 b , p 2 c , as well as the Mandelstam variables The conditions are imposed on three different Lorentz contractions of the amputated threepoint function at the non-symmetric point defined bỹ The conditions will be imposed on the following contractions: whereḡ µν denotes the metric tensor in D = 4 spacetime dimensions.
Quark two-point functionqq The quark two-point function is evaluated at the nonsymmetric point defined bỹ The renormalization conditions will be imposed on suitable contractions of the two-point function and derivatives with respect to q 2 , p 2 a , and the MS masses (again evaluated at the scale µ = Λ): Tr stands for the trace in Dirac space. The traces in flavor space are written explicitly with summed indices.
Quark-gluon three-point functionqqg The quark-gluon three-point function with momentum insertion again depends on six Lorentz invariants, e.g., p 2 a , p 2 b , p 2 c , as well as the Mandelstam variables (4.23) At the second non-symmetric point defined bŷ the renormalization conditions will be imposed on eight different contraction of the amputated three-point function: where Tr stands for the trace both in Dirac space and in color space.
Quark-gluon four-point functionqqg 2 Due to momentum insertion, the quark-gluon four-point function has five-point kinematics, i.e., there are 10 independent Lorentz invariants that can be chosen as p 2 a , p 2 b , p 2 c , p 2 d , as well as the Mandelstam variables At the non-symmetric kinematical point we use the following projections for the renormalization: (4.28) Tr again stands for the trace both in Dirac and color space.
Quark-photon three-point functionqqA For the quark-photon three-point function, we choose the same kinematical configuration as for the quark-gluon three-point function. The projections are: where Tr denotes the trace in Dirac space.

Renormalization conditions in the RI-SMOM scheme
The physical matrix element of the CP -odd three-gluon operator in the MS scheme, defined in (2.6), depends on the matching coefficients C ij and on the matrix elements of the gaugeinvariant RI-SMOM operators. One needs to provide renormalization conditions to define RI-SMOM operators on the lattice, but, as can be seen from (2.6), at O(α s ) it is not necessary to give the entries of the renormalization matrices Z MS ij and Z RI ij with i ≥ 2. In the following, we define the RI-SMOM scheme by providing explicit renormalization conditions for all physical operators. We impose 3 + 1 + 30 conditions on Green's functions with insertions of the gCEDM by requiring that the projections of two-, three-, and four-point functions defined in Sect. 4.2 agree with their tree-level values. For the renormalization of the additional physical operators, only a subset of these conditions is needed to fix all possible mixings.
The gauge-invariant operators at dimension six (3.27) include the qCEDM, O 2 , and the qEDM, O 3 . The five operators O (6) 4,6,7,9,10 are related to the divergence of the axial current and the QCDθ term, with additional powers of the external momentum or of the quark masses, which have little influence on the renormalization. As an alternative to the conditions provided below, these additional operators could be renormalized by using the conditions given in [47], where the case of generic flavor structure was discussed. With minor modifications the renormalization conditions of [47] could be adjusted to the case considered in this paper, where the flavor structure is determined by the mass M and charge Q matrices.

Conditions for the gCEDM
In order to renormalize the gCEDM, we need to impose 3 + 1 + 30 renormalization conditions. They are given as 2Nc . The coupling g on the RHS of these equations could be chosen as the renormalized coupling in any scheme. In order to simplify the matching between the MS and RI-SMOM schemes, we choose g MS (µ = Λ, ε = 0). Due to the contractions chosen in (4.28), no gluon-exchange diagrams survive in the projections of the quark-gluon four-point function and only the contact terms contribute at tree level in (4.30).

Condition for the qEDM
The qEDM operator renormalizes diagonally. Hence, it suffices to impose a single renormalization condition on the quark-photon three-point function: The coupling e to the external electromagnetic field is not renormalized in QCD.

Conditions for the remaining operators
The remaining physical operators that mix with the gCEDM are the dimension-four QCD θ term, O 1 , the dimension-five operator O 1 , and the dimension-six operator O 4 , a mass correction to the θ term. The remaining dimension-six operators are total derivatives and do not contribute for vanishing momentum insertion into the physical matrix element.
The operator O 5 . We need three renormalization conditions, one condition on the gluon two-point function and two conditions on the quark two-point function: 1 renormalizes diagonally. We impose the single condition The operator O As discussed in [47], these conditions define a renormalized GG operator that does not satisfy the singlet Ward identity. The Ward identity can be restored by a finite renormalization, as done in [47].

Matching at one loop
In this section we calculate the matching coefficients C 1j , defined in (2.4), at one loop in QCD. Since the RI-SMOM operators are independent of the chosen regulator, we can obtain the matching coefficients C 1j by calculating the n-point functions in dimensional regularization, and then imposing the MS and RI-SMOM renormalization conditions. Together with the nonperturbative definition of the RI-SMOM operators, ensured by the renormalization conditions discussed in Sect. 4, this will allow to convert lattice-QCD calculations of the nucleon EDM induced by the gCEDM to the MS scheme, up to O(α 2 s ) corrections. In Sect. 5.1, we discuss two different gauge fixing procedures: conventional covariant gauge and background-field gauge. In Sect. 5.2, we define our dimensional MS scheme. We present the results for the matching coefficients at one loop in Sect. 5.3.

Gauge fixing
We provide results with two gauge-fixing choices. First, we work in a generic covariant gauge, where the QCD Lagrangian in (3.1) is complemented by the gauge-fixing term and by the ghost Lagrangian given in (B.1). This family includes the Landau gauge ξ = 0 that can be easily implemented on the lattice. Second, we will employ the background-field method [60,75], which greatly simplifies the mixing structure. In the background-field method, all fields are split into a classical background fieldF and a quantum field F , The quantum fields are the integration variables in the functional integral. External fields and tree-level propagators are background fields, while internal loop propagators are quantum fields. For fermion fields, quantum and background fields need not be distinguished. The gauge of the background and quantum fields can be fixed independently. The backgroundfield method manifestly preserves gauge invariance with respect to the background fields, hence one only has to consider mixing with gauge-invariant operators in the classes I and IIa defined in Sect. 2, whereas no counterterms of class IIb are required. The gauge-fixing term for the quantum fields is given by whereD µ denotes the covariant derivative with respect to the background fieldĜ a µ , The background-field gauge-fixing term is simply given by where the gauge-fixing parameterξ is independent of ξ. As the background fields only appear at tree level, ghost terms can be ignored.

Dimensional regularization and renormalization
In dimensional regularization, we employ the 't Hooft-Veltman (HV) scheme [76,77]. The definition of γ 5 is where the Levi-Civita symbol µνλσ with 0123 = +1 strictly remains in four space-time dimensions. The commutation relations read In general, this scheme leads to spurious anomalies that break chiral invariance and require the introduction of symmetry-restoring counterterms [77,78]. The spurious anomalies can be traced back to higher powers of the anticommutator {γ µ , γ 5 }, which are matrices of rank D − 4 [79]. In the present case, we do not encounter these problems because QCD is a vector theory and we only consider single-operator insertions. For this reason, we do not work with chiral fields and use the D-dimensional Dirac matrix γ µ both in the QCD quark-gluon vertex and the quark propagator. External momenta and polarization vectors in S-matrix elements are treated in the HV scheme as having components only in D = 4. The same applies to the projectors that we introduced in Sect. 4.2 to define the renormalization conditions. We define the renormalization constants for the fields, coupling, and the quark masses by where γ E is the Euler-Mascheroni constant and µ denotes an arbitrary parameter with dimensions of mass, introduced to keep the renormalized coupling g dimensionless ( Note that g and α s := g 2 /(4π) depend on both µ and ε, so that dα s /d(log µ) = −2εα s + O(α 2 s ). In dimensional regularization with the MS scheme, the renormalization prescription is to subtract poles proportional to

Results
At one loop, the constants ∆ ij introduced in (2.5) can be defined as which implies that the matching coefficients are given by As can be seen from (2.7), the determination of Z ij requires the knowledge of the field renormalization. Furthermore, the renormalization conditions

Covariant gauge
We start by calculating the renormalization of the gCEDM operator in a generic covariant gauge. At one loop, it is sufficient to know the renormalization of the gluon field and of the strong coupling, where β 0 is the lowest order β function In the RI-SMOM scheme, we define the residue of the gluon propagator at p 2 = −Λ 2 to be equal to one, which determines the finite part of Z G in the chiral limit [80]: 1 . At the kinematic pointS 3 , one one-particle-reducible (1PR) diagram with the same topology as the last diagram in Fig. 5 contributes to the three-point function. The 1PR contributions to the other two legs vanish due to the renormalization condition imposed on the two-point function at the kinematic point S 2 .
The quark two-point function is zero at one loop. The quark-gluon three-and four-point functions are shown in Fig. 6. The 1PR contributions to the gluon leg of the three-point function vanish at one loop due to the renormalization conditions R 1 -R 4 . The 1PR contributions to the incoming quark leg vanish due to the renormalization conditions imposed on the quark two-point function at the non-symmetric kinematic pointS 2 . The 1PR contributions to the other quark leg contain no loop corrections, which start at two-loop order, but a coun-terterm contribution has to be taken into account since the kinematic configuration does not correspond toS 2 .
Imposing the conditions R 29 -R 31 greatly simplifies the calculation of the four-point function at one loop. Since the three-gluon vertex from the operator O (6) 1 vanishes when contracted with the gluon momentum, it is easy to see that the box diagrams in the second line of Fig. 6 do not contribute to R 29 -R 31 , leaving only the simpler triangle diagram. The same argument applies to the 1PR diagrams in the third line of Fig. 6. Of the other 1PR loop contributions to the four-point function, the gluon two-and three-point function with insertion of the gCEDM do not contribute to the projections R 29 -R 31 . The last topology shown in the second line of Fig. 6 does not receive a contribution from the loop, leaving only the contribution of the quark-gluon three-point function, contracted with the QCD three-gluon vertex. Several 1PR counterterm contributions need to be taken into account, as their kinematic configuration does not correspond to the renormalization point of the sub-amplitude.
We find the coefficients of the poles in ε to be z 1,n 8 = 1,n 11 = while all other coefficients z 1j vanish. Here, for clarity we introduced the notation z (1 − ξ) (21 + 10K + 7ψ + log 2) 1,n 16 = 1,n 20 = 59 608 with analogous notation as for the z 1j . All other c 1j vanish at one loop. The additional logarithmic term in c 1,1 is an artifact of using g MS in the renormalization conditions and disappears if the matching is performed at µ = Λ. Note also that a finite renormalization with the dimension-four operator O 1 is present. The triangle integrals in the three-point function depend on the two constants ψ and K, which are defined as where ψ (1) is the first derivative of the Digamma function. To assess the numerical impact of the conversion between the RI-SMOM and MS scheme, we can evaluate the coefficient C 11 . At the scale µ = Λ = 3 GeV, C 11 = 0.87 in Landau gauge and C 11 = 0.76 in Feynman gauge, indicating a 10% -25% correction, as to be expected at one loop.

Background-field method
We can avoid mixing with gauge-variant operators by working with the background-field method [60], with a gauge-fixing Lagrangian as specified in (5.3). In this case, the class-IIb nuisance operators in (3.29) can be disregarded and at dimension six, the operator basis reduces to the 10 operators in (3.27) and the 10 gauge-invariant EOM operators in (3.28). We can define the RI-SMOM scheme by selecting a subset of 24 conditions R k [O 1 ], with the understanding that these conditions are imposed on Green's functions of the background fieldĜ a , not of the quantum field G a . In the background-field method, we therefore replace the set of conditions (4.30) by 4, 8, . . . , 21, 24, 28, 32, . . . , 34} , In particular, the quark-gluon four-point function is no longer required. The field renormalization of the background field, which we still denote by Z G , is given by Because explicit gauge invariance is preserved, the divergent part of Z G and Z g satisfy z G = −2z g . The same relation should be retained for the finite pieces, r G = −2r g in order to preserve the Ward identity. However, we again remark that r g does not enter our matching relations, as we use the MS coupling g MS (µ = Λ, = 0) on the RHS of the conditions (5.19).
The result for the divergent pieces of the matching relations in the background-field gauge are given by 1,n 8 = It can be checked that the background-gluon two-and three-point functions respect the Ward identities implied by gauge invariance at one loop for any value of the quantum gauge parameter ξ. However, even in the background-field method, off-shell Green's functions are unphysical quantities and they are gauge-parameter dependent. Therefore, our matching relations depend on the quantum gauge parameter ξ (but of course not on the background gauge parameterξ), since the RI-SMOM conditions themselves are gauge dependent.

Conclusions
The CP -odd three-gluon operator gives the main contribution to the nucleon EDM in several beyond the Standard Model scenarios, especially when CP is violated in the interactions of heavy particles, such as the Higgs [48] or the Higgs and the top quark [81]. First-principle calculations, with controlled theoretical uncertainties, of the matrix elements of the gCEDM on the nucleon are necessary to derive the constraints of EDM experiments on this operator, and the implications for BSM physics. At the moment, the best estimates of the nucleon EDM from the gCEDM have been obtained with QCD sum rule calculations [82,83], which are however affected by large theoretical uncertainties, at the level of 50%-100%. While for this operator lattice QCD calculations are still in their infancy [53,84], this method can in principle provide fully nonperturbative results, in which all sources of systematic uncertainty can be quantified, controlled, and improved. LQCD and continuum calculations are interfaced via the definition of a renormalization scheme. In this paper, we have defined a RI-SMOM scheme for the renormalization of the gCEDM, and we have provided the conversion matrix to the MS scheme at O(α s ).
As a dimension-six, flavor-singlet operator, the gCEDM has a complicated mixing pattern in an off-shell scheme. On the lattice, insertions of O 1 induce power divergences, which under the assumption of good chiral symmetry can be absorbed by three dimension-four and one dimension-five operator, defined in (3.24) and (3.25). Both on the lattice and in the continuum, the gCEDM mixes into 10 dimension-six gauge-invariant operators that do not vanish by EOM, given in (3.27), 10 gauge-invariant nuisance operators (3.28), and 10 gauge-variant nuisance operators (3.29). In this work, we have provided 34 renormalization conditions that define our RI-SMOM scheme. In order to obtain enough independent conditions, it is necessary to compute the gluon two-and three-point functions, the quark two-point function, the quarkgluon and quark-photon three-point functions, and some projections of the quark-gluon fourpoint function. We have imposed the renormalization conditions at one loop and computed the conversion matrix between the RI-SMOM and MS schemes, both in a conventional covariant gauge and in background-field gauge.
The number of operators and renormalization conditions make the lattice implementation of the RI-SMOM renormalization scheme challenging, even though calculations of comparable complexity have been carried out for ∆S = 1 operators that contribute to K → ππ decays [85,86]. For this reason, we explored the definition of the RI-SMOM scheme in the backgroundfield gauge [60], which allows to discard gauge-variant operators. Using the background-field method, the definition of the RI-SMOM scheme involves only two-and three-point functions, a very noticeable simplification. While background-field methods have not extensively been used to study higher-dimensional operators on the lattice, there are no particular technical problems for the implementation of the background-field condition [87], and thus of the renormalization conditions enumerated in Sect. 5.3.2. It will be interesting to further explore the use of background-field methods in actual numerical simulations.

A Construction of gauge-invariant operator basis
In this appendix, we provide details on the construction of the basis of gauge-invariant operators. In App. A.1, we describe the symmetries of the building blocks. In App. A.2, A.3, and A.4, we construct a complete list of pure gauge operators, two-quark, and four-quark operators, respectively, which we summarize in App. A.5. Here, we disregard evanescent operators away from D = 4 dimensions, which will be discussed in App. C.

A.1 Symmetries and building blocks
The gauge-invariant class-I operators that are needed to renormalize the CP -odd three-gluon operator are constructed from the building blocks (3.21). The mass matrix has been promoted to a spurion field. The chiral transformations (3.9) assigned to spurion and external fields allow us to take into account explicit chiral-symmetry breaking.
In order to renormalize the three-gluon operator, the operators have to be chirally invariant in the spurion sense, Lorentz scalars, CP -odd, and P -odd. Since we are working at leading order in the electromagnetic coupling and the external photon field always comes together with a charge matrix Q and a coupling e, we only consider operators with at most one QED field-strength tensor.
The symmetry properties of the building blocks are listed in Table 3. The chargeconjugation matrix fulfills and can be written in the Dirac representation as C = iγ 2 γ 0 , hence Table 3. Properties of dynamical fields, spurion and external fields, and derivative operators. For simplicity, additional arbitrary phases in P -and CP -conjugation are neglected. η(a) is defined in (A.6). The Lorentz group is locally isomorphic to SU (2) L × SU (2) R .

The Dirac field transforms under charge conjugation as
where η c is a phase factor, which we put equal to 1 in the following. The electromagnetic gauge field transforms as whereas the non-abelian gauge field transforms under charge conjugation as We classify the operators according to the field content and mass dimension (up to O(e)): • pure gauge operators: • two-quark operators: • four-quark operators: where ψ 2 denotes a quark bilinear. In this list, we have already excluded classes that obviously contain no gauge-invariant operators, e.g., F G, F G 2 classes. In the following, we construct the explicit operators by hand. We use the Hilbert series techniques [64][65][66][67][68] as a cross-check to count the number of operators in each class, including total derivatives and EOM operators.
As the known Hilbert series method does not include the discrete symmetries, even in this cross-check we select by hand the operators that are P -odd and CP -odd.

A.2 Pure gauge operators
As a first class of operators, we consider the pure gauge operators. The building blocks are the quark-mass matrix, partial and covariant derivatives, and field-strength tensors. There are no operators at dimension two or three, hence we start at dimension four. Note that we need at least two field-strength tensors in order to have a non-vanishing trace. As we are only interested in operators up to O(e), we disregard the electromagnetic field-strength tensor in this section. Although in dimensional regularization mixing is only possible within operators of the same mass dimension (the mass matrix is treated as a spurion field), this is not necessarily true for other schemes. Therefore, we also look for P -odd, CP -odd operators of dimension smaller than six. dim = 4 At dimension four, we have two operators that consist only of gauge fields: The first term is the standard CP -even kinetic term for the gauge field, the second one is the CP -odd and P -odd QCD θ-term. It belongs to our basis: Here, we use the tilde to distinguish a preliminary set of operators from the final ones after having removed redundancies. dim = 5 At dimension five, there are no chirally invariant pure gauge operators. dim = 6 We reach dimension six by adding either two mass matrices or two derivatives to a dimension-four operator (adding one mass matrix and one derivative does not give a Lorentz scalar). The only way to add two mass matrices is within a trace: where now the first trace is in flavor space, the second one in color space.
We consider the addition of two partial derivatives. The six Lorentz indices can be contracted either with g µν g λσ g αβ or with µνλσ g αβ to form a Lorentz scalar. In total, there are only four different contractions: Furthermore, the Schouten identity g αρ µνλσ + g αµ νλσρ + g αν λσρµ + g αλ σρµν + g ασ ρµνλ = 0 (A.11) implies the relation This only leaves the following P -odd and CP -odd operator: Next, we consider the case where we add one partial and one covariant derivative. The possible contractions of the Lorentz indices are: However, we only have to consider contractions with the Levi-Civita tensor: they are P -odd and CP -odd, while the other contractions are even. By applying the Bianchi identity (3.19), we remove redundancies. Furthermore, we note the Leibniz rule (3.18) for the covariant derivative in adjoint representation: We find the relations hence, there is only one additional independent operator: Finally, we can build operators with two covariant derivatives and two field-strength tensors. The requirement that the operator be P -and CP -odd allows again only the contraction with µνλσ g αβ . If the two derivatives do not act on the same field-strength tensor, we can use (A.15) and obtain a linear relation to an operator where both derivatives act on the same tensor and an operator involving a partial derivative.
The possible contractions are We take some linear combinations to replace this set by Using the Jacobi and Bianchi identities, we can eliminate three elements of the set: The commutators of covariant derivatives in the adjoint representation can be expressed in terms of the field-strength according to Therefore, the only additional operator is the CP -odd three-gluon operator itself, which is the only P -odd and CP -odd operator that can be constructed with three field-strength tensors. This completes the construction of the set of pure gauge operators.

A.3 Two-quark operators
We continue with two-quark operators, which at least have mass dimension three.

dim = 3
There is no quark bilinear that is a Lorentz scalar and chirally invariant. dim = 4 In order to reach mass dimension four, we can add either one mass matrix or one derivative to a quark bilinear. The chirally invariant operators obtained by adding a mass matrix to a quark bilinear arē There is one Hermitian linear combination that is P -and CP -odd: Next, we consider the insertion of a derivative, which has to be contracted with a Lorentzvector, hence we need a vector quark bilinear. The possible gauge-invariant and chirally invariant contractions with a partial or covariant derivative are the following: The third and fourth operators are CP -even (they are the standard kinetic terms), while the first two are CP -odd. There is only one linear combination that is also P -odd, the divergence of the axial current:Õ 2q,(4) 2 whereγ µ is defined in (3.23).

dim = 5
We reach mass dimension five by inserting two masses, one mass and one derivative, two derivatives, or a field strength tensor into a quark bilinear. There are two chirally invariant operators obtained from the insertion of two mass matrices into a quark bilinear: which are possible due to the fact that for SU (3), the following tensor decompositions hold: hence a product of three (anti-)fundamental representations contains a singlet. If the spurions M and M † are fixed to a diagonal mass matrix, the above operators are flavor conserving as well and represent a correction to the mass terms themselves. We can again form one P -odd and CP -odd combination: For a diagonal mass matrix, the following relation holds [69,70]: Consider the insertion of a single mass matrix and a derivative into a quark bilinear. In order to contract the Lorentz index of the derivative, we need a vector bilinear. With an additional mass matrix, it is impossible to construct a chirally invariant operator. Finally, we consider the case of two derivatives. We have to start either with a (pseudo-)scalar or with a tensor quark bilinear, and add two derivatives. Also here, we cannot construct a chirally invariant operator. The same is obviously true for the insertion of a field-strength tensor in a two-quark operator. dim = 6 We obtain operators of dimension six by inserting either three mass matrices, two mass matrices and one derivative, one mass matrix and two derivatives, or three derivatives into a quark bilinear. Furthermore, a field-strength tensor can take the role of two derivatives.
We start with the insertion of three mass matrices. A basis for the chirally invariant operators is given by: The SU (N f ) Fierz identity hence t A ⊗ t A operators are linearly dependent of 1 ⊗ 1 operators. We can form the following Hermitian P -odd and CP -odd linear combinations: Next, we insert two mass matrices and one derivative into a quark bilinear. We find the following chirally invariant Hermitian operators: The four operators with covariant derivatives i ← → D µ are CP -even. The following linear combinations are P -and CP -odd: The next operator class consist of insertions of one mass matrix and two derivatives in a quark bilinear. We start with the following set of chirally invariant operators: The following Hermitian linear combinations are P -and CP -odd: whereσ µν is defined in (3.23). The next class of two-quark operators contains insertions of three derivatives. The three Lorentz indices of the derivatives {·} µνλ can either be contracted with g µν γ λ or with µνλσ γ σ . For the moment, we disregard evanescent operators (see App. C) and use the four-dimensional relation γ µ γ ν γ λ = g µν γ λ + g νλ γ µ − g µλ γ ν + i µνλσ γ σ γ 5 .
(A. 39) Due to the odd number of gamma matrices, all operators will be chirally invariant, hence we work directly in the parity basis. We start with the contractions with g µν γ λ . A γ 5 matrix is required for P -odd operators. The derivatives can be either covariant derivatives or partial derivatives of a gauge singlet. Note that due to the covariant derivatives acting on the left can always be put on the left-hand side of derivatives acting on the right. Furthermore, by using the relation left-acting covariant derivatives can be traded for partial derivatives of the gauge singlet. Hence, we find the following list of nine operators with three derivatives: (A.42) By taking linear combinations, we make them manifestly Hermitian: Four operators are CP -odd:Õ 2q,(6) 9 Next, we investigate the contractions of three derivatives with µνλσ γ σ . It is only possible to have three covariant derivatives or two covariant and one partial derivative: partial derivatives are commuting, hence two or three of them vanish upon contraction with the Levi-Civita tensor. Left-acting covariant derivatives can again be traded for right-acting and partial derivatives. If we choose two covariant and one partial derivative, we can immediately insert the commutator of the covariant derivatives. In order to have a P -odd operator, no γ 5 matrix is allowed. The only Hermitian operator with two covariant derivatives is thereforẽ O 2q,(6) 13 which indeed is CP -odd. Finally, consider the insertion of three covariant derivatives: Hermitian conjugation of this operator is identical to a CP conjugation. The CP -odd Hermitian component is again identical toÕ 2q,(6) 13 , which is therefore the only CP -odd operator. Finally, we consider the operator classes with field-strength tensors. Due to (3.14), we only need to take into account the external (electromagnetic) field-strength tensor: the QCD field-strength tensor can be written as a linear combination of the commutator of covariant derivatives and the external field-strength tensors. In the class ψ 2 F M, the chirally invariant operators arē One Hermitian linear combination is both P -odd and CP -odd: The last class of two-quark operators is ψ 2 F D. Here, we find the P -odd operators Two CP -odd Hermitian linear combinations exist:

A.4 Four-quark operators
A basis for the chirally invariant four-quark operators is given by the following six operators: where t a are the SU (3) generators in color space. Note that flavor-octet operators 1c8 f , and O V,RR 8c8 f are related to the above operators through Fierz identities in Dirac space as well as the SU (N f ) and SU (N c ) Fierz relations: On the other hand, the flavor-octet LR operators O V,LR 1c8 f , O V,LR 8c8 f are not chirally invariant. For the four-quark operators, Hermitian conjugation acts in the same way as a CPtransformation. All the operators O V,LL , O V,RR , and O V,LR are Hermitian and CP -even. We conclude that there is no four-quark operator that could mix with the CP -odd three-gluon operator.
Note that if SU (2) instead of SU (3) chiral symmetry is considered, there are additional chirally invariant operators [88] due to the absence of the symmetric structure constants d ABC .

A.5 Intermediate summary
Here, we summarize the P -odd, CP -odd, Lorentz-and gauge-invariant, chirally invariant Hermitian operators up to dimension 6.
Pure gauge operators At dimension four, there is the QCD θ-term: while at dimension six, we find four operators: Four-quark operators There are no four-quark operators that can mix with the gCEDM.

B BRST invariance and nuisance operators
In this appendix, we provide details on the construction of the nuisance operators, which vanish by the EOM. We follow the method of [58]. In App. B.1, we review the EOM. In App. B.2, we discuss the Slavnov-Taylor identities. The recipe for the construction of the nuisance operators is reviewed in App. B.3. The symmetry properties of the building blocks are discussed in App. B.4. We construct the seed operators in App. B.5 and list the resulting nuisance operators in App. B.6. In App. B.7, we derive redundancies in the preliminary operator set, and we relate the operators to our final basis.

B.1 Gauge fixing and equations of motion
The QCD Lagrangian including gauge fixing and Faddeev-Popov ghosts is given by where D ac µ = ∂ µ δ ac + gf abc G b µ is the covariant derivative in the adjoint representation. The gauge-fixing term can also be written in terms of an auxiliary field G a [89]: The EOM for G a is ξG a = ∂ µ G a µ , which, inserted into L GF , leads to the original gauge-fixing term plus a total derivative.
When the quark mass matrix is promoted to a spurion field, the mass term has to be replaced by The complete list of EOM reads: For notational convenience, we define the EOM fields: The definition of the EOM quark fields is chosen in such a way that they fulfil q † E =q E γ 0 .

B.2 Slavnov-Taylor identities
We add source terms for the fields and the (composite) BRST variations: The action including a set of sources for gauge-invariant ghost-free operators O i is defined as We introduce the BRST transformation: where δλ is an anticommuting infinitesimal parameter. With this transformation, we find δ(gc a t a q) = gδc a t a q + gc a t a δq = 1 2 and similarly δ(qgc a t a ) = 0 , (B.10) as well as where in both relations we have used the Jacobi identity for the SU (3) structure constants. Therefore, we see that for all fields φ ∈ {G a µ , q,q, c a ,c a , G a } i.e., the BRST transformation is nilpotent. Note that if the auxiliary field G a is not used, nilpotency for anti-ghosts only holds on-shell [90].
For the fields G a µ , q,q, the BRST transformation corresponds to an infinitesimal gauge transformation with the parameter a = c a δλ. Therefore, the physical part of the Lagrangian is invariant under BRST transformations. For the ghost and gauge-fixing part, one finds: The source terms for the BRST variations are obviously invariant as well, hence the only variant part of the Lagrangian are the source terms for the fields: δL 0 = H µ,a δG a µ +Lδq + δqL +N a δc a + δc a N a = H µ,a D ac µ c c − igLc a t a q + igqc a t a L + (B.14) The generating functional is invariant under the variable transformation (see [89] for the invariance of the measure) which implies (J generically denoting the sources) We introduce the effective action as the Legendre transform of W , which is the generating functional for one-particle-irreducible truncated Green's functions (we do not transform the auxiliary field, which is not propagating): where the "classical fields" are expectation values (the functional derivatives are understood to act from the left): The variations of Γ with respect to classical fields and sources are given by: where N are additional nuisance operators. Working to first order in the external sources Φ, Φ, one finds that the nuisance operators satisfŷ with the operator The BRST operatorŴ is nilpotent,Ŵ 2 = 0, and it carries ghost number 1. In [57] it was shown that the most general solution for the nuisance operators is given bŷ where F is a set of anti-Hermitian "seed operators" with the same Lorentz, chiral, and global SU (3) c properties as O, the same discrete symmetries and dimension, and ghost number −1.
Hence, in the following we construct systematically the set F and derive from it the (gaugevariant) nuisance operators: we act with the operatorŴ on the set F and afterwards we set the sources to zero. Let us work out the explicit form of the operatorŴ . With sources already set to zero, we have This leads tô Note that after acting with the term on the seed operators, the EOM quark field q E needs to be anticommuted to the right-hand side, which produces an additional minus sign.
In case that the seed operator contains derivatives of the sources J µ a , . . ., we have to use partial integration, e.g.,

B.4 Symmetry properties of sources and building blocks
In Table 4, we list the transformation properties of the various fields and sources, which are the building blocks for the seed operators. In particular, the given transformation properties ensure that the leading-order Lagrangian is Hermitian, 5 P -even, and CP -even. Note that we define the complex conjugate of the product of Grassmann variables c 1,2 as (c 1 We assign zero mass dimension to the ghost field and mass dimension 2 to the anti-ghost field. In this convention, the operator d 4 xŴ (x) does not change the mass dimension of the seed operators. If a mass dimension 1 is assigned to both ghost and anti-ghost fields, the operator d 4 xŴ (x) raises the mass dimension by one unit. The assignment of the mass dimensions is purely conventional and does not affect the results.
We assign the following chiral transformations:  Due to gauge fixing, the seed operators need not be gauge invariant under SU (3) c . Therefore, the gauge field G a µ is allowed as a separate building block and not only as part of the full covariant derivative D µ . However, the gauged SU (3) L × SU (3) R chiral symmetry, which also contains U (1) em as a subgroup, remains intact. Therefore, we define covariant derivatives with respect to the external fields only, and impose on the seed operators invariance under the local chiral group.

B.5 Seed operators
Let us now systematically construct the gauge-variant seed operators using the building blocks in Table 4. Due to the ghost EOM (B.23), the anti-ghost and the source J µ a only appear as a building blockĴ µ,a := J µ,a − ∂ µca . (B.38) After applying the operatorŴ , we will set the sources M L,R ,M L,R , K a ,K a to zero, hence we only need to take into account seed operators with at most one of these sources. The sourceĴ µ,a should be set to −∂ µca .
In order to construct SU ( in the case of four indices [91,92].
K a operators At dimension 4, the only operator c a K a is P -and CP -even. At dimension 5, there are no operators. At dimension 6, the following operators have to be considered: c a ( K a ), even, even, ( c a )K a , even, even, Therefore, the only anti-Hermitian P -and CP -odd seed operator at dimension 4 is No dimension-5 operator can be constructed. At dimension six, we find the following list of chirally invariant operators. Note that we already neglect operators that vanish when the spurions and external fields are fixed to their physical values: seed operator P CP i(q L,R σ µν t a M R,L +M R,L σ µν t a q L,R )∂ µ G a ν , L ↔ R, even, This results in the following list of 19 anti-Hermitian P -and CP -odd seed operators: , even, even, , even, even, q L γ µ t a q LĴ a µ ,q R γ µ t a q RĴ a µ , L ↔ R, even.
Dropping all operators that are not both P -and CP -odd, we are left with three anti-Hermitian seed operators:

B.6 Nuisance operators
After acting with theŴ operator on the seed operators F i , we obtain the list of nuisance operators. Furthermore, we perform a basis change: we write as many nuisance operators as possible in a manifestly gauge-invariant form and as total derivatives. In order to write the operators in a more compact form, we fix the spurion and external fields to their physical value, M, M † → M, l µ , r µ → eQA µ , and write everything in the parity basis. There is one nuisance operator at dimension four: At dimension six, we find the following nuisance operators:

B.7 Redundancies
The nuisance operators are constructed as BRST variations of a complete set of linearly independent seed operators. However, it turns out that the nuisance operators themselves are redundant. This can be understood as follows: acting with the BRST operator on the seed operators replaces the BRST sources by EOM fields, i.e., this operation reduces the degrees of freedom. Therefore, linear independence of the seed operators does not imply linear independence of the resulting nuisance operators. The remaining redundancies are most easily identified by considering the vertex rules for all the operators, see Sect. 4.1, which leave two linear combinations of nuisance operators undetermined. An explicit calculation then confirms the following linear relations: 12 , 0 = gÑ 9 . In a last step, we remove redundancies from the list of operatorsÕ in App. A.5: there are linear combinations that are identical to nuisance operators. Removing these redundancies leads to a minimal set of class-I operators O, i.e., gauge-invariant operators that do not vanish by the EOM.

C.1 Generalities
As is well known, dimensional regularization leads to the appearance of evanescent operators [59,[93][94][95]. These operators are present in D dimensions, but they vanish for D = 4. If bare evanescent operators are inserted into loop diagrams, the combination of poles in with the evanescent structure can lead to finite contributions. Evanescent operators can be renormalized by finite counterterms so that the renormalized evanescent operators have vanishing matrix elements. As shown in [94,95], the renormalized evanescent operators do not mix into physical operators. However, the counterterms affect the calculation of the anomalous dimension. Since the bare evanescent operators are ambiguous, their choice affects the anomalous dimension matrix of the physical operators and their definition is part of the scheme. In dimensional regularization, the evanescent operators are present both in the MS scheme and in the MOM scheme. Let us denote the relation between bare and renormalized operators as Since in general E tree = O(ε), one expects that either the evanescent counterterms z OE need to be determined separately, or the number of conditions to be imposed must match the number of physical plus evanescent operators. However, this can be avoided if the renormalization conditions and the set of evanescent operators is chosen so that R[E (0) ] tree = 0 instead of R[E (0) ] tree = O(ε). This allows us to consider only as many conditions as physical operators are present and to solve the system for the coefficients c OO that determine the conversion matrix.
In the following, we define the set of relevant evanescent operators and show that their tree-level insertions into the renormalization conditions vanish identically. Throughout, we use the HV scheme [76,77] to deal with the Levi-Civita symbol and the Dirac matrices in D = 4 spacetime dimensions. In defining our scheme, it is useful to divide the operators in our basis (up to and including dimension six) in two categories: (i) purely bosonic operators O B , involving one gluonic dual field strength (such as gCEDM); (ii) fermionic operators O F containing a quark bilinear with a Dirac structure involving one γ 5 , and possibly gluonic structures, the external electromagnetic field, and derivatives (such as the pseudoscalar density and the qCEDM).

C.2 Definition of evanescent operators
The bosonic operators at dimension four and six can be written schematically as where O B are Lorentz tensors of rank four and six, respectively, built out of ∂ µ , G a µ , and A µ . In the HV scheme, the indices of the Levi-Civita symbol are restricted to D = 4 dimensions. In addition, external momenta and polarization vectors in S-matrix elements are considered to be objects in D = 4 dimensions. However, the restriction to D = 4 of external momenta and polarizations can be performed after performing the loop calculation. As we are considering only QCD corrections, all vertices and propagators in loops are continued to D dimensions. Therefore, any metric tensor that appears in a loop calculation (either from propagators, tensor reductions of loop integrals, or the Dirac trace of closed fermion loops) is D-dimensional. In particular, an evanescent structure where the indicesα andβ are restricted to −2ε dimensions, cannot be independently generated in QCD with single insertions of the gCEDM operator. This implies that the only evanescent bosonic operator appears at dimension six due to the Schouten identity (A.11): where the indices µ, ν are in 4, and α in D dimensions.
For the fermionic operators, we make use of the HV definition (5.7) for γ 5 and the fact that in any spacetime dimension D a string of k Dirac matrices can be decomposed as a linear combination of the fully antisymmetric products Γ (n) α 1 ...αn = 1 n! γ [α 1 ...γ αn] . (C.14) In particular, the product of n gamma matrices is expressed as a combination of Γ (n) with n = k, k − 2, ..., with the remaining Lorentz indices provided by appropriate powers of the D-dimensional metric tensor g µν [96]. Note that all structures involving Γ (n) with n > 4 are evanescent.
The fermionic operators at dimension four, five, and six can be written schematically as ..λn q g µaν b g µc λ d g νeν f g νg λ h , k = 0, 1, 2, 3 , n = 4 − k, 6 − k, . . . , 4 + k , (C. 15) where O F is built out of ∂ µ , G a µ , A µ , color structures, and the charge and mass matrices. The highly symbolic product of metric tensors is needed to ensure the final result is a Lorentz scalar (note that µ a and λ b indices cannot be contracted among themselves due to the antisymmetry of the Levi-Civita symbol and Γ (n) ). As in the case of the bosonic operators, QCD loops only generate metric tensors in D dimensions, since we only consider single insertions of the gCEDM. Therefore, the only possible evanescent operators arise either due to the evanescent structures Γ (n) with n > 4, or due to the Schouten identity. An explicit list (with generic operators O F ) is given by objects, in the same way as it happens with S-matrix elements in the HV scheme. When the evanescent operators are inserted into the truncated vertex functions at tree level, the fields are removed and derivatives turn into four-dimensional external momenta. All indices are either contracted with the four-dimensional Levi-Civita symbol or the four-dimensional indices of the projectors, hence both the Dirac structures Γ (n) and the metric tensors in the evanescent operators are projected to D = 4 dimensions. Therefore, both the four-dimensional Dirac algebra and the Schouten identity apply and the tree-level insertions of the evanescent operators vanish identically.