Off-lightcone Wilson-line operators in gradient flow

Off-lightcone Wilson-line operators are constructed using local operators connected by time-like or space-like Wilson lines, which ensure gauge invariance. Off-lightcone Wilson-line operators have broad applications in various contexts. For instance, space-like Wilson-line operators play a crucial role in determining quasi-distribution functions (quasi-PDFs), while time-like Wilson-line operators are essential for understanding quarkonium decay and production within the potential non-relativistic QCD (pNRQCD) framework. In this work, we establish a systematic approach for calculating the matching from the gradient-flow scheme to the $\overline{\rm MS}$ scheme in the limit of small flow time for off-lightcone Wilson-line operators. By employing the one-dimensional auxiliary-field formalism, we simplify the matching procedure, reducing it to the matching of local current operators. We provide one-loop level matching coefficients for these local current operators. For the case of hadronic matrix element related to the quark quasi-PDFs, we show at one-loop level that the finite flow time effect is very small as long as the flow radius is smaller than the physical distance $z$, which is usually satisfied in lattice gradient flow computations. Applications include lattice gradient flow computations of quark/gluon quasi-PDFs, gluonic correlators related to quarkonium decay and production in pNRQCD, and spin-dependent potentials in terms of chromoelectric and chromomagnetic field insertions into a Wilson loop.


Introduction
The path-ordered phase factor P exp ig A µ dx µ , referred to as the Wilson line 1 , is a fundamental object in gauge theories, playing a crucial role in connecting local fields to construct gauge-invariant operators.The study of the renormalization of Wilson lines and Wilson loops was initiated by Polyakov [1] and Gervais, Neveu [2].The all-order proof of the multiplicative renormalizability of Wilson lines was made using diagrammatic methods in refs.[3,4] and functional approach in ref. [5] based on the one-dimensional auxiliary-field formalism [2,[6][7][8][9]; see ref. [5] for a review.
Recently, there has been increasing interests in the study of the quasi-PDFs and pseudo-PDFs defined as the hadronic matrix elements of the Wilson-line operators of the type O Γ (zv) = ψ(zv)ΓW (zv, 0)ψ(0), O µναβ (zv) = g2 F µν (zv)W (zv, 0)F αβ (0), where ψ is a quark field, F µν is the gluon field strength tensor, Γ represents a Dirac structure, v is a four-vector satisfying v 2 = −1 (or v 2 = 1) for space-like (or time-like) Wilson-line operators in Minkowski spacetime 2 , z is a real number and W (zv, 0) is a straight-line path-ordered Wilson line connecting the two fields, in which the fundamental representation and the adjoint representation of the gauge group are assumed for the quark case, and the gluon case, respectively.Without loss of generality, we choose z > 0 throughout this work.Notably, quasi-PDFs and pseudo-PDFs can be factorized in relation to the standard light-cone parton distribution functions (PDFs) based on the short-distance operator product expansion (OPE) and the large momentum effective theory (LaMET) [21][22][23].Employing hadronic matrix elements of off-lightcone Wilson-line operators as opposed to lightlike ones, offers the advantage of being considerably more suitable to lattice computations.However, this convenience comes with the trade-off of the highly non-trivial renormalization of the off-lightcone Wilson-line operators.For a comprehensive overview, we refer to ref. [24].Furthermore, the vacuum expectation value of gluonic Wilson-line operators, such as the one defined in eq.(1.2), are related to the studies of quarkonium decay and production in the framework of pNRQCD, and the QCD vacuum structure (see ref. [25] for a review).
The gradient flow was originally introduced as a method to regularize and suppress the ultraviolet (UV) fluctuations [26,27] in lattice calculations.It has proven useful in lattice QCD scale settings and in lattice calculation of local matrix elements and correlation functions [28][29][30][31][32][33][34][35][36][37][38][39][40].Alongside its noise-reducing effect, gradient flow exhibits the notable property that composite operators at positive flow time remain finite even in the continue limit, provided they are expressed using renormalized physical parameters and fields [29].This property renders gradient flow suitable as a matching scheme between lattice and perturbative calculations.It is worth noting that gradient flow does not modify the infrared behavior of operators in the small flow-time limit, allowing the use of dimensional regularization to regulate infrared (IR) divergences when matching from the gradient-flow scheme to MS-like schemes.
The initial proposal to explore quasi-PDFs with the incorporation of gradient flow was presented in ref. [41].Subsequently, a one-loop matching calculation was performed, considering full flow time dependencies, for the hadronic matrix element defined using the operator in eq.(1.1) with zero external momenta, which is related to the quark quasi-PDF [42].In this study, we present a systematic approach for calculating the matching from the gradient-flow scheme to the MS scheme for Wilson-line operators, such as those defined in eq.(1.1) and eq.(1.2), specifically in the limit of small flow time.Our approach is based on the one-dimensional auxiliary-field formalism, and we provide comprehensive details of the one-loop order matching calculations.For more recent applications of perturbative gradient flow, see refs.[43][44][45][46][47][48][49][50][51][52][53] and the references therein.
We structure the remaining content of the paper as follows.In section 2, we introduce the one-dimensional auxiliary-field formalism and elaborate on its application in the renormalization of Wilson-line operators.In section 3, we introduce the gradient-flow formalism and implement the small flow-time expansion to the local current operators, yielding the core formulas for our matching calculations.Section 4 is dedicated to the computation of the matching coefficients up to one-loop order.Our summary and conclusions can be found in section 5.In appendix A, we present the detailed matching calculation with full flow time dependence for the hadronic matrix element that is related to quasi-PDFs at zero external momenta.We compare these findings with existing results from ref. [42], as well as with the results of section 4 in the small flow-time limit.
2 One-dimensional auxiliary-field formalism and renormalization of the off-lightcone Wilson-line operators

The one-dimensional auxiliary-field formalism
The one-dimensional auxiliary-field formalism is defined by enlarging the QCD Lagrangian to include an extra term [5], where the subscript "0" indicates a bare quantity, h v,0 is an auxiliary "heavy" Grassmann (or complex) scalar field in either fundamental or adjoint representations of the SU(3) gauge group, D µ 0 = ∂ µ − ig 0 A µ,a 0 T a is the covariant derivative and T a are the SU(3) generators in the appropriate representation.Note that, for v = (1, 0, 0, 0) and h v,0 being in the fundamental representation, the effective Lagrangian in eq.(2.1) essentially coincides with the Lagrangian of the heavy-quark effective theory (HQET) in the rest frame [54], except that we choose h v,0 to be a scalar.
The effective Lagrangian L hv is amenable to renormalization through the redefinitions, where Z hv is the wave-function renormalization constant of the "heavy" auxiliary field, Z g and Z A are the usual QCD renormalization constants of the strong coupling and the gluon field, respectively.There will be a "mass correction" iδm for the "heavy" auxiliary field and the effective Lagrangian L hv expressed in terms of the renormalized fields and parameters is given by3 , where proper representations are implied for the fields and the generator T a .
The "mass correction" iδm is linearly divergent.The linear divergence vanishes in dimensional regularization, but persists, in lattice regularization and the gradient-flow formalism.Furthermore, the "mass correction" iδm could incorporate contributions of O(Λ QCD ), attributed to the renormalon ambiguities that are recognized to be present in HQET [55,56].
Based on the one-dimensional auxiliary-field formalism, up to a constant normalization factor, the connected two-point function of the "heavy" auxiliary fields can be related to the Wilson lines through [5,24,57], where ⟨...⟩ hv stands for integrating out the "heavy" auxiliary field h v , and In this way, the Wilson-line operators can be substituted with products of local current operators.For the Wilson-line operators defined in eq.(1.1) and eq.(1.2), we have where the superscript "B" indicates bare composite operators, the bare "heavy" auxiliary fields h v,0 in O B Γ (zv), O µναβ,B (zv) are in fundamental representation and adjoint representation, respectively.

Renormalization of the off-lightcone Wilson-line operators
The operator as defined in eq.(1.2) is not renormalized multiplicatively because, in the presence of the external four-vector v µ , the components parallel and transverse to v µ are renormalized differently [57].We introduce the following projectors, to project out the parallel and the transverse components of the field strength tensor F µν by defining, (2.9) Notice that we ignore the mixing with gauge-noninvariant operators since they do not contribute to gauge-invariant observables.Nevertheless, we will verify the mixing patterns found in refs.[58][59][60] whenever they are relevant for our computations.Obviously, when v = (1, 0, 0, 0) and µ ̸ = ν, F µν ∥⊥ , F µν ⊥⊥ are proportional to the chromoelectric field E and the chromomagnetic field B, respectively.We introduce the following gluonic Wilson-line operators, where the superscript "R" stands for renormalized composite operators, µ is the renormalization scale, Λ represents a UV cutoff (UV regulator), δm(Λ) encodes the "mass correction" of the Wilson line, with its its linear divergence, while Z q (Λ, µ), Z ∥⊥ (Λ, µ) and Z ⊥⊥ (Λ, µ) encode all logarithmic divergences originating from wave-function and vertex renormalizations.
Within the one-dimensional auxiliary-field formalism, after subtracting linear divergences, the renormalization of the off-lightcone Wilson-line operators defined in eq.(1.1), eq.(2.10), and eq.(2.11) is simplified into the renormalization of two local "heavy-to-light" or "heavy-to-gluon" currents4 , ) (2.17) The renormalization formulas for the above gauge-invariant local current operators are given by, leading to the equivalence (2.23) In the following text, we may omit indicating the arguments Λ and µ in the renormalization constants.
The renormalization constants for the local current operators can be further decomposed as where Z V J , Z V |⊥ , and Z V ⊥⊥ represent the corresponding renormalization factors arising from vertex corrections.Implicit in this notation is the dependence of the wave-function renormalization constant of the "heavy" auxiliary field from the color representations.
The Wilson-line operators that define the spin-dependent potentials in terms of chromomagnetic field insertions into a Wilson loop [67,68] are also of interest in this study.Within the auxiliary-field formalism, each insertion of a chromomagnetic field into a Wilson loop can be related to the following current operator as, where v = (1, 0, 0, 0) and the subscript F indicates that the "heavy" auxiliary fields are in the fundamental representation of the SU(3) group.The renormalization formula for this operator is , where Z V F represents the renormalization factor stemming from the vertex corrections of the current hv (x)gF µν ⊥⊥ (x)h v (x).
3 Gradient flow and small flow-time expansion

The gradient-flow formalism
In the following, we work in d-dimensional Euclidean spacetime with d = 4 − 2ϵ.The gradient flow extends the d-dimensional QCD to a (d + 1)-dimensional field theory by introducing an additional coordinate t (not to be confused with the time coordinate in Minkowski spacetime) of mass dimension −2, which is called the flow time [28,29].In the gradient-flow formalism, the flow time dependent gluon field B a µ (x, t) and quark field χ i α (x, t) are determined by the flow equations [28][29][30], where κ is a gauge parameter that drops out in physical observables, and with the boundary conditions In perturbative gradient flow, one extends the regular QCD Feynman rules by adding flow time dependent exponentials to the propagators, introducing extra "flow lines" and "flow vertices" [29]; see ref. [43] for more details about the gradient flow Feynman rules and higher order calculation techniques.In the following perturbative calculation, we adopt the convenient choice κ = 1.
With renormalized physical parameters, the flowed field B µ does not need further renormalization [28,29], but the flowed quark fields χ need to be renormalized [30], where C F = (N 2 c − 1)/(2N c ) with N c = 3 the number of colors.In order to avoid the complication due to the matching between Z χ in dimensional regularization and that in the lattice regularization, the ringed fermion fields χ(x, t), χ(x, t) were introduced in ref. [35], with for one quark flavor, where The next-next-to-leading order (nnLO) result of Zχ can be found in ref. [43].At next-to-leading order (NLO), the renormalization factor Zχ is given by [35] where the factor (8πt) −ϵ is introduced to balance the dimension of the ringed fermion field, which has dimension 3/2, while the flowed fermion field has dimension (d − 1)/2.

The small flow-time expansion
In this work, we consider the local current operators defined in eqs.(2.15) and eq.(2.27), which are flowed at the flow time t.We use the same symbols for both the flowed and the un-flowed current operators as long as it does not lead to confusions.In the context of Wilson-line operators, we flow the external fields and the Wilson lines at the flow time t.Specifically, this means that all the gluon fields directly connecting to the "heavy" auxiliary field lines are flowed at the flow time t, while the "heavy" auxiliary fields themselves remain un-flowed.In the regime of small flow time, we employ a "short-distance" OPE for the flowed current operators.Our objective is to derive matching coefficients between renormalized flowed current operators and the MS renormalized un-flowed current operators, where we have neglected operator mixing.We work in the massless limit for quark fields and set all the external scales to be zero whenever feasible.This strategic choice allows us to leverage established techniques from HQET for matching calculations.Specifically, we do not need to calculate the un-flowed operators because they yield scaleless integrals to all orders in perturbative calculations in dimensional regularization so that we can obtain the matching coefficients c O (t, µ) from the renormalized flowed operators simply by subtracting the IR poles.
Similarly to the renormalization of current operators, we can break down the matching coefficients c O (t, µ) for the local current operators into products of field matching coefficients and corresponding vertex matching coefficients.For the four local current operators ψh v , gF µν ∥⊥ h v , gF µν ⊥⊥ h v , and g hv F µν ⊥⊥ h v , the relationships between matching coefficients for the local current operators and those for the corresponding fields and vertices are given by Here ζψ (t, µ), ζ F hv (ζ A hv ) and ζ A (t, µ) are the matching coefficients of the ringed fermion field χ, the "heavy" auxiliary field h v in fundamental (adjoint) representation and the flowed gluon field B, respectively.Additionally, ζ ψhv (t, µ), ζ ∥⊥ (t, µ), ζ ⊥⊥ (t, µ) and ζ F (t, µ) denote matching coefficients stemming from vertex corrections for the corresponding flowed local current operators.We calculate these matching coefficients to one-loop order in section 4.
A key insight arising from our study in section 2 on the multiplicative renormalization of Wilson-line operators is that matching coefficients for Wilson-line operators in the transition from the gradient-flow scheme to the MS scheme can be directly derived from matching coefficients for the corresponding local current operators.This occurs because in the one-dimensional auxiliary-field formalism, combinations of the renormalized local current operators at different spacetime points do not introduce additional UV divergences, apart from the subtracted linear divergence, and because of eq.(3.13) that relates local operators in the small flow-time limit.Therefore, additional matching for these combined local current operators in different spacetime points is not needed in the small flow-time limit.

Computations and results
While the matching coefficients for the fields can be gauge-dependent, it is important to note that the combinations leading to the local current operators under consideration in this work are gauge invariant.Therefore, we have chosen to employ the Feynman gauge consistently throughout our calculations.Our computations are conducted in the context of Euclidean spacetime.In Euclidean spacetime, the renormalized Lagrangian L hv is given by where the superscript "E" indicates Euclidean spacetime.Consequently, the Feynman rules for the "heavy" auxiliary field propagator, with momentum k flowing along the "heavy" auxiliary field line, and the interaction vertex are given by 1 k•v−i0 and −gv µ T a , respectively.The −i0 prescription in Euclidean spacetime plays key role in obtaining correct linear divergences and "mass correction" iδm.
We decompose a vector into parallel and transverse components according to (v 2 = 1 in Euclidean spacetime in our settings) We use the following Schwinger parametrization (α-representation) for the QCD propagator denominators.For the eikonal propagator, we use another version of Schwinger parametrization with k • v real.

δm & ζ hv
The calculation of the self-energy of the "heavy" auxiliary field at one-loop in gradient flow is similar to that of the heavy quark field in HQET.Setting the external momentum p flowing from left to right, the one-loop self-energy diagram shown in figure 1 gives the amplitude where μ2ϵ = µ 2 e γ E 4π ϵ and we have chosen the "heavy" auxiliary field to be in fundamental representation.For the adjoint representation case, we just need to replace C F with C A in the above expression.
Summing up geometric series of the one-loop self-energy diagrams of the "heavy" auxiliary field propagator gives 1/(p • v − M (a) ), which indicates that the one-loop "mass correction" is given by (−M (a) ).For simplicity, we set p = 0 in M (a) , yielding Note that the above expression is not odd under k • v → −k • v and the k • v integration does not give a vanishing result because the sign of the i0 prescription matters.
To derive the matching coefficient of the "heavy" auxiliary field at one-loop, we differentiate M (a) with respect to p • v, and set p = 0 for simplicity, resulting in Now, we explicitly evaluate eq.(4.6) and eq.(4.7) using the Schwinger parameterization in eq.(4.3) and eq.(4.4).We find (we set ϵ = 0 whenever the integral is finite) ) Performing the k momentum integration gives ) Integrating out x, y, we obtain the results ) Repeating the same calculation using dimensional regularization with t = 0 gives vanishing M (a), t=0 | p=0 and This confirms that we just need to subtract the IR poles from gradient flow results for the matching from the gradient-flow scheme to the MS scheme when external momenta are set to zero.Therefore, we obtain the following results (notice that the "mass correction" in the renormalized Lagrangian is denoted as iδm) ) where the C R = C F for fundamental representation case and C R = C A for adjoint representation case.The logarithmic dependence on t in the above result of ζ hv is consistent with the renormalization constant of h v calculated in [4].
In the following matching calculations, we drop the calculations of un-flowed diagrams since they can be easily obtained from the amplitudes of the corresponding flowed diagrams by setting t = 0 from the beginning, which lead to results proportional to 1 ϵ UV − 1 ϵ IR and obviously match the IR poles from flowed diagrams.We have explicitly checked that the coefficients of − log(2µ 2 te γ E ) in the matching coefficients calculated below and above align with those of the UV poles from the corresponding un-flowed operators at one-loop level , which can be served as a check point of our calculation.

ζψ
The Feynman diagrams for the one-loop corrections of the quark self-energy in the gradientflow formalism are displayed in figure 2. A similar calculation was performed in ref. [44], here we repeat the calculation for completeness, adhering to our strategy of splitting the matching coefficients into several parts that can be calculated separately.
Diagram (c 1 ) is the same as the one in ordinary QCD in the small flow-time limit, which gives the contribution For diagrams (c 2 ) and (c 3 ), we have and Diagram (c 4 ) and its corresponding mirror diagram contribute proportionally to the external fermion momentum (with the gradient flow gauge parameter κ = 1) and subsequently vanish upon setting the external momentum to zero.Diagram (c 5 ) yields a null contribution due to its inherent oddness under the gluon loop momentum transformation k → −k.
Adding up and M (c 3 ) , removing the UV divergences by introducing the ringed fermion field, and subtracting the IR divergences through matching procedure, we obtain the one-loop matching coefficient for the matching from the ringed fermion field χ to the MS renormalized fermion field ψ,

ζ ψhv
The one-loop correction Feynman diagrams for the ψh v vertex in gradient-flow formalism are given in figure 3. Setting the quark external momentum to be zero, we have the amplitude for diagram (b 1 ), Obviously, the transverse components of / k give vanishing loop integration, thus we have Diagram (b 2 ) gives contribution proportional to the external fermion momentum (with the gradient flow gauge parameter κ = 1) and subsequently vanishes upon setting the external momentum to zero.Therefore, we have the one-loop result of ζ ψhv after matching to the MS scheme,

ζ A
The Feynman diagrams for the one-loop corrections of gluon self-energy in the gradientflow formalism are given by figure 4. For simplicity, we set the external momentum to be zero while extracting the corresponding contribution.Diagram (d 1 ) is just the conventional gluon vacuum polarization, which gives the contribution in the small flow-time limit, Diagrams (d 2 ) and (d 3 ) give the contributions and from which, we see that the UV divergences in diagrams (d 2 ) and (d 3 ) cancel.
Diagram (d 4 ) gives the contribution Diagram (d 5 ) does not contribute in the small flow-time limit because it gives contribution proportional to t in the small flow-time expansion.The results of M (d 2 ) , M (d 3 ) and M (d 4 ) are consistent with those given in ref. [35].
Adding up M (d 1 ) , M (d 2 ) , M (d 3 ) and M (d 4 ) , we obtain with β 0 = 11 3 C A − 2 3 n f .The IR poles are removed through matching procedure and the UV poles are absorbed into the renormalization of the strong coupling in the MS scheme, which is consistent with the general argument that the flowed gluon field in B µ does not need renormalization while the strong coupling in B µ needs to be renormalized [29].Finally, we have the matching coefficient for the matching from the flowed B µ field to the MS renormalized gluon field gA µ , Note that the logarithmic dependence in ζ A cancel that in ζ hv (see eq. (4.16)) in adjoint representation at one-loop level.However, the constant terms in the square brackets in eq.(4.30) are not canceled.

ζ ⊥⊥ & ζ F
The one-loop vertex type Feynman diagrams for the operators gF µν ∥⊥ h v , gF µν ⊥⊥ h v and g hv F µν ⊥⊥ h v in the gradient-flow formalism are shown in figure 5. Choosing the external gluon momentum to be q (flowing into the vertex), we have the leading-order amplitudes where we have suppressed the color indices and dropped the external "heavy" fields but kept the external gluon fields for the convenience of matching procedure.
Since we are only interested in obtaining the matching coefficients, which are independent of the external momenta, we can set q • v = 0 for the amplitudes of the operators gF µν ⊥⊥ h v and g hv F µν ⊥⊥ h v , but not for the amplitudes of the operator gF µν ∥⊥ h v .The computations of the one-loop vertex corrections for the operators gF µν ⊥⊥ h v and g hv F µν ⊥⊥ h v are straightforward and similar.These three operators share the same diagrams (e 1 )−(e 5 ).We will only give the amplitudes of diagrams (e 1 ) − (e 5 ) for the operators gF µν ∥⊥ h v and gF µν ⊥⊥ h v , supplied with the amplitude of diagram (e 6 ) for the operator g hv F µν ⊥⊥ h v and the amplitude of diagram (e 7 ) for the operator g hv F µν ∥⊥ h v .The calculation of the vertex corrections for gF µν ∥⊥ h v is not so straightforward due to non-trivial mixing with the gauge non-invariant operator iv µ A ν ⊥ (iv • D − iδm)h v .We will deal with gF µν ∥⊥ h v in the next subsection in details.The amplitudes of diagrams (e 1 ) − (e 6 ) are given by Extracting the LO expressions and setting q = 0 thereafter, we obtain the contributions from each diagram for the operator gF µν ⊥⊥ h v , ζ ⊥⊥, (e 2 ) = 0, ( Summing up all the contributions from diagrams (e 1 ) − (e 5 ) and discarding the IR poles through matching procedure, we obtain the vertex matching coefficient for the operator gF µν ⊥⊥ h v , Similarly, we have the contributions from each diagram for the operator g hv ζ F, (e 2 ) = 0, (4.47) which lead to

ζ ∥⊥
Diagrams (e 4 ), (e 5 ) have no mixing and give the same contributions for gF µν ⊥⊥ h v and gF µν ∥⊥ h v , Complications arise from diagrams (e 1 ), (e 2 ), (e 3 ) and (e 7 ), whose amplitudes are given by These amplitudes do not give expressions that are exactly proportional to the LO expression of the operator, gF µν ∥⊥ h v even after tensor reduction for the loop integrations.They lead to mixing with the gauge non-invariant operators iv µ A ν ⊥ (iv • D − iδm)h v , whose LO expression is proportional to q µ ∥ A ν ⊥ , and thus it is not straightforward to disentangle the mixing from those that give contributions to the operator gF µν ∥⊥ h v .Fortunately, we can extract the contributions to the operator gF µν ∥⊥ h v through the results that are proportional to q ν ⊥ A µ ∥ .To check the mixing pattern and cancellation of linear divergences that were found in refs.[58][59][60], we will also keep the results that are proportional to √ t in the small flow-time limit.
If we set t = 0 in M ∥⊥ (e 1 ) and M ∥⊥ (e 7 ) , we can easily see that the linear divergences coming from M ∥⊥ (e 7 ) and the second term of M ∥⊥ (e 1 ) cancel with d = 3 − 2ϵ.However, in gradient flow scheme, the cancellation of linear divergences are much more complicated.Apparently, the linearly divergences arising from M ∥⊥ (e 7 ) and the second term of M ∥⊥ (e 1 ) do not cancel with d = 4 − 2ϵ and t > 0. The cancellation of linear divergences in gradient-flow scheme will be restored when diagrams (e 2 ) and (e 3 ) are included.
Great simplification can be achieved for our matching calculation by adopting the welldeveloped method used in the matching calculation in HQET.That is, we focus on the large loop momentum region k >> q and expand the amplitudes up to linear order of the external momentum q before performing the loop momentum integration.This method is valid in our case because, in the small flow-time limit, the flowed operators share the same IR behavior with the corresponding un-flowed operators and contributions to the matching coefficients come from the UV region.
Expanding the amplitudes according to k >> q up to linear order of q and completing the loop integrations, we obtain where From the above results, we can clearly see the cancellation of the linear divergences.After subtracting the contributions to the operator gF µν ∥⊥ h v , the leftover logarithmic dependence on t also matches the UV divergence of the mixing presented in refs.[58][59][60].
Summarizing all the contributions from diagrams (e 1 ) − (e 5 ) and (e 7 ), we obtain the matching coefficient contributed from the vertex correction of the operator gF µν ∥⊥ h v , which cancels the constant term in the square brackets of ζ A in eq.(4.30) and thus makes the one-loop contribution to the matching coefficient of the operator gF µν ∥⊥ h v vanish5 .

Results and applications for Wilson-line operators
We are now ready to summarize the final results of the matching coefficients for the four local current operators: ψh v , gF µν ∥⊥ h v , gF µν ⊥⊥ h v and g hv F µν ⊥⊥ h v .The relationships between the matching coefficients of the local current operators and the matching coefficients of the fields and vertices are provided in eqs.(3.14) We have verified that the logarithmic dependences on t in eq.(4.65), eq.(4.66) and eq.(4.67) align with the anomalous dimensions given by refs.[5,58,70] for the current operators ψh v , gF µν ∥⊥ h v and gF µν ⊥⊥ h v , respectively (see ref. [57] that uses similar notation for these operators for comparison).For the current operator g hv F µν ⊥⊥ h v , its anomalous dimension was also calculated in refs.[71,72] using background field method (see ref. [73] for a review of background field method), which is consistent with the logarithmic dependences on t in eq.(4.68).It is noteworthy that while gA does not contribute to the anomalous dimension of g hv F µν ⊥⊥ h v within the background field method, it does so using the method adopted in this work.The UV poles in the MS scheme for diagrams (e 1 ) and (e 4 ) differ depending on whether the calculations are conducted with or without the background field method, owing to the distinct Feynman rules of three-gluon vertex.We have checked that these differences compensate for those arising from the anomalous dimensions of gA with or without the background field method using Feynman gauge.With correct logarithmic dependences on t for the matching coefficients of local current operators, "heavy" auxiliary field h v , quark field ψ and gA, the corresponding logarithmic dependences on t for the vertex matching coefficients are ensured to be correct.While the results of c ⊥⊥ (t, µ) and c F (t, µ) coincide at one-loop level, it is intriguing that the anomalous dimensions of gF µν ⊥⊥ h v and g hv F µν ⊥⊥ h v differ by a term proportional to π 2 at two-loop level (see the two-loop anomalous dimensions given in refs.[57,74,75]).In addition, for the operator gF µν ∥⊥ h v , its anomalous dimension starts to contribute at two-loop level [57], the two-loop and three-loop anomalous dimensions of the operator ψh v can be found in refs.[76,77] and ref. [78], respectively.
As mentioned in the end of section 3.2, the matching coefficients for the Wilson-line operators can be directly obtained from the corresponding matching coefficients of the local current operators in the small flow-time limit within the one-dimensional auxiliaryfield formalism.In the small flow-time limit, expressing the matching equations for the Wilson-line operators defined in eq.(1.1), eq.(2.10) and eq.(2.11) as where we have omitted the Lorentz indices for simplicity, we establish relations These relations are independent of z and external states of the matrix elements constructed using these Wilson-line operators.
Regarding the spin-dependent potentials in terms of chromomagnetic field insertions into a Wilson loop, each such insertion will lead to a factor of the matching coefficient c F (t, µ) for the matching from the gradient-flow scheme to the MS scheme.

Summary and outlook
We have developed a systematic approach for matching from the gradient-flow scheme to the MS scheme for Wilson-line operators in the small flow-time limit.This approach extends to various applications, including quasi-PDFs, gluonic correlators that appear in quarkonium decay and production in the framework of pNRQCD and spin-dependent potentials in terms of chromomagnetic field insertions into a Wilson loop.
Within the one-dimensional auxiliary-field formalism, the matching from the gradientflow scheme to the MS scheme for the Wilson-line operators discussed in this work can be reduced to the matching for the corresponding local current operators in the small flow-time limit.We match these local current operators from the gradient-flow scheme to the MS scheme by comparing perturbative computations in these two schemes.Great simplification for the matching procedure can be achieved by using small flow-time expansion and setting external momenta to be zero whenever possible.In this way, the well-developed method of matching calculation in HQET can be applied.After completing the matching for the local current operators, one automatically obtains the matching coefficients for the corresponding Wilson-line operators in the small flow-time limit, because the combinations of the renormalized local current operators at different spacetime do not introduce further UV divergences, apart from the subtracted linear divergences, due to the multiplicative renormalizability of the Wilson-line operators.
The matching for the local current operators can be further decomposed as the matching for the fields and vertices due to the multiplicative renormalizability of the local current operators.We thus compute the matching coefficients at one-loop level for auxiliary "heavy" auxiliary field h v , massless quark field ψ, ψh v vertex, gluon field gA and the vertices of the local current operators gF µν ∥⊥ h v , gF µν ⊥⊥ h v and g hv F µν ⊥⊥ h v .The one-loop result of the linearly divergent "mass" correction δm is also obtained.Provided with these results, we give the matching coefficients for the local current operators and the Wilson-line operators in the small flow-time limit.
It is worth mentioning that this method is also applicable for matching from lattice regularization schemes to the MS scheme in the small lattice spacing limit, albeit with added complexity due to UV and IR divergence treatment [79].
In the appendix, we have shown that, at one-loop level, for the case of hadronic matrix element defined in eq.(A.1) at zero external momenta, the finite flow time effect is very small as long as the flow radius is smaller than physical distance z, which is usually satisfied in lattice gradient flow computations.
With the systematic approach we have developed for the matching from the gradientflow scheme to the MS scheme for the Wilson-line operators, it is straightforward to extend the current study to two-loop order, which would be more helpful for lattice gradient flow computations of matrix elements constructed using the off-lightcone Wilson-line operators.We leave it for future research.

Note added
Before submission of this work, we noticed that the one-loop result of the matching coefficient for the gluonic correlator G B was recently published in ref. [80], which was the same as our one-loop result of c 2 F (t, µ), because the matching calculation for the gluonic correlator G B could be reduced to the matching calculation for the local current operator g hv F µν ⊥⊥ h v .A Quark quasi-PDF operator in gradient flow In this appendix, we consider the hadronic matrix element constructed using the operator defined in eq.where P is the hadronic state with momentum P z , Γ = γ • v or Γ = γ α ⊥ .The one-loop matching calculation for h(z, P z ) from the gradient-flow scheme to the MS-like schemes was done in ref. [42] with P z = 0 and full flow time dependence.
However, the matching coefficient given by eq. ( 44) of ref. [42] does not lead to z independent result in the small flow-time limit, which contradicts with our argument that we only need to do matching calculations for the local current operators stemming from the Wilson-line operators based on the one-dimensional auxiliary-field formalism.Therefore, we re-calculate the matching coefficient in great detail and compare with the calculation done in ref. [42] 6 .
The matching relation for h(z, P z ) can be expressed as h R (z, t) = C q (t, z, µ)h MS (z, µ), (A.2) where we have dropped the exponentiated "mass" correction δm in the matching relation and choose P z = 0 for simplicity and comparison.The Feynman diagrams that contribute to the one-loop corrections for h(z, P z ) in gradient-flow formalism are shown in figure.where the subscript χ indicates that we have used the ringed fermion field to remove the UV divergences.Cq, t→0(t,µ)−1 .The renormalization scale µ is chosen so that log 2tµ 2 e γ E = 0.
To see how large the finite flow time effect will be, we plot the ratio of the one-loop corrections with full flow time dependence and the one-loop corrections in the small flowtime limit in figure.7.As it is shown in figure.7, the finite flow time effect is very small as long as the flow radius r F = √ 8t is smaller than the distance z (z > 1), which is usually satisfied in lattice gradient flow computations.

Figure 1 .
Figure 1.One-loop self-energy Feynman diagram for the "heavy" auxiliary field in gradient-flow formalism.The dashed arrow line denotes the "heavy" auxiliary field, the filled squares represent the flowed B µ field at flow time t.

Figure 2 .
Figure 2. One-loop quark self-energy Feynman diagrams in the gradient-flow formalism (mirror diagrams are implicit).The solid line is the fermion line, the double solid line is the flowed fermion line, the open circle represents the flow vertex.

(b 1 s(b 2 )Figure 3 .Figure 4 .
Figure 3. Feynman diagrams of the one-loop corrections for the ψh v vertex in gradient-flow formalism.The cross circle represents the ψh v vertex.The fermion field connected to the cross circle is flowed at flow time t.

(f 1 )(f 2 ) (f 3 )Figure 6 .
Figure 6.Feynman diagrams that contribute to the one-loop corrections for the hadronic matrix element defined as eq.(A.1) in gradient-flow formalism.The dashed line is the Wilson line, the filled squares represent the flowed B µ field at flow time t, the cross circles represent the spacetime points that connecting the Wilson line and the fermion fields, which are all flowed at the flow time t.Diagrams that give contribution proportional to the external momentum and thus vanish in our calculation are dropped.The external quark one-loop self-energy diagrams are not listed here, they are shown in figure. 2 (1.1),h(z, P z ) = ⟨P | ψ(zv)ΓW (zv, 0)ψ(0)|P ⟩, (A.1) 2 and figure.6.The contribution coming from the quark self-energy diagrams in figure. 2 can be summarized as 2 + γ E − log(432) , (A.3)

Figure 7 .
Figure 7.The ratio of the one-loop correction result with full flow time dependence and the oneloop correction result in the small flow-time limit: Ratio = Cq(t,µ,z)−1 .33) Figure 5. Vertex type Feynman diagrams that contribute to the one-loop matching for the operators gF µν ∥⊥ h v , gF µν ⊥⊥ h v and g hv F µν ⊥⊥ h v in the gradient-flow formalism (mirror diagrams are implicit).The gluon fields at the cross circles and the filled squares are flowed at flow time t.For g hv F µν ⊥⊥ h v , the fields hv are not drawn in diagrams (e 1 ) − (e 5 ) for simplicity.The "heavy" auxiliary fields in g hv F µν ⊥⊥ h v are in fundamental representation and the "heavy" auxiliary fields in gF µν ∥⊥ h v , gF µν ⊥⊥ h v are in adjoint representation.Diagrams (e 6 ) and (e 7 ) only give contributions to the operators g hv F µν ⊥⊥ h v and gF µν ∥⊥ h v , respectively.