One-loop matching of CP -odd four-quark operators to the gradient-ﬂow scheme

: The translation of experimental limits on the neutron electric dipole moment into constraints on heavy CP -violating physics beyond the Standard Model requires knowledge about non-perturbative matrix elements of eﬀective operators, which ideally should be computed in lattice QCD. However, this necessitates a matching calculation as an interface to the eﬀective ﬁeld theory framework, which is based on dimensional regularization and renormalization by minimal subtraction. We calculate the one-loop matching between the gradient-ﬂow and minimal-subtraction schemes for the CP -violating four-quark operators contributing to the neutron electric dipole moment. The gradient ﬂow is a modern regularization-independent scheme amenable to lattice computations that promises, e.g., better control over power divergences than traditional momentum-subtraction schemes. Our results extend previous work on dimension-ﬁve operators and provide a necessary ingredient for future lattice-QCD computations of the contribution of four-quark operators to the neutron electric dipole moment.


Introduction
The search for sources of CP violation beyond the Standard Model (SM) is primarily motivated by the baryon asymmetry of the universe and has resulted in a very active program adressing both leptonic [1][2][3][4] and hadronic electric dipole moments (EDMs), see Refs.[5,6] for reviews.The current best experimental bound on the neutron EDM (nEDM) [7] |d n | < 1.8 × 10 −26 e cm (90% C.L.) (1.1) was obtained by the nEDM collaboration at PSI.Further improvements in the experimental sensitivities are expected in the near future [8][9][10][11][12].The SM contribution to the nEDM due to the CP -violating phase in the CKM matrix is several orders of magnitude smaller than the current experimental bound [13][14][15][16].Therefore, the measurement of the nEDM is an interesting probe of CP violation beyond the SM.
Given the absence of clear direct signals of physics beyond the SM at the LHC, new particles need to be either very weakly coupled or very heavy, with masses well above the electroweak scale.In the second case, their indirect low-energy effects can be described in terms of effective field theories (EFTs), in particular the Standard Model EFT (SMEFT) above the electroweak scale [17,18] and the low-energy EFT (LEFT) below the weak scale [19].The matching of models for new physics to the SMEFT is currently being automated [20,21] and the complete renormalization-group equations (RGEs) and the matching of the SMEFT and LEFT have been derived at one loop [22][23][24][25][26], enabling a treatment that avoids large logarithms in each step of the calculation.
The EFT approach is ideal to obtain constraints on new physics at a high scale from low-energy precision observables, such as the nEDM.However, the calculation of the observable itself within the LEFT involves matrix elements of effective operators between neutron states, schematically1 where L i (µ) denotes renormalized LEFT Wilson coefficients, and the sum runs over all renormalized operators O MS i that are not excluded by symmetry principles: due to the nature of the strong interaction at low energies, the hadronic matrix elements are non-perturbative and all CP -odd and flavor-neutral operators contribute.Due to the running and mixing effects of the RGEs, it is not trivial to turn the experimental bound (1.1) into a strong constraint on heavy new physics: in order to avoid possible cancellations, the uncertainties on the non-perturbative operator matrix elements in Eq. (1.2) should be reduced to a level of 10 − 25% [6,27].In addition, disentangling different sources of CP violation will require not only the nEDM as a single observable, but rather a whole portfolio of experiments [5,6].
Ideally, the non-perturbative matrix elements of effective operators relevant for the nEDM should be obtained from lattice-QCD computations, which provide a first-principles approach with controlled systematic uncertainties, see Ref. [28] for a recent review.Since the EFT framework is based on dimensional regularization and renormalization by (modified) minimal subtraction (MS or MS), the EFT description of the observable (1.2) involves matrix elements of MS operators.Therefore, the use of lattice-QCD input necessitates a matching calculation to a different renormalization scheme.Traditional schemes amenable to lattice-QCD computations are momentum-subtraction schemes (MOM), and the matching between MOM and MS has been worked out for the dimensionfive operators contributing to the nEDM [29], as well as the dimension-six CP -odd three-gluon operator [30].A more modern scheme is provided by the gradient flow [31,32], which promises a better control of power divergences [33,34].Recently, the one-loop matching between the gradientflow and MS schemes was worked out at dimension five [35].In the present paper, we extend this work to dimension-six four-quark operators.The gradient-flow matching was previously performed for the left-chiral current-current operators arising in the Fermi theory of weak interactions [36,37], using the naive dimensional regularization (NDR) and dimensional-reduction schemes.Here, we consider instead the CP -odd and flavor-neutral four-quark operators that contribute to the nEDM.
The article is structured as follows: in Sect.2, we define our operator basis up to dimension six.In addition to the physical operators, we define the unphysical nuisance operators that appear in the off-shell matching calculation, as well as evanescent operators in two different schemes.In Sect. 3 we briefly review the gradient-flow formalism and comment on a modification of the flow equations in the presence of an electromagnetic field.In Sect.4, we discuss the short-flow-time expansion and present our results for the one-loop matching coefficients of the flowed four-quark operators to MS operators up to dimension six.We conclude in Sect. 5 and provide some details on conventions and Feynman rules in the appendices.

Operator basis
The indirect effect of heavy physics beyond the SM on observables below the electroweak scale is described by the LEFT with the Lagrangian where the sum runs over all operators with mass dimension d ≥ 5 that respect the SU (3) c × U (1) em gauge symmetry.Although dipole operators in the LEFT appear at dimension five, their contribution beyond the SM arises in the SMEFT only at dimension six due to SU (2) L gauge invariance.Hence, in a scenario of heavy new physics where the LEFT is matched to the SMEFT at the weak scale, a consistent treatment should involve both dimension-five and -six effects in the LEFT.We start by identifying the LEFT operators at the hadronic scale of a few GeV that can contribute to the neutron EDM up to dimension six.Having an accuracy goal for the hadronic matrix elements of about 10 − 25% in mind [6], we are not interested in higher-order QED corrections and we treat the photon as a static external field.At the hadronic scale, we consider either three or four active quark flavors, collected in a vector q = (u, d, s) or q = (u, d, s, c).Due to non-perturbative effects, any operator with the right symmetry properties has to be taken into account: the operators need to be CP -odd and flavor neutral.
We start by listing the physical operators in Sect.2.1.In a second step, we extend the operator basis to unphysical operators that appear in the loop calculation.In Sect.2.2, we list the nuisance operators, which vanish by the equations of motion (EOM).Evanescent operators related to dimensional regularization are discussed in Sect.2.3: their definition is part of the renormalization scheme.

Physical operators
From the LEFT operator basis up to dimension six as classified in Ref. [19], we select the CP -odd and flavor-neutral operators that can contribute to the nEDM at leading order in QED.We work in D = 4 − 2ε Euclidean space-time dimensions and largely follow the conventions of Ref. [35], see App. A. The Lagrangian of Euclidean QCD is given by where the covariant derivative includes the external electromagnetic field A µ .The field-strength tensors are defined by The quark-mass matrix is active quark flavors, respectively.The gauge-fixing and ghost terms are with the covariant derivative in the adjoint representation We include a trivial leptonic Lagrangian with D µ = ∂ µ + A µ , but we disregard dynamical photons, i.e., we approximate the full LEFT by the QCD and leptonic Lagrangian, supplemented by a tower of effective operators: Without dynamical photons, the leptonic interactions are restricted to the effective operators as well as the coupling to the external electromagnetic field.In the context of the nEDM we focus on CP -odd and flavor-neutral operators up to dimension six that involve quarks or gluons.At dimension three, there is the pseudoscalar density, or CP -odd mass term where p is a fixed quark-flavor index.It is always possible to switch to a basis where the mass matrix is real and diagonal, which removes the CP -odd mass term.In general, this field redefinition involves an anomalous axial rotation, which affects the only CP -odd operator at dimension four, the theta term, or topological charge density: where the dual field-strength tensor is G µν = 1 2 ϵ µναβ G αβ .At dimension five, there are the electric and chromo-electric dipole operators, where we again only keep the flavor-diagonal contributions relevant for the nEDM and we do not sum over the flavor index p.We are using the definition [35] σHV depending on the scheme for γ 5 , with At dimension six we encounter the CP -odd three-gluon operator which will be left for future studies [38], as well as a larger number of four-fermion operators, which are the focus of this article.Leptonic four-fermion operators only contribute at higher orders in α QED , hence we restrict ourselves to semileptonic and non-leptonic operators (baryon-and leptonnumber-violating operators do not contribute to the nEDM at dimension six).Schematically, they have the form where Γ 1,2 denote Dirac structures.The condition that the operators be flavor neutral means that in the case of semileptonic operators we are interested in flavor indices (p = r) ∧ (s = t), while in the case of four-quark operators we need to consider either (p = r) ∧ (s = t) or (p = t) ∧ (s = r).

Semileptonic operators
To leading order in α QED , the contribution of semileptonic operators to the nEDM can be written as where p and r are fixed flavor indices.The nEDM is determined by the terms linear in the momentum of the external photon.Due to Lorentz and gauge invariance, to leading order in QED the nEDM receives a contribution from the first term only in the case of semileptonic tensor operators with Γ 1 = σ µν and Γ 2 = σµν or vice versa.The non-perturbative hadronic matrix element is the matrix element of a tensor quark-bilinear operator.The matching of quark bilinears to the gradient-flow scheme is known [39] and we have reproduced these results.
The second contraction in Eq. (2.15) vanishes unless Γ 1 = 1 and hence Γ 2 = γ 5 .In this case, the semileptonic operator contributes as a renormalization of a CP -odd quark-mass term.Again, this contribution can be shifted into the theta term by an anomalous axial field redefinition and the calculation of the hadronic matrix element of the pseudoscalar density is equivalent to the one of the nEDM induced by the topological charge [40].

Four-quark operators
The non-redundant set of CP -odd four-quark operators that contribute to the nEDM has been identified previously in Ref. [41].We write them as O S1 p = (q p γ 5 q p )(q p q p ) , O S8 p = (q p γ 5 t a q p )(q p t a q p ) , O S1 pr = (q p γ 5 q p )(q r q r ) , p ̸ = r , O S8 pr = (q p γ 5 t a q p )(q r t a q r ) , p ̸ = r , (q p σµν q p )(q r σ µν q r ) + (q r σµν q r )(q p σ µν q p ) , p ̸ = r , (q p σµν t a q p )(q r σ µν t a q r ) + (q r σµν t a q r )(q p σ µν t a q p ) , p ̸ = r , where no implicit sums over flavor indices are performed.The tensor operators are symmetric in the two flavor indices, in contrast to the scalar operators.There are in total n q (3n q − 1) CP -odd and flavor-neutral four-quark operators, i.e., 24 operators for n q = 3 quark flavors or 44 for n q = 4 quark flavors.Tensor operators with identical quark flavors in the two bilinears have been reduced to scalar operators using the Fierz relations [41].In the parity basis, the relevant Fierz relations in where round and square brackets denote Dirac indices and the sign from the anticommutation of fermion fields is not included.Together with the SU (N c ) Fierz relation for the anti-Hermitian SU (N c ) generators t a , the four-quark operators in Eq. (2.16) can be related to LEFT operators in the chiral basis as follows: O T 1 pr = −4 (q Lp q Rp )(q Lr q Rr ) − (q Rp q Lp )(q Rr q Lr ) + 2 N c (q Lp q Rr )(q Lr q Rp ) − 2 N c (q Rp q Lr )(q Rr q Lp ) − 4(q Lp t a q Rr )(q Lr t a q Rp ) + 4(q Rp t a q Lr )(q Rr t a q Lp ) , O T 8 pr = −4 (q Lp t a q Rp )(q Lr t a q Rr ) − (q Rp t a q Lp )(q Rr t a q Lr ) where the chiral fields are defined by (2.20) However, the relations (2.19)only hold in D = 4 space-time dimensions and in order to do a proper matching to the LEFT operator basis of Ref. [19], one needs to take into account evanescent operators in the Fierz relations [42][43][44][45].In the present context of the nEDM we find it convenient to work in the parity basis given in Eq. (2.16).

Nuisance operators
Since we will perform an off-shell matching calculation, we will encounter not only physical operators, but in addition unphysical "nuisance operators."They can be split into two classes [46][47][48][49][50]: on the one hand, we need gauge-invariant operators that vanish by the classical EOM, known as class-IIa operators.On the other hand, additional gauge-variant nuisance operators appear as the solutions of the Ward-Slavnov-Taylor identities, known as class-IIb operators.They can be constructed as BRST variations of operators with ghost number −1.

Equation-of-motion operators
From the Euclidean QCD Lagrangian (2.2), we obtain the quark EOM as well as the gluon EOM where the index p runs over the quark flavors.Operators that are proportional to the classical EOM can be removed from the operator basis by field redefinitions that effectively shift their effects to higher orders in the EFT power counting.Furthermore, S-matrix elements of EOM operators vanish [49].However, in an off-shell matching of Green's functions unphysical operators appear as counterterm contributions.The gluon EOM is not relevant if we restrict our attention to operators contributing to the nEDM: a CP -odd operator involving D µ G µν needs to have the form and hence it vanishes for flavor-conserving operators with p = r.
For the quark EOM operators, it is convenient to define the fields where We work with the following set of EOM operators: where no implicit sum over flavor indices is understood.

Gauge-variant operators
The gauge-fixing term in Eq. (2.5), which is introduced in perturbation theory, breaks SU (3) c gauge symmetry down to BRST invariance.This implies that we encounter not only gauge-invariant counterterms, but also nuisance operators that are only BRST exact but not gauge invariant [46][47][48][49][50].These nuisance operators could be avoided by working with the background-field method.A backgroundfield formulation of the gradient flow has been established in Ref. [51], but in the present paper we will work with conventional R ξ gauge.Since we do not consider dynamical photons, we do not need to fix the QED gauge.We will work in a formalism with manifest U (1) em invariance, which however also requires flow equations that respect electromagnetic gauge invariance, as will be discussed in Sect.3.This allows us to include only class-IIb nuisance operators that are U (1) em gauge invariant.They have been classified in Ref. [30] in the context of the three-gluon operator.Imposing momentum conservation, we drop total-derivative operators and relabel the remaining relevant nuisance operators from Ref. [30] as follows: Because the four-quark operators are not singlets under chiral transformations, we find it convenient not to include explicit factors of the quark masses in the operators (2.26) and (2.27), apart from the quark EOM fields q E and qE .

Evanescent operators and γ 5 schemes
In addition to the physical and nuisance operators discussed in Sects.2.1 and 2.2, in dimensional regularization in D = 4 − 2ε space-time dimensions we encounter evanescent operators that disappear in four dimensions.The convention for evanescent operators affects the finite matching coefficients for the physical operators at one loop and hence is part of the scheme definition.In order to avoid a mixing of the unphysical evanescent sector into the physical sector, evanescent operators should be renormalized to have vanishing matrix elements in D = 4 dimensions [42][43][44].However, in our case the evanescent operators are generated only at one loop and their renormalization does not affect the matching coefficients of the physical operators at one-loop accuracy.The definition of the evanescent operators depends on the scheme chosen for γ 5 in dimensional regularization.The simplest scheme is NDR, which treats γ 5 to be anticommuting with all Dirac matrices in D space-time dimensions: (2.28) In connection with the NDR scheme, we replace σµν in the effective operators by σNDR µν defined in Eq. (2.12).As is well known, the NDR scheme leads to algebraic inconsistencies with γ 5 -odd fermion traces [52].We will use it only for the self-matching of the four-quark operators, since the calculation of the matching of flowed tensor operators O T 1 pr and O T 8 pr to lower-dimension operators leads to ill-defined Dirac traces.
We will use the original scheme by 't Hooft and Veltman (HV) [53,54] for the complete matching of the four-quark operators.In the HV scheme, we replace σµν in the effective operators by σHV µν defined in Eq. (2.12), where the indices of the Levi-Civita symbol only run over four dimensions.
In the HV scheme, global chiral symmetry is violated by the regulator.This implies that our renormalized MS operators in general do not fulfill the chiral Ward identities.The Ward identities could be restored by finite renormalizations, as discussed, e.g., in Ref. [29].However, in vector-like theories such as QCD, gauge symmetry remains unaffected by the regulator and the scheme remains consistent without symmetry-restoring finite counterterms.Here, we do not perform symmetryrestoring finite renormalizations: this corresponds to a scheme choice and all symmetry-breaking contributions will cancel in relations between observables.A comprehensive treatment of these finite renormalizations will be provided in Ref. [55].

Evanescent operators in the NDR scheme
In the calculation of loop integrals, we encounter higher tensor products of gamma matrices, which in D = 4 space-time dimensions can be reduced to the tensor structures of the physical four-quark operators.In D = 4 − 2ε dimensions, such relations do not hold but require the introduction of evanescent structures.For NDR, we choose the same scheme for evanescent four-fermion structures as Ref. [26], but we translate it from the chiral to the parity basis.Explicitly, we define where our evanescent structures are related to the ones of Ref. [26] by RL . (2.31) The evanescent structures depend on the parameters ãev , . . ., fev , which can be expressed in terms of the evanescent-scheme parameters of Ref. [26] as ãev = a ev , dev = 4f ev − 3d ev , ẽev = e ev , fev = 4 7 f ev + 3 7 d ev . (2.32) Note that the tensor operators are defined in Eq. (2.16) in a manifestly symmetric way, since in NDR the following structure is evanescent but non-vanishing in D dimensions: where σµν is interpreted as σNDR µν .Finally, the evanescent four-quark operators are defined by inserting the corresponding evanescent Dirac structures (2.31) into two antiquark-quark bilinears.The operators with two different flavors are given by E (2),1 pr = (1 + 2ã ev )ε (q p γ 5 q p )(q r q r ) − (q r γ 5 q r )(q p q p ) + 1 2 (q p σµν q p )(q r σ µν q r ) − (q r σµν q r )(q p σ µν q p ) , E (2),8 pr = (1 + 2ã ev )ε (q p γ 5 t a q p )(q r t a q r ) − (q r γ 5 t a q r )(q p t a q p ) + 1 2 (q p σµν t a q p )(q r σ µν t a q r ) − (q r σµν t a q r )(q p σ µν t a q p ) , E (4),1 pr = (q p γ µ γ ν γ λ γ σ γ 5 q p )(q r γ µ γ ν γ λ γ σ q r ) − 8(5 + 2 dev ε)(q p γ 5 q p )(q r q r ) − 8(3 − 14 fev ε)(q r γ 5 q r )(q p q p ) + 4(2 − ẽev ε) (q p σµν q p )(q r σ µν q r ) + (q r σµν q r )(q p σ µν q p ) , E (4),8 pr = (q p γ µ γ ν γ λ γ σ γ 5 t a q p )(q r γ µ γ ν γ λ γ σ t a q r ) − 8(5 + 2 dev ε)(q p γ 5 t a q p )(q r t a q r ) − 8(3 − 14 fev ε)(q r γ 5 t a q r )(q p t a q p ) + 4(2 − ẽev ε) (q p σµν t a q p )(q r σ µν t a q r ) + (q r σµν t a q r )(q p σ µν t a q p ) , where p ̸ = r.In the case of a single quark flavor, the evanescent structures (2.31) do not show up at one loop.As mentioned above, tensor operators with identical flavor indices are not included in the physical operator basis, since they can be reduced to scalar operators using the Fierz relations (2.17) and (2.18).Since the Fierz relations only hold in D = 4 space-time dimensions, they give rise to the following evanescent operators with single quark flavors: where

.36)
Fierz-evanescent operators would also show up in the one-loop basis change from Eq. (2.16) to the LEFT operators in chiral basis.
In the HV scheme, we use the same form of Fierz-evanescent operators as in the NDR scheme, where we replace in Eq. (2.35) σµν by σHV µν .For the matching to the lower-dimension operators, which we only perform in the HV scheme, we need additional evanescent quark-bilinear operators.A generic classification can be found in Ref. [30].We work with operators of the form where O F µ1...µn are built out of ∂ µ , G a µ , A µ and color structures.In this scheme, we can easily project to the non-evanescent sector by keeping the momenta and polarization vectors of external photons and gluons in D = 4 space-time dimensions.

Gradient flow
In this section, we briefly review the gradient-flow formalism as established in Refs.[31,32], largely following the conventions of Ref. [35].The gradient flow is a gauge-covariant D + 1-dimensional extension of D-dimensional Euclidean QCD.The additional dimension is called flow time t with mass dimension [t] = −2.The theory agrees with Euclidean QCD at the boundary t = 0. Instead of working with the D + 1-dimensional theory including Lagrange-multiplier fields [56], we can directly use D-dimensional Euclidean QCD and work with flowed quark and gluon fields, defined through the gradient-flow equations.
In the present context, the final goal is to obtain hadronic matrix elements of LEFT operators contributing to the nEDM from lattice QCD.In the first place, we need a non-perturbative and regularization-independent definition of renormalized effective operators, which will be provided by the gradient flow.In a second step, we establish the matching to operators renormalized in the MS scheme, using the short flow-time operator product expansion (SFTE), together with one-loop perturbation theory.Since we are interested in the nEDM, we require matrix elements including an external electromagnetic field, which in the LEFT can couple through the quark electromagnetic current from the renormalizable part of the Lagrangian, or directly through effective operators.Conventionally, the gradient flow is applied to the pure QCD sector of the theory.However, in this case the flow equations explicitly break U (1) em invariance.At dimension six, this leads to matching contributions to unphysical MS operators in addition to the ones discussed in Sect.2, which are not invariant under U (1) em .Here, we propose to avoid this problem by using modified quark flow equations that include the (static) external electromagnetic field.
The flowed gluon and quark fields are defined by the flow equations [31,32] and the boundary conditions at vanishing flow time In Eq. (3.1), α 0 denotes a gauge parameter.The gluon flow equation depends on the flowed field-strength tensor, defined by and the flowed covariant derivative, given by for a field X c (x, t) transforming in the adjoint representation of SU (3) c .The flowed covariant derivative appearing in the flow equations for quark and antiquark fields is given by where X(x, t) denotes a field transforming in the fundamental representation of SU (3) c × U (1) em .
Note that here we include a coupling of the quark fields to the external photon field in order to preserve local U (1) em gauge invariance.The external photon field itself remains unflowed.The gradient flow directly maps the gauge fields to smooth renormalized fields, provided that the D-dimensional theory has been renormalized [56], whereas quark fields require an additional multiplicative renormalization [32].At one loop in the MS scheme, the flowed bare (anti-)quark fields χ (0) , χ(0) are related to renormalized fields χ, χ by A regularization-independent renormalization condition for the quark fields is provided by [57,58] ⟨0| χ(x, t) where In dimensional regularization, the so-called ringed fields χ and χ differ from renormalized MS fields by a finite renormalization, given at one loop by [35,58] If the theory is expressed in terms of renormalized coupling and quark masses, operators built from the flowed gluon and renormalized flowed quark fields are automatically UV finite [56].Therefore, the gradient flow provides a non-perturbative and regularization-independent definition of renormalized operators.At short flow times, the flowed operators can be related to MS operators via the SFTE, see Sect.4.1.If one can find a scale where lattice-QCD computations are not too expensive, but where α s is still small enough, perturbation theory can be used to determine those coefficients of the SFTE that are not affected by power divergences.Typically, this is the case at an MS renormalization scale of μ ≈ 3 GeV, which is related via to an MS scale of µ ≈ 1.13 GeV, where γ E is the Euler constant.In the case of power divergences, a reliable determination of the matching coefficients must use non-perturbative methods [34].In contrast to the standard approach [59], power divergences show up as 1/t n singularities.Therefore, the gradient flow disentangles power divergences from the continuum limit, which can be taken for any fixed finite flow time.
The differential flow equations have the form of modified heat equations and can be written in integral form [56]. Upon rescaling the gauge field B µ → g 0 B µ , an expansion of the integral equations in powers of the bare coupling g 0 leads to additional Feynman rules compared to perturbative QCD.We use the same conventions as Ref. [35], apart from the new interaction terms involving the external electromagnetic field that emerge from the quark flow equations.The corresponding Feynman rules are listed in App.B.

Short flow-time expansion
In the following, we apply the SFTE to the case of the four-quark operators given in Eq. (2.16).We define flowed renormalized versions of the (physical) operators by replacing the quark fields with the ringed flowed quark fields: The SFTE is an operator-product expansion that relates flowed operators to MS renormalized operators, integrands before integration in all scales apart from the flow time t of the inserted operator.In this case, the flowed one-loop integrals become very simple single-scale integrals and both momentum and flow-time integrals can be performed with standard methods.No loop integrals with insertions of (unflowed) MS operators need to be calculated, since the expansion of the loops results in scaleless integrals, which vanish in dimensional regularization.
A large degree of automation allows us to perform the calculation efficiently: we use qgraf [61] for the generation of the Feynman diagrams and our own Mathematica routines for the evaluation of the diagrams, partially making use of FeynCalc [62][63][64].
We perform the whole calculation with generic gauge parameters ξ and α 0 : their cancellation in the final result for the matching coefficients provides a strong consistency check.

Results for the matching coefficients
In this section, we present the results of the one-loop matching calculation of flowed four-quark operators to MS operators up to dimension six.We explicitly show the contribution of the finite renormalization ζ χ that arises from the relation (3.8) between flowed MS fields and ringed fields.For the matching onto four-quark operators, we present results in both the NDR and HV schemes, whereas for the matching to the lower-dimension operators we only use the HV scheme.For a compact notation, we define 1 , in the NDR scheme, 0 , in the HV scheme.(4.5) The logarithmic dependence of the matching coefficients on the matching scale µ is predicted by the anomalous dimensions of the MS operators.Applying the change to the chiral basis, we have checked that our results are compatible with the anomalous dimensions of Ref. [25].

Matching coefficients of the pseudoscalar density
The coefficients of the pseudoscalar density in the matching equations are obtained by computing the insertion of flowed four-quark operators into antiquark-quark two-point functions.The diagram topologies are shown in Fig. 1.We find the following results in the HV scheme: The matching coefficients of the pseudoscalar density contain at least one mass insertion, which follows from chirality arguments.Since we are matching to a dimension-three operator, the coefficients still contain power divergences in the flow time, which require a non-perturbative subtraction.Therefore, the perturbative result can only be trusted in the case of the m 3 contributions.The parts of the coefficients containing a Kronecker δ pr only contribute in the case of single-flavor operators in Eq. (4.4): since there are no single-flavor tensor operators in our basis, there remains no matching contribution from tensor four-quark operators to the pseudoscalar density.We still list these contributions in Eq. (4.6), because this enables the reconstruction of matching coefficients for the case of general flavors: at one loop, a flowed four-quark operator O X prst (t) with generic flavor indices matches onto generic pseudoscalar operators O P pr with a coefficient proportional to δ st , as well as onto pseudoscalars O P pt and O P sr with identical coefficients times δ sr and δ pt , respectively.

Matching coefficients of dipole operators
The four-quark matching coefficients of the dimension-five dipole operators are obtained by inserting flowed four-quark operators into Green's functions with a quark and antiquark and either a gluon or an external photon, see Fig. 2. We only provide results for the HV scheme, since in the NDR scheme the insertion of tensor operators leads to problematic γ 5 -odd fermion traces.The qq two-point and qqG three-point functions do not fix the coefficients of all the nuisance operators, but they allow us to uniquely determine the coefficient of the physical dipole operator.We checked that in combination with the qqGG four-point function, we obtain a consistent set of matching equations that uniquely determines all coefficients.The topologies for the qqGG four-point function are shown in Fig. 3. 4The results for the matching coefficients of the chromo-electric dipole operator read Figure 3: Topologies relevant for the qqGG four-point function: we do not show crossed diagrams or diagrams with reversed fermion flow.In diagrams with gray blobs, some quark propagators can be replaced by flow lines, according to the gradient-flow Feynman rules: the gray blobs denote either QCD or flow vertices.In total, there are 76 diagrams.
whereas for the coefficients of the electric dipole operator, we obtain The coefficients of the dipole operators are always proportional to a mass insertion, which again follows from chirality.Therefore, these matching coefficients do not contain any power divergences.As in the case of the pseudoscalar density, we list contributions that only appear for single-flavor operators (and hence are of no relevance in the case of tensor operators), in order to enable the reconstruction of the case of generic flavors: as before, a flowed four-quark operator O X prst (t) matches onto generic dipole operators O E,CE pr with a coefficient proportional to δ st , onto dipole operators O E,CE st with a coefficient proportional to δ pr , as well as onto dipoles O E,CE pt and O E,CE sr with identical coefficients times δ sr and δ pt , respectively.We note that in the case of general flavors, a gluon EOM operator as in Eq. (2.23) is generated, leading to a penguin contribution to the four-fermion matching, which is absent for flavor-conserving operators relevant to the nEDM.

Matching coefficients of four-quark operators
The q p q r → q p q r four-point function allows us to extract the coefficients of the MS four-quark operators in the SFTE.In the matching between operators of equal mass dimension, the expansion in the IR scales amounts to setting external momenta and quark masses to zero.The matching of the four-point function involves 72 diagrams with insertions of the flowed operators, shown in Fig. 4 (in the case of scalar operators, there are two different insertions in each diagram).A subset of the diagrams can be related to the matching of flowed quark-bilinear operators, shown in Fig. 5.
Collecting the results for the NDR and HV schemes in one expression, we obtain for the matching coefficients in the SFTE of the flowed scalar singlet operator The matching coefficients in the SFTE of the flowed scalar octet operator are given by For the flowed tensor singlet operator, we find whereas for the coefficients in the SFTE of the flowed tensor octet operator, we obtain At one-loop order, the finite renormalization ζ χ only affects the diagonal matching coefficients, which are equal to 1 at tree level.Due to the SU (N c ) Fierz relation (2.18), our results only hold for the color gauge group SU (N c ): we use Eq.(A.3) to write the results in a compact form. 5

Anomalous axial rotations
As noted above, the pseudoscalar density (2.9) can be removed by an anomalous axial field redefinition.This amounts to replacing the pseudoscalar density by an EOM operator [29]  one finds the following relation of bare operators: N ′ p := qEp γ 5 q p + qp γ 5 q Ep = 2m p O P p + qp ← → / D γ 5 q p − ∂ µ (q p γµ γ 5 q p ) .(4.13) Our matching results are given in the basis specified in Sect.2, i.e., we do not include the EOM operator N ′ p .Instead, one could eliminate the pseudoscalar density in favor of the EOM operator N ′ p .
This would result in a shift in the evanescent anomaly operator qp ← → / D γ 5 q p , which is only visible in the HV scheme.The requirement that renormalized evanescent operators have vanishing matrix elements induces then a finite renormalization of the theta term.This is a manifestation of the anomaly in dimensional regularization, where the determinant of field redefinitions in the path integral is always trivial.
In the presence of CP -odd operators that also violate chiral symmetry, physical CP violation arises only if complex phases remain in the Lagrangian after vacuum alignment [65].Here, we assume that the dominant source of chiral symmetry breaking is given by non-vanishing mass matrices.A detailed discussion of vacuum alignment can be found in Ref. [29].

Impact of higher-order effects
In the following, we evaluate the matching coefficients numerically and estimate their perturbative uncertainty.The matching coefficients depend on the flow time t and the MS renormalization scale µ through logarithms log(8πµ 2 t), which are dictated by the anomalous dimensions of the operators.For a given flow time t 0 , the matching should be performed at an MS scale close to in order to avoid large logarithms in the matching.As pointed out in Sect.3, this typically corresponds to a reference scale of μ0 ≈ 3 GeV in the MS scheme, or µ 0 ≈ 1.13 GeV in the MS scheme.An additional scale dependence of the matching coefficients arises through the strong coupling constant α s (μ), which however is beyond the control of our one-loop calculation.As in Ref. [35], we use this residual scale dependence as an estimate of higher-order corrections to the matching coefficients.To this end, we evaluate the matching coefficients at µ 0 = 1.13 GeV and around t 0 in the range t ∈ [(1/4)t 0 , 4t 0 ].The results for the diagonal matching coefficients c X,X (t, µ 0 ) are shown in Figs. 6 and 7 for the HV and NDR schemes, respectively.The blue curves are for fixed α s = α s (μ 0 ), whereas for the red curves the coupling is evaluated at

.15)
The scale dependence of the coupling is determined with the two-loop QCD β-function [66][67][68][69][70] and the input value at the weak scale α s (M Z ) = 0.1179 [71].The maximal difference between blue and red curves in the considered range illustrates the residual scale dependence and gives a very rough estimate of higher-order corrections to the matching coefficients: we expect the one-loop contribution to the matching coefficients to have a relative perturbative uncertainty of This estimate is based solely on the logarithmic terms of higher order in α s that arise from the RG evolution of the coupling.For the matching coefficient, one is mostly interested in the finite contributions, as all logarithmic contributions vanish at the matching scale (4.14).Obtaining the actual higher-order matching corrections requires a proper two-loop matching calculation, which is beyond the scope of this paper.We observe that compared to the tree-level matching the one-loop contributions provide corrections of the order of 30% − 60% to the diagonal matching coefficients, while the off-diagonal matching coefficients of course only start at one loop.The numerical values for the four-quark matching coefficients are listed in Table 1, using the two-loop QCD β function.We note that in some cases the HV scheme leads to much larger off-diagonal contributions than the NDR scheme.It would be interesting to see if this effect can be traced back to the spurious breaking of chiral symmetry by the regulator [55].

Conclusions
In the present paper, we have calculated the one-loop matching for flavor-neutral CP -odd four-quark operators between the gradient-flow scheme and the more familiar MS scheme used in EFT analyses.The matching coefficients are obtained by inserting flowed four-quark operators into two-, three-, and four-point Green's functions and applying the method of regions to extract the coefficients of the MS operators.We provide the coefficients of four-quark operators both in the NDR and HV schemes.For the calculation of the coefficients of lower-dimension operators, we have only used the HV scheme in order to avoid problematic γ 5 -odd traces in NDR.The matching coefficients to lower-dimension operators are provided in a way that allows one to reconstruct the case of generic quark flavors.
For the matching to the dimension-five EDM operator, we have employed modified quark-flow equations that involve the static external electromagnetic field and manifestly respect U (1) em invariance.Working instead with pure QCD flow equations would result in matching contributions to unphysical dimension-six operators that are not U (1) em gauge invariant.
The gradient flow is a promising scheme for the treatment of CP -odd operators that contribute to the neutron EDM, providing a regularization-independent definition of renormalized operators.In lattice-QCD implementations, the gradient flow disentangles the power-divergent mixing with lowerdimension operators from the continuum limit, which can be taken for any fixed non-vanishing flow time.Our results extend previous work on the gradient-flow matching of dimension-five operators [35] and they provide a necessary ingredient for future lattice-QCD computations of the contribution of four-quark operators to the neutron EDM.With a forthcoming study of the CP -odd three-gluon operator [38], the gradient-flow matching for the complete set of operators up to dimension six that contribute to the neutron EDM will be available at one loop.As our studies show, in some cases the desired accuracy goal motivates the calculation of the matching at two loops [37,72].q 2 , β µ, q 1 q 4 , α ν, q 3 t = q 2 , β µ, q 1 q 4 , α ν, q 3 t = 2δ µν δ αβ ∞ 0 dt , (B.3) q 2 , β µ, q 1 q 4 , α a, ν, q 3 t = q 2 , β µ, q 1 q 4 , α a, ν, q

B.2 Operator insertions
Here, we list all the required vertex rules for effective operators.In contrast to the convention in Ref. [35], we regard the operators as part of the Lagrangian (2.8) and include in the Feynman rules both the Wilson coefficient and the minus sign from the exponential in the generating functional The rules for the quark electric and chromo-electric dipole operators are given by (all momenta q i are outgoing) r, j, β q 1 , µ p, i, α = −L E pr 2iδ αβ (σ µν ) ij q 1ν , (B. where Γ 1,2 are the Dirac structures of the two bilinears.For the color-singlet operators, the SU (3) c generators are replaced by the identity.We explicitly distinguish the two bilinears and hence in general require four instead of two separate insertions of the four-quark vertex into each topology, which doubles the number of diagrams mentioned in the main text.

Figure 1 :
Figure 1: Diagrams contributing to the four-quark matching to the pseudoscalar density.In the case of scalar operators, which are not symmetric under exchange of the quark bilinears, there are two different insertions for each diagram.

Figure 2 :
Figure2: Diagrams contributing to the four-quark matching to the dipole operators.In the case of scalar operators, there are two different insertions for each diagram.With our modified flow equations, the same diagrams appear in the photonic case, with the gluon replaced by a photon, whereas with pure QCD flow equations only the first and fourth diagrams would be present.

Figure 4 :Figure 5 :
Figure4: Topologies relevant for the qq → qq four-point function: in the first two diagrams, the dark-gray blob denotes the sum of all sub-diagrams contributing to the insertion of a flowed quark bilinear, shown in Fig.5.In the remaining diagrams, some quark propagators can be replaced by flow lines, according to the gradient-flow Feynman rules: the gray blobs denote either QCD or flow vertices.We do not show crossed diagrams.In total, there are 72 diagrams.

Figure 6 :
Figure6: Scale dependence of the diagonal four-quark matching coefficients in the HV scheme, evaluated at µ0.The solid blue curves are for a fixed value of αs(μ0) while the dashed red curves are for αs evaluated according to Eq. (4.15).Flowed operators are defined in terms of ringed fields.As input for αs we use the solution of the QCD β function at two loops.
. In the HV scheme,

Table 1 :
Numerical values of the four-quark matching coefficients in the HV and NDR schemes, evaluated at µ0, t0 and using αs(μ0)| 2-loop .The errors correspond to the (higher-order) scale dependence of the coupling for t ∈ [(1/4)t0, 4t0].The index i runs over the flowed operators (in terms of ringed fields), while j runs over the MS operators.