One-loop matching for quark dipole operators in a gradient-flow scheme

The quark chromoelectric dipole (qCEDM) operator is a CP-violating operator describing, at hadronic energies, beyond-the-standard-model contributions to the electric dipole moment of particles with nonzero spin. In this paper we define renormalized dipole operators in a regularization-independent scheme using the gradient flow, and we perform the matching at one loop in perturbation theory to renormalized operators of the same and lower dimension in the more familiar MS scheme. We also determine the matching coefficients for the quark chromomagnetic dipole operator (qCMDM), which contributes, for example, to matrix elements relevant to CP-violating and CP-conserving kaon decays. The calculation provides a basis for future lattice QCD computations of hadronic matrix elements of the qCEDM and qCMDM operators.


Introduction
The baryon asymmetry of the universe cannot be explained by known sources of charge (C) and parity (P) violation in the standard model of particle physics [1][2][3][4]. CP violation in the standard model occurs through the Cabibbo-Kobayashi-Maskawa (CKM) quark-mixing matrix, the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) neutrino mixing matrix, and, potentially, could be generated by the quantum chromodyamics (QCD) θ term. Electric dipole moments (EDMs) capture the distribution of positive and negative charge within different systems, and the existence of a permanent EDM in neutral particles, such as the neutron, can only occur in the presence of CP-violating interactions. Experimental constraints on the neutron EDM, d n = (0.0 ± 1.1 stat ± 0.2 sys ) × 10 −26 e cm at the 90% confidence level [5], leave open the possibility of unknown sources of CP violation several orders of magnitude larger than those generated by the standard model [6][7][8][9]. For a full review of EDM searches in a wide range of systems, see Ref. [10].
Under the assumption that these unknown sources of CP violation are due to heavy physics beyond the standard model (BSM), their indirect low-energy effects can be described in terms of effective field theories. Above the electroweak scale, this is the standard model effective field theory (SMEFT), provided that electroweak symmetry is linearly realized [11,12], while below the weak scale one should use the low-energy effective field theory (LEFT), which is invariant under the SU (3) c × U (1) em gauge group. The complete running and matching in these effective theories up to dimension six has been worked out to one-loop accuracy [13][14][15][16][17][18][19]. At hadronic scales, the CP-violating BSM effects are described by higher-dimensional operators within the LEFT. The impact of these operators on the neutron EDM must be disentangled from the standard model contributions generated by the CKM matrix and the QCD θ term. One therefore needs to determine individual contributions to the nucleon EDM and to CP-odd pion-nucleon and nucleonnucleon couplings, which determine nuclear EDM and Schiff moments, at hadronic energy scales. Experimental bounds on the neutron EDM are expected to improve by an order of magnitude with the next generation of experiments [20]. Current determinations of the relevant hadronic matrix elements are primarily based on model estimates from QCD sum rules [21][22][23][24][25][26] or from chiral perturbation theory [27][28][29][30][31], which have large, O(50% − 100%), uncertainties and preclude rigorous reduction of systematic uncertainties in theoretical predictions. First-principles calculations of the hadronic matrix elements with controlled uncertainties are required to exploit fully the improved bounds from upcoming experiments. Lattice QCD, in which QCD is formulated on a Euclidean hypercubic spacetime lattice, and correlation functions are determined stochastically, provides the best approach for first-principles calculations of QCD at hadronic energies.
Lattice calculations of the matrix elements relevant to the nucleon EDM are challenging for two reasons. First, in Euclidean space the QCD θ term is complex, which prevents the efficient use of Monte Carlo methods. This difficulty can be circumvented by expanding the Euclidean action in θ, justified by the current experimental bound on the neutron EDM, which implies θ ∼ 10 −10 . This approach is theoretically well-defined, but faces a challenging signal-to-noise problem associated with the insertion of the θ term in correlation functions on the lattice. This problem can be mitigated through signal-to-noise-improved ratios [32]. Second, renormalization of higher-dimensional composite operators on the lattice is difficult. For example, the quark chromoelectric dipole moment (qCEDM) operator mixes under renormalization with the lower-dimensional pseudoscalar density. On the lattice, this mixing is proportional to an inverse power of the lattice spacing and diverges in the continuum limit. This power-divergent mixing must be removed nonperturbatively to extract meaningful results [33]. For a recent review on lattice-QCD results for the nucleon EDM, see Ref. [34].
We proposed applying the gradient flow [35,36] to renormalize both the QCD θ term and the BSM CP-violating operators, such as the qCEDM [37,38] 1 . At finite flow time the qCEDM is multiplicatively renormalizable, but to relate this operator to physical matrix elements, one must perform an operator-product expansion at short flow times [40]. We studied the leadingorder short flow-time expansion of the qCEDM and the CP-odd three-gluon operator at one-loop in perturbation theory in Ref. [41] and nonperturbatively in Ref. [42], and our nonperturbative implementation of the hadronic matrix elements required for this program is ongoing [32,37,38,[43][44][45]. Alternative nonperturbative methods for higher-dimensional operators are regularizationindependent momentum-subtraction schemes, which have been defined for the qCEDM [46] and the three-gluon operator [47], as well as coordinate-space methods [48]. The mixing of these operators under renormalization was studied first in [49][50][51], then calculated at two loops for the qCEDM in [52] and at two and three loops for the three-gluon operator in [53].
Power-divergent mixing with lower-dimensional operators hampers the renormalization of the quark chromomagnetic dipole moment (qCMDM) as well [54]. The flavor-changing qCEDM and qCMDM operators describe low-energy effects of heavy SM and BSM particles on flavor observables, such as the CP-conserving long-distance contributions to K 0 − K 0 mixing [55], direct CP-violation in hyperon decays [56], / and the ∆I = 1/2, K → ππ transition [57], or the CP-violating part of the K → 3π decay [55]. Furthermore, matrix elements of the flavor-conserving qCMDM can be used to extract CP-odd pion-nucleon couplings [31,58,59], which contribute to nuclear Schiff moments. A first lattice determination of the matrix element relevant to K 0 − K 0 mixing has been obtained by ETMC using twisted mass fermions [60]. Here we propose to use the same strategy adopted for the qCEDM to renormalize the qCMDM operator with the gradient flow.
In [41] we determined at one loop in perturbation theory the leading contributions to the short flow-time expansion of the qCEDM, which are generated by the dimension-three pseudoscalar density operator and the dimension-four topological charge density. Here we extend this calculation to determine the dimension-five contributions to the short flow-time expansions of the qCEDM and the related qCMDM operator. We include the complete set of operators up to dimension five that mix with qCEDM and qCMDM operators and extract the corresponding short flow-time expansion coefficients at one-loop order in perturbation theory. These coefficients are necessary to relate the hadronic matrix elements of the qCEDM and qCMDM operators, determined nonperturbatively from lattice QCD, to their counterparts in the MS scheme, which provide inputs into the phenomenological analysis of nucleon EDM measurements.
We organize the rest of this paper as follows. In Sec. 2, we introduce the gradient flow and notation relevant for our discussion of the short flow-time expansion in Sec. 3. We then determine the matching coefficients to the MS scheme in Sec. 4. In Sec. 5 we discuss the scale dependence of the matching coefficients, and we summarize our results and conclusions in Sec. 6. In the Appendices A and B we list our conventions and Feynman rules.

Gradient flow
The gradient flow introduces an additional coordinate t of mass dimension [t] = −2, the so-called flow time (not to be confused with the Minkowski time coordinate-in the following, t refers to the flow time) [35,36]. Euclidean QCD (see App. A and B for our conventions) can be regarded as the boundary theory of a D + 1-dimensional field theory at t = 0 [61]. Integrating out suitable Lagrange-multiplier fields in the D + 1-dimensional action is equivalent to imposing the following gradient-flow equations on gauge fields, B µ , and quark fields, χ, χ, in D dimensions: together with the boundary conditions Here α 0 is a free (gauge) parameter, required for perturbative calculations. 2 We use the same symbol for the flowed field-strength tensor as for the field-strength tensor at zero flow time.
The (differential) flow equations of (2.1), together with boundary conditions, can be rewritten as integral equations [61] where the heat kernels are 5) and the interaction terms are given by The flow equations in integral form (2.4) can be solved iteratively, which corresponds to an expansion in powers of g 0 upon rescaling the gauge field B µ → g 0 B µ . This allows one to express the flowed (bulk) fields in terms of the fundamental fields at the boundary t = 0. From the expansion of the kernel, one obtains propagator-like structures called flow lines. The interaction terms in (2.6) induce interaction vertices with three and four fields, the flow vertices. Our conventions for the Feynman rules and diagrams are given in App. B.

Short flow-time expansion
In the following, we consider Green's functions with operator insertions at finite flow time t. The goal is to extract the relation between renormalized flowed operators and MS operators at zero flow time: This "short flow-time expansion" (SFTE) is an operator-product expansion (OPE) that is valid at small flow time t, where the hard scale is proportional to t −1/2 . We will take into account operators up to dimension five on both sides of the matching equation (3.1).
To extract the flow-time dependent coefficients c ij (t, µ), we consider insertions of the operators O R i (t) in off-shell amputated one-particle irreducible (1PI) Green's functions. We work in the massless limit and consider the matching to one-loop accuracy.
The relation between amputated Green's functions of the bare operators and renormalized Green's functions is schematically given by where we allow a generic number of n ψ external fermion, n ψ antifermion, and n G gauge fields at zero flow time. 3 Using standard procedures we renormalize each field with the corresponding Z −1/2 renormalization factor and denote with Z MS ij the matrix renormalizing O

Dirac algebra
We perform the calculation both in the 't Hooft-Veltman (HV) scheme [62,63] as well as in the scheme with anticommuting γ 5 (NDR). For the CP-odd operators, we define the following Dirac structures:σ where as usual σ µν = i 2 [γ µ , γ ν ]. As the Levi-Civita symbol is a purely four-dimensional object, in the HV schemeσ µν only contains four-dimensional components. In order to compare the HV and NDR schemes, we introduce δ HV = 1 in the HV scheme, 0 in the NDR scheme.

Quark-field renormalization
We express the SFTE in terms of operators in the MS scheme. Therefore, we renormalize the parameters of the boundary theory in the MS scheme and define the renormalized coupling as One can easily switch to the MS scheme by replacing the MS scale µ with the MS scaleμ according to While the flow equations regulate most of the UV singularities, the fermion fields (as well as the gauge coupling and quark mass) require renormalization. In the MS scheme, we define renormalized flowed fermion fields according to where the superscript (0) marks the bare fields and The quark-field renormalization leads to a counterterm contribution to the two-point function of renormalized quark fields ofS ct (p, where the propagator is defined in (B.9). Summing tree-level and one-loop diagrams leads to the following finite result for the two-point function in the massless limit: Note that the finite part of the two-point function depends on the gauge parameters, while the flowed quark-field renormalization Z χ in (4.6) is independent of ξ and α 0 . To make contact with lattice calculations, it is necessary to implement a renormalization scheme that is regularization independent. This is achieved by imposing the following regularizationindependent renormalization condition on a gauge-invariant composite operator (for one quark flavor) [64,65] 0|χ(x; t) In dimensional regularization, this implies an additional finite renormalization compared to MS. The "ringed fields" are related to the MS renormalized fields by The prefactors (8πt) ε/2 are introduced in dimensional regularization because the renormalization condition fixes the dimension of the fieldsχ andχ to be equal to 3/2 instead of (D − 1)/2. The next-to-leading order (NLO) contribution to the vacuum expectation value in (4.8) is obtained from vacuum two-loop diagrams [65], leading to the following finite renormalization ζ χ : We have performed the calculation of ζ χ for generic ξ and α 0 , confirming its gauge-parameter independence.

Expanding loops
The matching coefficients c ij (t, µ) only depend on the flow-time t and the MS renormalization scale µ and can be expanded in the renormalized MS coupling α s = g 2 /(4π) as The coefficients c (1) ij are independent of the soft scales-we include powers of the quark mass explicitly in the operators. When solving the matching equation for the coefficients c (1) ij , the non-analytic dependence on the soft scales cancels between the MS and flowed loop diagrams. Therefore, one can apply standard techniques for matching calculations [66][67][68] and expand the integrands of the loop integrals in all scales apart from the flow time t, before integration: although this alters the analytic structure of the loop integrals and distorts the infrared (IR) structure, these IR modifications drop out in the difference between MS and flowed loop diagrams. Expanding the loop integrals leads to scaleless integrals for the operator insertions on the LHS of (3.3), which vanish in dimensional regularization, hence ultraviolet (UV) and IR singularities of the expanded loops are identical. Insertions of the flowed operators are free from UV singularities (apart from the renormalization of the gauge coupling and the quark-field renormalization Z χ ). The IR singularities of the expanded loop integrals on the RHS of (3.3) exactly match the UV MS counterterms. The finite matching coefficients c ij can then be most easily obtained from the expanded integrals of insertions of the flowed operators, which are single-scale integrals and straightforward to calculate. Even the inclusion of generic gauge parameters ξ and α 0 does not lead to major complications in the calculation of the integrals, hence we perform all calculations for generic ξ and α 0 , which provides a useful check: the coefficients c ij of gauge-invariant operators in the SFTE need to be independent of ξ and α 0 .
When expanding the loop integrands in all soft scales before performing the loop integrals, one potential pitfall arises in the calculation due to the fact that the ringed fields are renormalized through the condition (4.8) and not in the MS scheme, as we will explain in the following.
In the dimensionally regularized theory, we can relate both the renormalized MS and flowed operators to the bare operators at zero flow time by where the relation for the flowed operators involves the short flow-time OPE and hence an infinite sum and where we assumed for notational simplicity that the operators contain in total n fermion and antifermion fields. We have introduced an arbitrary mass scaleμ, which compensates the mismatch of mass dimension between the bare operators and the flowed operators in terms of ringed fields. If we chooseμ = (8πt) −1/2 , at one loop the renormalization factors Z R ik can be written as For a generic choice ofμ, this changes to where we have included the O(α 0 s ) evanescent structure linear in ε. The matching coefficients follow from (4.12): which is the continuation of (3.1) to D dimensions. The coefficient c ik (t, µ,μ) depends onμ through evanescent terms. We are only interested in the limit → 0: which is indeed independent ofμ, but only if the linearly evanescent O(α 0 s ) term is correctly taken into account.
Including the finite renormalization and the factor (8πt) ε/2 that compensates the mass dimension of the ringed fields in the dimensionally regularized theory, the matching equation (3.3) reads We have kept all renormalization factors in the equation, so that both sides are UV finite. In particular, the factor (8πμ 2 t) −nε/2 multiplies a UV finite quantity. Its evanescent component matches the evanescent term in c ij (t, µ,μ). Therefore, the matching equation can be simplified to (4.18) At this stage, the integrands on both side of the equation can be expanded in the soft scales before performing the loop integrals. The loops on the left-hand side are transformed into scaleless integrals, where UV and IR divergences are identical. The same IR divergences are generated on both sides of the matching equation, so that they cancel when solving for c ij . However, since after the expansion both sides of the equation contain divergences, it is important that evanescent terms on both sides are consistently taken into account and that the factor (8πμ 2 t) −nε/2 is cancelled by the evanescent part of c ij (t, µ,μ). Indeed, we have confirmed explicitly that all IR poles do cancel on either side by resumming to all orders the contributions from the relevant soft scales.

Chromo-EDM
We define the flowed qCEDM operator in terms of renormalized (ringed) fields as We want to extract the coefficients of the short flow-time OPE up to dimension five: where the MS operators are minimally subtracted versions of the following operators: in terms of renormalized fields with F µν the field-strength tensor of the external U (1) gauge field. The dual gluonic field-strength tensor is defined as G µν = 1 2 µνρσ G ρσ . The five coefficients can be extracted by computing insertions of the flowed operator into suitable 1PI Green's functions (including wave-function renormalization).
The mixing with the pseudoscalar density is obtained from the 1PI matrix element with external quark-antiquark states (we can also insert momentum −q into the operator in order to obtain the mixing with total-derivative operators): In total, one obtains [41] c where c m 2 P is a new result. Due to chiral symmetry (for a discussion of chiral symmetry at finite flow time see Ref. [69]), mixing with the topological charge density requires an insertion of a mass factor, which we include in the operator O mθ . The SFTE coefficient can be extracted by calculating the Green's function where momentum −q is necessarily inserted into the operator [70]. The loop calculation leads to γ 5 -odd Dirac traces that are not well-defined in NDR. The result in the HV scheme reads [41] c mθ (t, µ) = i 4π 2 1 + log(8πµ 2 t) . (4.25) In order to extract the mixing with the qEDM operator, we consider a quark-antiquark matrix element with an external U (1) field: Only two diagrams contribute, shown in Fig. 1. As we are not interested in total derivative operators, we can send q → 0 and obtain for the diagrams expanded in all soft scales: Up to the contribution of equation-of-motion operators, this results in There is no tree-level contribution to this coefficient, so at the one-loop level it is immaterial whether the flowed fermion fields are renormalized in MS or through the ringed fermion renormalization condition in Eq. (4.8).
Finally, the coefficient of the MS qCEDM operator is obtained from the matrix element with external quark-antiquark-gluon states: The list of Feynman diagrams is shown in Fig. 2. There are 19 additional diagrams, which follow from crossing the quark-and antiquark legs and inverting the fermion-flow direction. We also do not show the diagrams needed for the calculation of the quark-field renormalization. For ξ = α 0 = 1, several diagrams are of second order in the soft scales and can be discarded, but they contribute in the calculation with generic gauge parameters. In total, we obtain for q = 0 (i.e., without additional momentum insertion into the operator) (4.30) Up to the contribution of equation-of-motion operators, this results in the following matching t t t Figure 3: Feynman diagrams for the matching calculation of the flowed qCMDM operator onto the identity. coefficient, including the finite renormalization imposed by (4.8): The divergences of the expanded flowed diagrams cancel in the matching equations against the counterterms on the MS side, which are determined by the anomalous dimension of the qCEDM operator. We have again checked that the result for the matching coefficient is independent of the gauge parameters ξ and α 0 .

Chromo-MDM
Similarly to the qCEDM operator, we also define the CP-even flowed qCMDM operator in terms of renormalized (ringed) fields as where the renormalized MS operators are the minimally subtracted versions of Due to chiral symmetry, mixing with the identity or with the gluon kinetic term requires an insertion of a mass factor, which we include in the definition of the operators. The mixing with the identity starts at two loops and is determined by the vacuum diagrams shown in Fig. 3: only the first diagram gives a non-vanishing contribution.
The remaining coefficients can be calculated in analogy to the CP-odd case. The results agree with the CP-odd sector for δ HV = 0: In the case of the qCMDM, flavor off-diagonal components are also of interest [54,60], because they can mediate BSM contributions to K → ππ and ε /ε. The matching coefficients c S (t, µ), c M (t, µ) and c CM (t, µ) are the same for flavor diagonal and off-diagonal components, while only the diagonal components contribute to c mG (t, µ), c m , c m 3 and c m 5 . For the flavor-changing components of the qCMDM, the factor m 2 in O m 2 S is replaced by m 2 → (m 2 i + m 2 f )/2, with i and f the flavors of initial-and final-state quarks. With this replacement, c m 2 S is unchanged.

Scale dependence and perturbative uncertainty
The SFTE connects renormalized operators at positive and vanishing flow time t. Operators at positive flow time defined in terms of the flowed gauge field and the ringed quark fields are independent of the renormalization scale µ, which implies that the scale dependence of the matching coefficients has to be cancelled by the renormalization-scale dependence of the renormalized MS operator at vanishing flow time, up to higher-order corrections in the perturbative expansion.
In Sec. 4, we have defined appropriately renormalized flowed operators and determined the matching coefficients in the SFTE. The matching coefficients depend on the matching scale and the flow time in the combination log(8πµ 2 t). The additional scale dependence of the coupling α s is beyond the accuracy of the one-loop matching calculation.
For illustration, we numerically evaluate the matching coefficient c CE (t, µ). We take as input the MS coupling at the weak scale α s (M 2 Z ) = 0.1179 [71] and evolve it down to an MS scalē µ = 3 GeV using the one-, two-, or three-loop QCD β-function [72,73]. This corresponds to an MS scale of µ 0 = 3 GeV × e γ E /2 (4π) 1/2 ≈ 1.13 GeV . (5. 2) The deviation from the tree-level result c CE = 1 illustrates the impact of the one-loop corrections, which are larger in the HV scheme than in NDR. The blue curves show the results for α s evaluated at the fixed scale µ 0 , while the red curves show the result for α s evaluated at a MS scale µ = (8πt) −1/2 . In principle, we could resum the leading logarithms (α s log(8πµ 2 t)) n in the matching coefficient by solving the renormalization-group equations. However, the logarithms in the matching are smallfor a matching performed at the flow time t 0 , they vanish. We checked that the resummation has a small impact in the range of t shown in Fig. 4. The mild logarithmic t-dependence of the blue curves simply reflects the scale dependence of the MS operators as dictated by the one-loop anomalous dimensions. The red curves lead to a result that differs from the blue curves by with constant coefficients A i and where the subleading logarithm dominates numerically and is beyond the accuracy of our calculation. We take the maximal difference between the blue and red curves in the range t ∈ [t 0 /4, 4t 0 ] as an estimate of genuine O(α 2 s ) corrections, which require a two-loop matching calculation: this points to a relative uncertainty at these scales of about ∼ 10%-20%.

Non-perturbative effects and scheme dependence
A comment is in order about perturbative estimates of power divergences in the flow time. We have calculated the matching coefficients including the ones involving lower-dimensional operators: the coefficients c P (t, µ), c m (t, µ), c m 3 (t, µ), and c S (t, µ) are power divergent for vanishing flow time. Similarly to what happens with power divergences with the lattice spacing [33], perturbative estimates of these matching coefficients are not sufficient to enable a reliable extraction of the matrix elements of the MS dipole operators. Non-perturbative effects leave unsubtracted divergent or finite terms. In Ref. [42], a first non-perturbative estimate of the power divergent mixing of the qCEDM operator with the pseudoscalar density has been provided, and a smooth transition, toward small coupling, to the perturbative result for c P (t, µ) ∼ t −1 obtained in Ref. [41] has been observed. The perturbative estimates of c P (t, µ), c m (t, µ), c m 3 (t, µ), and c S (t, µ) obtained in this work thus provide a useful guide, at small coupling, for a non-perturbative determination of the same coefficients.
A possible strategy to determine the renormalized MS dipole operators at vanishing flow time is to define the flowed dipole operators in terms of the ringed fields as in Eq. (4.19) and then to provide a non-perturbative lattice definition of the matching coefficient of the power-divergent term. In this way the matching coefficient will depend on the scheme defined by the ringed fields, as well as on the renormalization scheme used to renormalize the lower-dimensional operator at vanishing flow time. However, since this operator renormalizes multiplicatively, the dependence on the scheme at vanishing flow time cancels in the product of the matching coefficient and operator matrix element. Therefore, the scheme dependence of the power divergence is only through the renormalization scheme at finite flow time (e.g., the definition of the ringed fields), which can be chosen in a regularization-independent way and subtracted non-perturbatively. Once the non-perturbative subtraction of the power divergence is performed, in the continuum limit, the subtracted operator where the superscript S denotes the non-perturbative scheme at zero flow time (which can be defined even in a renormalization-group invariant way), will now have a SFTE with contributions only from operators of the same or higher dimensions. The expansion can then be studied using standard techniques. In particular the target qCEDM operator at vanishing flow time needs to be renormalized in the scheme used to determine the matching coefficients. Similar subtractions can be applied to the CP-even operators.

Summary and outlook
The next generation of experimental searches for a permanent nucleon electric dipole moment (EDM) are expected to improve the precision of current constraints by an order of magnitude or more. The neutron EDM offers a unique window into new sources of charge and parity (CP) violation generated by phenomena beyond the standard model of particle physics. The effects of heavy BSM particles on measurements at hadronic scales can be parametrized by effective, higher-dimensional CP-violating operators. In order to constrain the coefficients of these effective operators from data and ultimately obtain constraints on BSM theories, one needs to know the nonperturbative hadronic matrix elements of the effective operators. We have calculated the short flow-time expansion coefficients for the quark chromo-EDM and the quark chromo-MDM operators. We have determined the complete set of one-loop matching coefficients up to dimension five, extending the work on the leading coefficients in Ref. [41]. These finite one-loop matching corrections enter at next-to-leading-log accuracy. For the qCMDM operator, we also included the two-loop matching to the identity operator.
This work is part of ongoing efforts to calculate the hadronic matrix elements of CP-violating operators with lattice QCD by the SymLat collaboration [32,37,43,45] and others [74][75][76][77][78][79][80]. Our calculation provides the matching relations necessary to relate results extracted from lattice QCD to the MS (or MS) scheme, which is the one used in perturbative effective-field theory calculations. They are important to obtain robust constraints on BSM operator coefficients from the phenomenological analysis of experimental data. In addition, the perturbative calculation of the short flow-time coefficients will help to constrain the nonperturbative determination of the matching coefficient, first carried out in [42], by constraining it in the weak-coupling regime.

Acknowledgments
We thank V. Cirigliano and J. de Vries for useful discussions and collaboration at an early stage of the project. We further acknowledge several interesting discussions with T. Bhattacharya

A Conventions
We use traceless and anti-Hermitian SU (3) generators t a , where λ a are the Gell-Mann matrices. The generators fulfill

A.2 Dirac algebra and dimensional regularization
We use the Dirac algebra in D = 4 − 2ε Euclidean dimensions either in the HV or NDR scheme: The Dirac matrices are Hermitian, γ † µ = γ µ and fulfill In the HV scheme, we decompose the metric tensor into a part projecting onto 4 and −2ε dimensions, respectively, The projections of gamma matrices (or four-vectors in general) are defined as The fifth gamma matrix is defined as with the purely four-dimensional antisymmetric Levi-Civita tensor, 1234 = +1. The fifth gamma matrix is Hermitian, γ † 5 = γ 5 , it anticommutes with the four-dimensional gamma matrices, and it commutes with the gamma matrices in the D − 4-dimensional subspace: The HV scheme is algebraically consistent, and since QCD is a vector theory, spurious anomalies that arise in chiral gauge theories are absent in the present context (as long as we only consider single-operator insertions).

B Feynman rules
We start from QCD in Euclidean space in D = 4 − 2ε dimensions. We absorb the gauge couplings into the gauge fields and treat the electromagnetic field as an external static field. The Euclidean QCD Lagrangian is given by where the covariant derivative is when acting on fields which take values in the fundamental representation of SU (N ), or when acting on objects in the adjoint representation. The field-strength tensors are related to the commutator of the covariant derivative by [D µ , D ν ] = G µν + F µν , The Feynman rules are obtained from the generating functional