Short flow-time coefficients of CP-violating operators

Measurements of a permanent neutron electric dipole moment (EDM) potentially probe Beyond-the-Standard Model (BSM) sources of CP-violation. At low energy the CP-violating BSM interactions are parametrized by flavor-conserving CP-violating operators of dimension higher than four. QCD calculations of the nucleon matrix elements of these operators are required to fully reconstruct the sources and magnitudes of the different CP-violating contributions to the nucleon EDM. Herein we study the quark-chromo electric dipole moment (qCEDM) operator and the three-gluon Weinberg operator. The non-perturbative determination, using lattice QCD, of the nucleon matrix elements of these CP-violating operators is hampered by their short-distance behavior. Under renormalization these operators mix with lower dimensional operators, which induces power divergences in the lattice spacing, as the continuum limit is approached. We study the short-distance behavior of the qCEDM and the Weinberg operators using the gradient flow. We perform a short flow time expansion and determine, in perturbation theory, the expansion coefficients of the linearly-divergent terms stemming from the mixing with the pseudoscalar density and the topological charge, confirming the expectations of the operator product expansion. We introduce a new method to perform calculations at non-zero flow-time for arbitrary values of the external momenta. This method allows us to work in four dimensions for most of the calculations described in this paper, avoiding the complications associated with defining $\gamma_5$ in generic d dimensions. We show that leading contributions in the external momenta can be reproduced by defining $\gamma_5$ using the 't Hooft-Veltman-Breitenlohner-Maison scheme.

Measurements of a permanent neutron electric dipole moment (EDM) potentially probe Beyondthe-Standard Model (BSM) sources of CP-violation. At low energy the CP-violating BSM interactions are parametrized by flavor-conserving CP-violating operators of dimension higher than four. QCD calculations of the nucleon matrix elements of these operators are required to fully reconstruct the sources and magnitudes of the different CP-violating contributions to the nucleon EDM. Herein we study the quark-chromo electric dipole moment (qCEDM) operator and the three-gluon Weinberg operator. The non-perturbative determination, using lattice QCD, of the nucleon matrix elements of these CP-violating operators is hampered by their short-distance behavior. Under renormalization these operators mix with lower dimensional operators, which induces power divergences in the lattice spacing, as the continuum limit is approached. We study the short-distance behavior of the qCEDM and the Weinberg operators using the gradient flow. We perform a short flow time expansion and determine, in perturbation theory, the expansion coefficients of the linearly-divergent terms stemming from the mixing with the pseudoscalar density and the topological charge, confirming the expectations of the operator product expansion. We introduce a new method to perform calculations at non-zero flow-time for arbitrary values of the external momenta. This method allows us to work in four dimensions for most of the calculations described in this paper, avoiding the complications associated with defining γ5 in generic d dimensions. We show that leading contributions in the external momenta can be reproduced by defining γ5 using the t Hooft-Veltman-Breitenlohner-Maison scheme.

I. INTRODUCTION
The nucleon electric dipole moment (EDM) is a physical quantity that, once measured, will provide a unique opportunity to detect and investigate beyond-the-standard model (BSM) sources of charge and parity (CP) violation. In principle, there are multiple sources for a non-vanishing nucleon EDM, including the Cabibbo-Kobayashi-Maskawa (CKM) quark-mixing matrix, the quantum chromodyamics (QCD) θ term, higher-dimensional CP-violating operators, or any combination of these. The current experimental limit for the neutron EDM [1,2], |d n | ≤ 1.8 × 10 −26 e cm at 90% confidence level, leaves open the possibility of a dominant BSM source of CP-violation, which could be several orders of magnitude larger than Standard Model sources (see Ref. [3] for a recent review).
In addition to the Standard Model contributions to the nucleon EDM from the CKM matrix [4] and from the θ term [5], BSM theories that contain complex CP-violating couplings can induce a non-vanishing EDM at the one loop level. At low energies the BSM degrees of freedom are heavy enough that one can parametrize their effects through effective, higher-dimension CP-violating operators. In this paper we consider two such operators, the quarkcolor EDM (qCEDM) operator and the CP-violating three-gluon operator, i.e. the Weinberg operator. To constrain couplings in BSM theories at high energies, one needs to determine the QCD contribution to the EDM at low energy. Broadly speaking, there are three approaches to determining the relevant matrix elements: QCD sum rules [6,7]; chiral perturbation theory [8,9]; and lattice QCD.
Lattice QCD provides the most systematic method to calculate individual contributions from different CP-violating sources to the nucleon EDM in terms of the QCD fundamental degrees of freedom, quark and gluons. There is a long history of attempts to determine the nucleon EDM from lattice QCD [10][11][12][13][14][15][16][17][18], and several technical difficulties have been encountered.
The first difficulty arises from the fact that in Euclidean space the θ term renders the QCD action complex, which prevents the use of stochastic methods. The current experimental bound on the neutron EDM implies a very small value for θ ∼ 10 −10 , justifying a perturbative expansion in θ. Correlators that include an insertion of the θ term, once the topological charge has been properly renormalized, are theoretically well-defined. Despite the very poor signal-to-noise ratio it is possible to determine the nucleon EDM induced by the θ term using signal-to-noise improved ratios [5,19].
The second difficulty arises from the renormalization of the relevant composite operators. In Ref. [20] we proposed using the gradient flow [21][22][23][24] to renormalize the θ term and the BSM CP-violating operators. We are currently pursuing this program and in Refs. [5,13,19,25] we investigated and calculated the nucleon EDM from the θ term.
Renormalization schemes based on the gradient flow include non-perturbative step-scaling approaches [55,56], removing power divergences in nonlocal operators relevant to hadron structure [57,58], and defining regularizationindependent quark-bilinear currents [59,60]. Perturbative calculations of the gradient flow have been carried out to three loops for certain quantities using automated perturbation theory routines [43,61,62] and to two loops via numerical stochastic perturbation theory [63,64].
Analytic loop-order calculations with the gradient flow often introduce some difficulties related to dimensional regularization. One method to avoid these complications employs an expansion in the external momentum p to some desired order. This can induce extraneous, non-physical infrared poles at fairly low orders in the external momentum. (In the calculation of the Wilson coefficient c CP below, for example, these appear as early as O(p 2 ).) We have used a novel combinatorial scheme to track the external momentum at all orders, which maintains finiteness at positive flow time throughout all of the calculations in this paper with the exception of those related to the renormalization of the flowed fermion propagator in App. C.
In this paper we focus on the renormalization of the higher dimensional CP-violating operators using the gradient flow. First results appeared in [65][66][67] and in this paper we determine the leading contribution to the short flow-time expansion (SFTE) coefficients of the CP-violating operators defined using the gradient flow. The renormalization and mixing of these operators was first studied in the MS scheme in [68][69][70], the determined at two loops for the qCEDM in [71], and at two and three loops for the Weinberg in [72]. After describing our perturbative strategy for determining these coefficients, we focus on the leading linearly-divergent expansion coefficients and some logarithmic terms.
The paper is organized as follows. We first introduce the gradient flow and some technical details relevant for our perturbative expansion in Sec. II. We calculate the expansion coefficients of the quark-chromo operator, parameterizing the mixing with the pseudoscalar density and the topological charge density, in Sec. III, and the corresponding coefficient of the Weinberg operator, induced by the mixing with the topological charge density, in Sec. IV. We summarize our results and our conclusions in Sec. V.
In Appendix A we detail our notations and conventions including the d−dimensional Dirac gamma matrices. In Appendix B we list Feynman rules for the flowed vertices and for the relevant operators. In Appendix C we use the calculation of the quark propagator as an example to elucidate the computational techniques for finite flow time.

II. THE GRADIENT FLOW
In this section we give a brief introduction to the gradient flow, emphasizing the technical details needed for our perturbative expansion. The gradient flow equations define the evolution of the bulk gauge and fermion fields, B µ (x; t) and χ(x; t) respectively, as a function of the flow time, t [22,24]: where and the covariant derivatives are The bulk fields are related via Dirichlet boundary conditions to the boundary fields, that is, the integration variables of the functional integral defining the theory, through The generalized gauge-fixing terms proportional to α 0 remove some technical complications associated with perturbation theory [22][23][24]. The solutions of the flow equations for α 0 > 0 are related to the solutions at α 0 = 0 by a flow-time dependent gauge transformation. We work in Feynman gauge and take α 0 = 1 throughout this work. We solve the flow equations (1) and (2) in d-dimensions by casting them into the integral forms Here K µν (x; t) and J(x; t) are the heat kernels and the interaction terms are We can solve the integral form of the flow equations, Eqs. (10) and (11), by iteration, generating a tree expansion of the bulk fields in powers of the boundary fields. Bulk vertices are then connected by "flow lines", which are flow-time ordered and governed by the heat kernel. We give explicit expressions for the relevant Feynman rules in Appendix B.
In pure Yang-Mills theory, all correlation functions are finite at finite flow time, provided the boundary theory is renormalized [23]. Fermions, however, require an additional wave-function renormalization at finite flow time, generally denoted by Z χ [24]. The pole contribution to this additional fermionic wave-function renormalization first appeared in [24] and was reproduced in [42], through a next-to-leading-order perturbative calculation of χ(x; t)γ µ ← → D µ χ(x; t)) , and in [58], in the context of nonlocal Wilson-line operators. In Appendix C we calculate the finite contributions to this extra wave-function renormalization that, to our knowledge, have not appeared in the literature. The calculation in Appendix C also serves as a sample calculation with flowed fermions fields.
Once the fermions have been renormalized, composite operators composed of fields at finite flow time are therefore finite and all scale dependence carried by these operators can be related to the flow time. In particular, any potential power divergence in the cutoff of the theory is removed. At small flow times, a short flow-time expansion (SFTE) can be used to relate these composite operators to linear combinations of local renormalized operators at vanishing flow time. The SFTE is an operator product expansion in the neighborhood of vanishing flow time, with coefficients, calculable in perturbation theory, that carry the flow time dependence [73]. The SFTE provides a perturbative understanding of the way in which power divergences are removed and the form of the flow-time dependence for which the power divergences are traded.
On the lattice, correlation functions involving higher-dimension operators can be plagued by power-divergent mixings with lower-dimension operators. In large volume calculations, the only accessible energy scale is the inverse lattice spacing ∼ 1/a, so the regularization and renormalization of correlation functions may depend only on the lattice spacing. Disentangling the dual roles of the lattice spacing, as cutoff and as energy renormalization scale, can be arduous, particularly in the presence of power divergences, which must be removed non-perturbatively.
The gradient flow provides a workaround: the flow renders all operators finite, and, in the continuum limit, the scale of all flowed correlators is parameterized by the flow time, µ 2 ∝ 1/t. The SFTE then provides a method to extract renormalized operators evaluated at t = 0 from finite operators calculated on the lattice at finite flow time, t > 0. In other words, we calculate correlation functions of local operators at non-vanishing flow time and relate them to physical correlation functions of boundary operators via a SFTE. The challenge associated with the renormalization of the correlators at t = 0 is traded for the difficulty of determining the expansion coefficients in the SFTE. One advantage of the SFTE, however, is that we can perform the analysis in the continuum, thus avoiding spurious chiral-symmetry breaking effects. In addition, the SFTE connects operators at several values of the flow time in a gauge-invariant way. This is a significant advantage compared to standard techniques, based for example on RI-MOM schemes, where determining the coefficients of the power divergent terms requires a non-perturbative gauge-fixing procedure [74][75][76][77][78]. An alternative gauge-invariant way to study power divergences is to use coordinate space renormalization methods [18,[79][80][81], although this does not provide a continuous probe of the fields, in practice.
We consider our theory in continuum Euclidean 4-dimensional space-time. For some gauge-invariant and local operator O i (t) in an associative operator algebra with basis B, defined at flow time t, the SFTE is [73] where the label R denotes a renormalized operator. Here, the Wilson, or expansion, coefficients c ij (t) have absorbed all flow time dependence, and the SFTE connects renormalized operators in the bulk and on the physical boundary. The SFTE is valid only if all fields are renormalized and all operators appearing in the SFTE are evaluated in correlation functions at non-zero physical distances to avoid spurious and additional contact terms. If the renormalized operators at vanishing flow time do not share the symmetries of the flowed operator, their expansion coefficients vanish. More specifically, the form of the SFTE and the operators contributing to the SFTE are dictated by the symmetries of the regulated theory. Thus, if our regulator breaks certain symmetries, those symmetries cannot be used to classify all the operators (O j ) R (0) contributing to the right-hand side of the SFTE in Eq. (19). The leading contributions in the SFTE stem from the lowest dimension operators and the renormalization group equation satisfied by the expansion coefficients dictates their asymptotic behavior at short flow time. In general, OPE's are linear and gauge-independent, so we are free to study the expansion in an arbitrary correlation function. Hence we are able, with the appropriate choice of external probes, to study the SFTE termwise, order-by-order. Moreover, the Wilson coefficients of the SFTE are universal, that is, the coefficients are insensitive to our choice of external states. This universality ensures that, once the Wilson coefficients are determined using one particular choice of external state, the resulting coefficients can be used with any other choice of external state.

III. QUARK-CHROMO ELECTRIC DIPOLE MOMENT
The effects of BSM physics at high energies can generate a set of effective, dimension-six, CP-violating operators at the electroweak scale. The five-dimensional qCEDM operator, which induces the nEDM at low energies, arises from the effects of electroweak symmetry breaking on the CP-violating Gluon-Higgs-Fermion operator [82]. We define the bare qCEDM to be where is a generalization of σ µν γ 5 that preserves Hermiticity in d dimensions [77]. All operators in this paper carry an arbitrary normalization factor, to simplify comparison to other results; in this case, k C is a complex number normalizing O C . The calculation of a renormalized qCEDM matrix element on the lattice is plagued by the presence of mixing with the other CP-violating operators [77]. In particular the mixing with the lower-dimensional pseudoscalar density generates power divergences in the lattice spacing a. A second lower-dimensional operator that mixes with the qCEDM is the topological charge density (TCD) The chirality of the TCD, opposite to the qCEDM, ensures the mixing is proportional to the quark mass. Our perturbative results confirm our expectations for the form of the power divergence and the mass dependence of the pseudoscalar density and the TCD, respectively. We remark that, if the lattice QCD calculation is performed with chiral symmetry breaking terms in the lattice action, chirality no longer protects the mixing of the qCEDM and the TCD and therefore a linearly divergent term in the inverse lattice spacing 1/a can arise. Although other operators of the same dimension mix with the qCEDM, in this work we focus on the calculation of the SFTE coefficients of the lower dimensional pseudoscalar density and TCD operators. The SFTE for the qCEDM reads where we have retained only the lowest-dimension operators in the expansion. For the remainder of this section we will not include contributions from higher dimensional operators, such as the renormalized qCEDM itself. The study of this logarithmic mixing will be considered in future work. In perturbation theory it is possible to extract the lowest dimensional operator contributions by selecting appropriate external sources [83]. Working at O(g 2 ), we can extract the expansion coefficients of the pseudoscalar and TCD by selecting two-fermion and two-gluon sources, respectively. With these choices of external states, and at O(g 2 ) in perturbation theory, higher dimensional operators do not contribute to the expansion coefficients.

A. Mixing with the pseudoscalar density
We start by extracting the coefficient for the pseudoscalar density. Choosing a two-fermion external state, we define, for any operator O, the connected correlation functions [84] Γ ψOψ (x, y; t) = ψ(x)O(t)ψ(y) .
We may then distribute over Eq. (24), so that where R once again denotes a renormalized quantity. Expanding both the correlation functions and the Wilson coefficients in powers of the renormalized coupling g, we find where the first term in the expansion of the left hand side of Eq. (26) vanishes because the correlator Γ R ψO Cψ (x, y; t) has no tree-level contributions, that is, the first term in the expansion of this correlator is the one-loop contribution proportional to g 2 .
Equating terms order-by-order and neglecting higher dimensional operators we obtain, up to O(g 4 ), The TCD vanishes at tree-level with two external quarks, Γ ψOqψ (x, y; 0) = 0, and we obtain The tree-level of the pseudoscalar density with two external quarks does not vanish Γ To extract c Cq (t) at leading order, we choose an external state with two gluons and define in analogy to Eq. (25). Applying the methods and results from above, because the tree-level of the qCEDM with 2 external gluons vanishes, Γ AO C A (x, y; t) = 0, and c AOqA (x, y; 0) = 0, from which we deduce that the leading order of the expansion coefficient c We are now in a position to extract c CP (t). There are three one-loop graphs that contribute to the left-hand-side of Eq. (33a), which we show in Fig. 1), and the correlator on the right-hand side is simply the tree-level for the pseudoscalar density.
We calculate these graphs to all orders in the external momenta and flow time. The inclusion of the mass results in a particularly cumbersome asymptotic analysis that lies outside the scope of this paper; a nonzero external momentum is sufficient to regulate all infrared divergences. The Feynman rules and mathematical details can be found in Apps. (A, B, and C). Additional mathematical details can be found in ref. [85]. We expand in powers of the quark mass and There are symmetric counterparts for diagrams (a) and (b), so the sum of these contributions is where we have omitted higher order corrections in flow time and quark mass. The final expression for the expansion coefficient reads We confirm the general expectation, based on symmetry and dimensional considerations, that the dominant contribution to the SFTE of the qCEDM is the pseudoscalar density, which has a corresponding expansion coefficient that diverges linearly in flow time.

B. Mixing with the topological charge density
To calculate the expansion coefficient c Cq (t), following Eq. (33b), we need to calculate the one-loop contribution Γ (1) AO C A (x, y; t), stemming from the three Feynman diagrams shown in Fig. 2. The graphs displayed in both 2b and 2c vanish under the traces of the fermion loops, so we are again left to calculate a single Feynman graph. To calculate the d-dimensional traces over fermion loops one could employ the 't Hooft-Veltman-Breitenlohner-Maison (HVBM) scheme [86][87][88]. Our conventions and details on the HVBM scheme can be found in Appendix A. This is, however, only necessary when these calculations are performed by expanding near p 2 = 0. Starting at O(p 2 ), this removes an essential IR regulator, the momentum, and introduces spurious divergences. The correlators listed below have been calculated by applying a new method that includes all orders in momenta, so our results are IR safe. The flow further removes all UV divergences, and the diagrams are finite in four dimensions. Following the methods outlined for the pseudoscalar density, we obtain We therefore find and We again confirm, following general chiral symmetry considerations, that the expansion coefficient for the TCD has a logarithmic dependence on the flow time. Chiral symmetry enforces the presence of a quark mass factor multiplying the TCD and this factor arises naturally in our calculation. Then, at non-zero mass, the qCEDM behaves, to leading-order, as where the dots indicate contributions from renormalized higher-dimensional operators.

IV. WEINBERG OPERATOR
Among the higher dimensional CP-violating operators obtained by integrating out heavy quarks and Higgs bosons, there is a dimension six gluonic operator, Weinberg's three-gluon operator [68], The Weinberg operator could potentially generate a large contribution to the nucleon EDM because it is purely gluonic and therefore not supressed by any small quark mass factor or by a small CKM phase.
To determine the SFTE of the Weinberg operator we need to isolate the lower dimensional CP-violating operators with the same symmetry properties. In principle, the pseudoscalar density, multiplied by a mass factor, could As with the qCEDM operator, we do not consider the contributions of operators with the same dimension as the Weinberg operator. The operators that could potentially contribute to the SFTE of the Weinberg operator originate from terms proportional to mO C and the Weinberg operator itself. By choosing external states of two quarks or two gluons, we can ensure that the leading contributions appear only at higher order in the external scales, such as momentum and flow-time, or at higher order in the coupling.
Expanding the Weinberg operator at short flow time, in a manner similar to the qCEDM, we obtain where we have considered only operators contributing to the expansion coeffcient c W q (t). These considerations confirm that the expansion coeffcient contribution from the qCEDM to the SFTE of the Weinberg operator starts at O(g 2 ).
We choose two gauge bosons as the external state and expand in powers of the coupling, leading to Equating order-by-order in the coupling, we obtain AOqA (x, y; 0) + c AOqA (x, y; 0) .
Thus the leading contribution to the expansion coefficient c W q vanishes, c W q = 0. The next order in the coupling expansion reads which allows us to determine c W q (t) once we have determined the one-loop contribution Γ AO W A (x, y; t). There are, once again, three Feynman graphs that contribute, which we show in Fig. 3. There are a large number of equivalent permutations of the fields of the Weinberg operator, so to simplify our calculations we employ a relation valid for any alternating 2-tensor which slightly generalizes the corresponding relation with Minkowski metric [69,70]. This relation decouples the indices of A, so that the permutations of any fields that may be contained in A become well-defined permutations on the indices within the trace. It should be noted that this formula is available in d-dimensions, but upon evaluation we reproduce exactly the four-dimensional trace in the HVBM scheme, so it may only contract nontrivially with other four-dimensional structures. This leaves only those pieces of a dimensionally-regularized integral that take values in the four-dimensional subalgebra.
In the calculation of the correlators involving the Weinberg operator, the flow automatically regulates the UV modes of the the bulk gauge field, and the external momentum controls infrared divergences. Thus all integrals are finite in four dimensions. Inserting the field tensor G in place of A, we find a simple expression for the Weinberg operator conducive to perturbative calculations: The Feynman rules for this operator are derived in Appendix B. The calculation of the Feynman diagrams in Fig. 3 leads tõ The second diagram has no logarithmic divergence in the flow time; a kernel line appears in place of the gauge boson propagator, which generates two additional powers of the loop momentum. The third diagram vanishes, because two of the legs on the Weinberg operator are contracted, and the Weinberg operator is antisymmetric with respect to its fields. Summing these contributions and factoring out the tree-level structure for the TCD, we isolate the Weinberg operator's leading-order divergent behavior: Our calculation again confirms the expectation that the leading contribution to the SFTE of the Weinberg operator stems from the lowest dimensional operator with the same symmetry properties; the TCD generates the linear divergence of the Weinberg operator at short flow time.

V. SUMMARY AND CONCLUSIONS
The nucleon electric dipole moment (EDM) provides a unique opportunity to probe of sources of charge and parity (CP) violation in the Standard Model and beyond (BSM). BSM theories that contain complex CP-violating couplings can induce a non-vanishing EDM, and at low energies one can parametrize the effects of the BSM degrees of freedom through effective, higher-dimensional CP-violating operators.
We have calculated, at one loop in perturbation theory, selected Wilson coefficients of the short flow time expansion (SFTE) for two CP-violating operators: the quark-chromo EDM (qCEDM) and the Weinberg operator. We have studied the leading contributions generated by the pseudoscalar density and the topological charge density, and confirmed the general expectation that the lowest-dimensional operators generate the dominant contributions at short flow time.
For the qCEDM, the Wilson coefficient of the pseudoscalar density is proportional to the inverse of the flow time, 1/t, and we have calculated the corresponding coefficient. In addition, we have calculated the logarithmic contribution to the qCEDM proportional to the topological charge density. Our calculation confirms the general expectation that chiral symmetry forces the contribution of the topological charge density to be proportional to the quark mass.
For the Weinberg operator, the leading contribution, which is proportional to the inverse of the flow time, stems from the topological charge density. We have determined both the coefficient of this 1/t term and additional logarithmic terms.
Further, we have introduced a method of evaluation for flowed loop-integrals, which permits, in many applications, the calculation of correlation functions in a natural four-dimensional setting. We fully avoid artificial divergences related to the zero-momentum or zero-mass calculations, while latently allowing for the study of these correlation functions at any or all positive values of momentum or mass. This also sidesteps the various problems that arise in continuing the spacetime algebra to any arbitrary dimension. This is particularly useful for our considerations, since the source of potential technical difficulties, γ 5 , is pervasive in CP-odd calculations yet well-defined only in four dimensions.
Our calculation is intended to provide a new framework to study the ultraviolet behavior of CP-violating operators contributing to the electric dipole moment. Ideally, the Wilson coefficients should be determined non-perturbatively and work in this direction is in progress [89]. Alternative strategies to pursue the same goals have been recently proposed based on coordinate space methods [18] and the RI-MOM scheme [77,78]. The one-loop calculation of the linearly divergent coefficients is also of practical importance for the non-perturbative determination of the Wilson coefficient, by constraining the perturbative behavior at small values of the gauge coupling.
We consider this calculation a first step toward the non-perturbative renormalization of all CP-violating operators contributing to the EDM. The next steps in our program are the non-perturbative determination of the linear divergence in the Wilson coefficients and a perturbative analysis that includes higher dimensional operators and their corresponding Wilson coefficients.

ACKNOWLEDGMENTS
We thank the members of the SymLat collaboration, Jack Dragos, Jangho Kim, Thomas Luu, Giovanni Pederiva, and Jordy de Vries for very useful discussions and a most enjoyable collaboration. In particular, we thank Jordy for his valuable insight in discussions regarding the renormalization properties of CP-violating operators and a careful reading of this manuscript. First, we define the set of generators for the gauge group, SU (N ), to be traceless and skew-Hermitian, so that the algebra is defined by for the N 2 − 1 generators t a ∈ su(N ), and for structure constants f abc . For any representation ρ : SU(N ) → GL(C), the trace over any two generators provides a natural Killing form for su(N ), normalized by the Dynkin index, We now turn our attention to two particular representations, the fundamental (F ) and the adjoint (A) representations, which have dimensions N and N 2 − 1, respectively. In these cases, our Casimir elements are C 2 (F ) = (N 2 − 1)/(2N ) and C 2 (A) = N , so the Dynkin indices become T F = −1/2 and T A = −N . Further, we can obtain an explicit set of generators for the adjoint representation by defining Clearly this definition is traceless and skew-symmetric, and it is trivial to prove that f abc must be real. Moreover, the Jacobi identity for f abc implicitly satisfies (A1), so that the N 2 − 1 matrices defined above indeed generate SU (N ). This allows for quick computations of objects such as

Quantum Chromodynamics
We work in d dimensions with a Euclidean metric, taking the d → 4 limit at the end. For all momentum integrals, we adopt the shorthand notation where µ is the energy scale introduced in dimensional regularization. We also define Fourier transforms so that the factor of (2π) d appears only in the momentum space measure: All calculations are performed on a QCD background, so that for any local operator O, correlation functions are given by with the gauge-fixed Lagrangian The generators of SU (N ) were chosen to be skew-Hermitian, so the covariant derivative is simply when acting on objects in the fundamental representation, where the coupling has been absorbed in to the fields, A a µ . When acting on objects in the adjoint representation, it assumes the form (A10) Then the field strength-tensor is
By definition, inner products between the two subalgebras vanish: and the metric tensors have a trace equal to the dimension of the subspace to which they belong: In the absence of γ 5 , this simply reduces to the natural d-dimensional generalization of the Dirac algebra. With γ 5 , however, there are some complications. In four dimensions, γ 5 is completely characterized by three properties: from which we find that, in d-dimensions, Since this prohibits a smooth limit for d → 4, we conclude that one of the above properties must be sacrificed to continue analytically to an arbitrary dimension. Our choice, introduced by 't Hooft and Veltman and systematized by Breitenlohner and Maison, relaxes the first condition (A17a), so that γ 5 anticommutes with the four-dimensional subspace and commutes with the (d − 4)-dimensional subspace. Thus Furthermore, the trace in (A17c) is taken to be fundamental, and the Levi-Civita symbol µνρσ is strictly fourdimensional, containing no evanescent components. As such, it is best to algebraically reduce expressions containing µνρσ after the d → 4 limit is taken. As a form of dimensional regularization, this scheme is manifestly Lorentz invariant, so that the reduction of tensor integrals is fairly straightforward. Moreover, the HVBM scheme maintains algebraic consistency in our applications; we have at most one instance of γ 5 in any correlation function. Finally, to maintain Hermiticity in all dimensions, we generalize the "pseudo-tensor" σ µν γ 5 = i 2 [γ µ , γ ν ]γ 5 to [77,78] Note that the tilde here does not signify a four-dimensional object as in the HVBM scheme; rather it is an unfortunate artifact of the literature. This modified version is central to the calculation of any correlation functions including the quark chromo-electric dipole moment operator.

Appendix B: Feynman Rules
We adopt the standard Feynman rules for QCD in d Euclidean dimensions. Below we describe in more detail the Feynman rules for gauge bosons and fermions at non-vanishing flow time. Some Feynman rules for flowed fields, and similar details relevant to perturbative calculations, have appeared already in the literature [22-24, 40, 42, 43, 58-62]. To keep this paper self-contained and provide a future reference, we list all the Feynman rules for flowed fields that we have used in these calculations, along with the relevant vertices arising from our operators. We note that all vertices with n-interacting fields are defined with inward-directed momenta p 1 , . . . , p n and that, unless stated otherwise (see Sec. B 2), there is an implicit factor of (2π) d δ (d) (p 1 + · · · + p n ) that ensures momentum conservation.

Gradient Flow
The nonlinearity of the flow equations produces extra vertices, which must be included in perturbation theory. For bosons, the vertices X (n,0) appear in the solutions of the flow equation, where n is the number of gluon fields involved. These flow vertices must always be connected to a kernel line. Kernels, called so for their role as the integral kernel of the solution to the flow equation, appropriately carry the information within a bulk field to its higher-order corrections. Diagrammatically, a kernel line may initiated at any vertex at positive flow time, replacing a bulk field leg, and terminating at a flow vertex. Thus for any interaction involving bulk fields with some functional form ∆(t), we will have corrections starting at O(g 0 ) attached with a kernel line. Let Γ(s) represent the associated flow vertex and all relevant subsidiary interactions involving all attached bulk fields. Then, representing a bosonic kernel as a double curly line, we define the Feynman rule: whereK is the bosonic kernel. Observe that it collapses to a simple Gaussian in the "generalized Feynman gauge," α 0 = 1. For clarity, note also that the ordering of the structures Γ and ∆ above is only restricted by the ordering of the fermionic fields contained within them. Turning our attention to the vertices, we have 1 2 X (2,0) at first order: The fields radiating out of this and all other flow vertices are bulk fields at some positive flow time, which in Eq. B1 we denote as s, whereas the kernel is generated by a bulk field at a flow time that, in Eq. B1, we denote t. The second-order vertex is 1 6 X (3,0) : The factors of 1/n! are placed within the vertex rules above so that the kernel line has the same Feynman rule regardless of the flow vertex to which it is attached. There are no intrinsic higher-order vertices, but these vertices may be nested to the desired order, ensuring that proper symmetry factors are included. For example, in the calculation a two-point Green's function at positive flow time and at one-loop order, we must account for all combinations up to O(g 2 0 ). Both vertices will contribute, along with the (at least) second-order structure: or, pictorially: The second line in (B5) is simply the NLO contribution to either of the two fields attached to the vertex X (2,0) (p, q, −p− q) abc µνρ . The initial factor of 2 accounts for the symmetry in choosing which of the B fields to expand. Since both fields include the same nonlinear corrections, either may be expanded, so long as the result is summed over all of these redundancies.
Fermions have similar rules. The fermionic kernels, produce Feynman rules analogous to the bosonic kernel. Letting ∆ and Γ be defined as before, and representing the fermionic kernel line by a double straight line, we have where the first rule applies to the flow-time evolution of the χ field while the second rule to the χ field. The distinction between J andJ is purely formal; J acts from the left on χ, while andJ acts from the right on χ. In the same manner as the fermion propagator, the direction of the arrow indicates the flow of fermion number from χ to χ. Analogously to what happens for the gauge bosons, the flow equations for the fermion fields (12,18) can be solved in an iterative manner, generating higher-order vertices containing one fermion field and n gauge fields, Y (1,n) . The term linear in B in the fermion flow equation produces Y (1,1) : while the analogous term in the adjoint fermion flow equation producesȲ (1,1) : where the first diagram refers to the perturbative expansion of the χ field and the second to the expansion of the χ field. The vertex Y (1,2) is thoroughly simpler: Since this term is quadratic in B, there is no sign change with respect to the direction of fermion flow, andȲ (1,2) is identical to Y (1,2)

Operators
In this section we list the Feynman rules for the CP-violating operators. The Feynman rules are flow-time independent, but the fields connected to these vertices may be flowed. The Feynman rules arising specifically from the perturbative expansion of the flowed fields are described in the previous subsection; only the tree-level fields enter our operator Feynman rules.
There is some subtlety in the implementation of these operators in perturbative QCD. A naïve calculation of any correlator with an odd number of CP-violating operators will always vanish. This should actually be expected; all correlation functions are calculated within a QCD background, so there may be no expectation values that violate CP. We circumvent this problem by temporarily ignoring momentum conservation; equivalently, we calculate all such correlations functions pointwise in coordinate space, integrating the point of interaction for our CP-violating operators over all spacetime only after we subtract off the desired quantities [69,70,90]. If momentum were to be conserved throughout these calculations, all operators would project to zero momentum at the onset, and structures like αβµν p µ p ν would contract to zero identically, trivializing the entire calculation. This trick allows us to break translational symmetry, giving the in and out states different total momenta and subsequently different transformations under the Lorentz group. After identifying the Wilson coefficients, we dynamically restore the conservation of momentum by integrating over all spacetime. In so doing, we also restore the appropriate discrete symmetries. We are simply keeping track of the various structures that vanish perturbatively.

a. Topological Charge Density
Appendix C: Sample calculation: One-Loop Fermion Propagator In this appendix we discuss in some detail the one-loop calculation of the fermion propagator for flowed fermion fields. Results for the one-loop calculation of the flowed fermion propagator have appeared in the literature [24,42,58] with varying degree of detail. We use this calculation as an example to elucidate features of a one-loop calculation at non-vanishing flow time and to collect all the relevant tools for a perturbative calculation with flowed fermion fields. For a more complete discussion of flowed perturbative calculations, we refer to [85].
The fermion propagator S(x, y; t, s) = χ(y; s)χ(x; t) = p e ip(x−y) S(p; t, s) , can be expanded in powers of the bare coupling with a tree-level expression The one-loop corrections can be calculated evaluating the Feynman diagrams depicted in Eqs. (C5a-C5h). There are eight nontrivial contributions to the flowed fermion propagator, of which only five are topologically distinct [24]. The diagrams involving flow kernels present some new features compared to standard perturbative calculations in QCD. While the standard one-loop diagram in Eq. (C5a) has the usual structure with tree-level propagators on the external lines, the flowed diagrams cannot truncated as easily, because they occur with one or two external kernel lines. For this reason we write the decomposition of the fermion propagator as follows i,a (p; t)S (0) (p; 0, s) +S (0) (p; t, 0)Γ (2) i,b (p; s) + Γ (2) 5 (p; t, s) .

(C4)
The functions Γ (2) i,a (p; t) and Γ (2) i,b (p; s) correspond to the first-order expansions of the external fields χ(x; t) and χ(y; s), respectively, though they are otherwise all but formally identical. The contribution Γ 5 includes the first-order expansion of both external fields. We list the individual contributions from each Feynman diagram in Eqs. (C5a-C5h) together with their evaluation, ignoring external propagators for brevity: where R(m 2 0 /p 2 ) is a remainder that vanishes for m 2 0 p 2 . The calculation of the first diagram Σ 1 (p) is identical to the standard QCD quark self-energy with tree-level external quark propagators carrying the flow-time dependence. We regulate the divergent integral with dimensional regularization with d = 4 − 2 and > 0.
The next contribution, proportional to Γ 2,a (p, t), contains a flow kernel and vertex. Following the Feynman rules we outline in Appendix B it is straightforward to write In standard perturbation theory, the integrand would next be recast with Feynman parameterization, shifted, decomposed into scalar integrals, and brought to a spherically-symmetric form for integration in d dimensions. Specifically, the integrand must be isotropic, so that the (d − 1)-dimensional surface may be integrated separately from the radial portion. This luxury is not afforded to us, however, as in this case, the gluon propagator introduces an exponential factor, e −(2p·q)t , which is only linear in the momentum q. No Feynman parameterization and corresponding shift in the integration variable will fix this; the exponential is neither even nor odd. Our solution is to reparameterize the propagatorà la Schwinger and to study the MacLaurin series of the cross-term: (−2(u + z)) n n! p µ1 · · · p µn q µ1 · · · q µn , where the sum over all µ n is implied. The symmetry of this structure is now manifest; that is, terms of even n are even, and terms of odd n are odd. We now let m 0 → 0, so that Indeed, in the complete calculation of the flowed diagrams of Eqs. (C5b-C5h), the mass only contributes at O(t). This allows for a concise demonstration of the techniques used in this article. In general the kernel diagrams do not contribute to all orders in the same way as the standard QCD diagrams (C5a). The full renormalization requires a coalescence of four semi-independent resummations. For this reason we only consider the leading Og 2 0 ) corrections and how they affect the wave function renormalization of the flowed fields. The above integral employs the multi-index I n = (µ 1 , µ 2 , . . . , µ n ). Note that the multi-index above is a 2n-tuple, because we neglect the mass and therefore the only term remaining outside of the gluon propagator, iq 2 , is even, and we may drop all odd n through the reindexation n → 2n. We have also rearranged the order of integration. In order to justify this, we invoke Tonelli: if the four integrals (including the sum, an integral with respect to the counting measure) in (C8) converge in some order, then we are free to choose any order, since the full integrand is strictly non-negative, and all domains of integration are clearly measure spaces with σ-finite measures. With this in mind, we freely reorder the integrals, and impose a posteriori restrictions on the integrals as we discover them. The momentum integral may now be calculated. Due to Lorentz invariance, the only available structure with the total indicial symmetry of the q I2n is the appropriately normalized sum over all (2n − 1)!! products of n metric tensors, where the indices are distributed according to all possible pairings. For example, for n = 2, we find for some smooth function f . In general, we have q f (q 2 )q I2n = q f (q 2 )q I4 = 1 (d) n,2 S (2n) I2n q f (q 2 )(q 2 ) n , where (d) n,2 = 2 n Γ(d/2+n) is a Pochhammer k-symbol, and the tensor S (2n) is the generalization of the structure in (C9). Each σ i is a permutation of the set [2n] ⊂ N corresponding to one of the (2n − 1)!! partitions without ordering of [2n] into n two-element subsets. For clarity, inspect the indices in (C9); each term splits the set {1, 2, 3, 4} into two unordered pairs, but the pairings are never the same. Indeed, any permutation of the indices simply permutes the summands. Thus the commutativity under addition of the terms in S (2n) I2n reproduces exact symmetry of the product of vectors q I2n . Further, we integrate over the (d − 1)-sphere to isolate the radial integral: The radial part is a simple gamma function, and the momenta p I2n saturate S To renormalize the propagator, following ref. [24], we define the renormalized flowed fermion fields as so that the renormalized propagator reads S R (x, t; y, s) = Z χ S(x, t; y, s) .
If we impose the family of conditions we obtain Z χ · 1 − g 2 0 C 2 (F ) (4π) 2 3 + log (8πµ 2 ) 2 st + log 4πµ 2 p 2 − γ E + 1 s=t,p 2 =µ 2 =1/(8π √ st) Expanding Z χ in powers of the bare coupling we find We note that if we choose the MS scheme we obtain the same result already obtained in Ref. [24], and that pole contribution matches the results of [42,58]. The finite terms, which depend on the choice of renormalization condition, have not, to our knowledge, appeared in the literature.