Improvement, generalization, and scheme conversion of Wilson-line operators on the lattice in the auxiliary field approach

Nonlocal quark bilinear operators connected by link paths are used for studying parton distribution functions (PDFs) and transverse momentum-dependent PDFs of hadrons using lattice QCD. The nonlocality makes it difficult to understand the renormalization and improvement of these operators using standard methods. In previous work, we showed that by introducing an auxiliary field on the lattice, one can understand an on-axis Wilson-line operator as the product of two local operators in an extended theory. In this paper, we provide details about the calculation in perturbation theory of the factor for conversion from our lattice-suitable renormalization scheme to the MS scheme. Extending our work, we study Symanzik improvement of the extended theory to understand the pattern of discretization effects linear in the lattice spacing, a , which are present even if the lattice fermion action exactly preserves chiral symmetry. This provides a prospect for an eventual O ð a Þ improvement of lattice calculations of PDFs. We also generalize our approach to apply to Wilson lines along lattice diagonals and to piecewise-straight link paths.


I. INTRODUCTION
Calculating parton distribution functions (PDFs) of hadrons using lattice QCD is challenging. The most direct definition using a bilocal light cone operator is not accessible because lattice QCD is formulated in Euclidean space. The traditional approach is to compute Mellin moments of PDFs, which are obtained from matrix elements of local twist-two operators, but higher moments are problematic because of mixing with lower dimensional operators and increasing noise. In recent years, there has been a renewed interest in alternative approaches that use matrix elements of nonlocal operators that can be related to PDFs via perturbatively computable factorization formulas [1,2]. Of these, the most widely studied has been quasi-PDFs, proposed in Ref. [3], which use matrix elements of the nonlocal equal-time operator, O Γ ðx; ξ; nÞ ≡ψðx þ ξnÞΓWðx þ ξn; xÞψðxÞ; ð1Þ where ψ andψ are spatially separated by distance ξ in direction n and connected by a straight Wilson line W.
Some of the difficulties in the quasi-PDF approach arise from the use of a nonlocal operator. In Refs. [28,29], the auxiliary field approach [30,31] was used to represent the nonlocal operator as the product of two local operators in an extended theory. 1 This makes it possible to understand the renormalization properties of O Γ using standard techniques applied to the local operators. Specifically, the auxiliary field ζðξÞ, which is defined only along the line x þ ξn, is given the Lagrangian and O Γ is obtained using the local operator ϕ ≡ζψ: for m ¼ 0 and ξ > 0. Following the lattice theory for a static quark [35,36], in Ref. [28] we also defined a lattice action for the auxiliary field with n ¼ AEμ pointing along one of the lattice axes, ðx þ ξnÞ½∇ n þ m 0 ζðx þ ξnÞ; where ∇ and ∇ Ã are the forward and backward lattice covariant derivatives, respectively, and a is the lattice spacing. This enabled us to identify that the operator mixing observed in one-loop lattice perturbation theory [37] is caused by mixing between ϕ and the operator = nϕ, which is allowed when the lattice fermion action breaks chiral symmetry.
In addition, we presented in Ref. [28] a scheme for nonperturbative renormalization of O Γ , called RI-xMOM, which proceeds by renormalizing the auxiliary field action and the local composite operator ϕ. In Sec. II, we present the calculation using perturbation theory of the scheme conversion from RI-xMOM to MS, the result of which was reported in Ref. [28].
In Sec. III we supplement our previous work by applying the Symanzik improvement program [38] to analyze OðaÞ lattice artifacts. Finally, in Sec. IV we present generalizations of the auxiliary field approach to operators with gauge connections that are not straight lines and for off-axis gauge connections. Conclusions are presented in Sec. V. In Appendix A we relate our results on improvement to previous work done for the static quark theory and in Appendix B we compare with the improvement program based on the whole-operator approach for O Γ .

II. SCHEME CONVERSION
In this section we summarize the approach for nonperturbative renormalization that we introduced in Ref. [28] and provide details of the perturbative calculation of the conversion factor to the MS scheme. We stress that the results in Sec. III, which deals with the improvement of the auxiliary theory and of O Γ , and Sec. IV, which generalizes our approach to a broader range of operators, do not depend on the use of a particular renormalization scheme and that the auxiliary field framework can be used quite broadly for understanding Wilson-line operators.
The RI-xMOM scheme is based on a variant of the Rome-Southampton method [39] for nonperturbative renormalization. The essential feature is the definition of renormalization conditions that can be imposed both nonperturbatively on the lattice and in dimensionally regularized perturbation theory, which allows for conversion to the MS scheme. In Landau gauge, we make use of the positionspace bare ζ propagator, the momentum-space bare quark propagator, and the mixed-space bare Green's function for ψ: These renormalize as For the quark field renormalization, we adopt the standard RI 0 -MOM condition, The remaining conditions are imposed at momentum p 0 and distance ξ 0 , which define a family of renormalization schemes at scale μ 2 ¼ p 2 0 that depend on the dimensionless quantities y ≡ jp 0 jξ 0 and z ≡ p 0 · n=jp 0 j: Note that the last two conditions have been formulated to eliminate dependence on the linearly divergent counterterm m.
For one-loop conversion to the MS scheme, we work in d-dimensional Euclidean space with dimensional regularization and use general covariant gauge with gauge parameter λ. At one-loop order, S ζ and S ψ are depicted in Fig. 1.
We generically define conversion factors as C X ≡ Z MS X =Z RI X . In the case of the quark field, this has been computed e.g., in Ref. [40]: The free gluon propagator takes the form For S ζ , we will use it in position space, which is given in Ref. [31]: Together with the tree-level ζ propagator, the loop integral for the auxiliary-field propagator is straightforward. We obtain the following in MS at scale μ: where γ E is the Euler-Mascheroni constant. This agrees with the Oðα s Þ term in Ref. [41]. Equation (13) implies For the one-loop vertex function (Fig. 2), we use the free quark and gluon propagators in mixed space: position space parallel to n and momentum space for the d − 1 orthogonal dimensions. We decompose a general vector as v ≡ v n n þ v ⊥ , where v n ¼ n · v. The gluon propagator takes the formD AB μν ðx n ; p ⊥ Þ ≡ Z dp n 2π e ip n x n D AB μν ðpÞ where p Ã ≡ p ⊥ þ isgnðx n Þjp ⊥ jn. Similarly, for the quark propagator we obtain S tree ψ ðx n ; p ⊥ Þ ≡ Z dp n 2π e ip n x n −i= p We write the quantity in Eq. (14) as This is similar to an amputated vertex function, except that the ζ-leg "amputation" is done in position space rather than the usual momentum space. At one-loop order, this takes the form We simplify the calculation by restricting ourselves to the kinematics p ⊥ ¼ 0 and p n > 0 (i.e., z ¼ 1), which leads to Eventually, we obtain where Ci is the cosine integral function, CiðzÞ ≡ − R ∞ z cosðtÞ t dt. The ϕ conversion factor can be evaluated using with the appropriate kinematics. In Landau gauge (λ ¼ 0), we obtain We also find an anomalous dimension consistent with the leading term obtained first in the auxiliary field theory [30,31] and later also for the static-light current [41][42][43][44].
To match the auxiliary field mass, one can use the higherorder results from Ref. [41]. Comparing conventions, our S ζ corresponds to their coordinate-space iS r and our Z ζ corresponds to theirZ −1 Q . That reference gives where μ ξ ≡ 2 ξ e −γ E , α ξ ≡ α s ðμ ξ Þ, and f ζ ðxÞ is a series in x given up to Oðx 3 Þ. Together with the anomalous dimensioñ which is independent of the scale μ 0 . Perturbatively, this is limited byγ Q , which has been computed to order α 3 s [41,45].

III. IMPROVEMENT
In the Symanzik approach [38,46], the lattice theory is described by a continuum effective theory that is an expansion in powers of the lattice spacing, where each term has the same symmetries as the lattice theory. The same is also done for composite operators. The idea of improvement is to add higher dimensional operators in order to tune the parameters of the Symanzik theory such that the leading [e.g., OðaÞ] term in the continuum extrapolation of every correlation function is eliminated.
For the lattice fermion action we take Wilson twisted mass fermions, working in the twisted basis with the fermion doublet χ, although we will also consider chiral fermions. The leading continuum Lagrangian is twisted mass QCD, At the next order L ð1Þ χ contains terms that can be absorbed into the parameters of the leading Lagrangian along with one nontrivial one, the well-known "clover" term [46][47][48] associated with chiral symmetry breaking. For the auxiliary field in the continuum, we have Exact symmetries of the lattice theory include the little group of hypercubic rotations that preserve n and the Uð1Þ charge symmetry for the auxiliary field. There are also several spurionic symmetries [31,48,49]: (1) General hypercubic rotations, together with the rotation of n.
(2) Parity with respect to the axis n, i.e., the negation of the part of space-time vectors orthogonal to n, together with the negation of μ q . (3) Time reversal with respect to the axis n, together with the negation of μ q and n. (4) Charge conjugation, together with the negation of n.
Specifically, the auxiliary field transforms as ζ n ðxÞ →ζ −n ðxÞ T ;ζ n ðxÞ → ζ −n ðxÞ T : ð31Þ (5) Flavor SUð2Þ, together with a rotation of the twisted mass term μ q τ 3 → μ q e iα a τ a τ 3 e −iα a τ a . (6) For a chiral fermion action, SUð2Þ axial transformations together with a rotation of the total mass term m q þ iμ q γ 5 τ 3 . Following [50], the next-order auxiliary field Lagrangian has the form where b ζ and b 0 ζ are Oðg 4 0 Þ. These terms produce a quarkmass dependence in the auxiliary field mass counterterm.
This effect could be significant if charm quarks are included in the lattice action.

A. Local bilinear operator
To simplify the study of improvement, we consider a bilinear defined in the twisted basis, ϕ ¼ζχ. The operator = nϕ transforms in the same way under all of the above symmetries except for SUð2Þ axial. Therefore at leading order the renormalized operator is given by where r mix ¼ 0 for a chiral action.
For the OðaÞ contributions, we follow the study of the static-light axial current in Ref. [50]. However, because of the equation of motion ð= D þ m q þ iμ q τ 3 γ 5 Þχ ¼ 0, there is some freedom in defining the on-shell improvement terms. In particular, we could use a derivative orthogonal to n, = D ⊥ ¼ ðδ μν − n μ n ν Þγ μ D ν as in Ref. [50]. However, for the case of quasi-PDFs it may be better to use a derivative along n: this would keep the improved operator from extending across more than one time slice and could potentially allow OðaÞ improvement to be applied to existing data. Furthermore, we use the equations of motion 2 for ζ to writeζn · Dχ ¼ n · ∂ðζχÞ. At OðaÞ, we obtain the improved operator, For a chiral action,b ϕ ,b 0 ϕ , andc ϕ vanish and b ϕ ¼ b 0 ϕ but there still exist terms at OðaÞ that cannot be excluded. Note that this expression assumes only a single doublet of fermions is present; for additional nondegenerate fermions there can be additional terms involving, for example, the trace of the mass matrix, as discussed in Ref. [51]. Forφ ¼χζ, we apply charge conjugation and then relabel n → −n to obtain One further consideration is if the lattice gauge links used for the Wilson line are obtained using an anisotropic smearing, which breaks some of the hypercubic rotations and allows additional contributions. This is used in calculations by ETMC, where (for n in a spatial direction) the gauge links are smeared only in spatial directions and not in time [9]. There is again some freedom in defining the on-shell improvement due to the equations of motion for χ. If smearing is not performed in the t direction, then one possible form for the additional term is ðc 0 ϕ þc 0 ϕ = nÞγ tζ D t χ. We will not explicitly consider this anisotropic case below, but the generalization is straightforward.
For specific choices of Γ and τ, this expression simplifies. Calculations of quasi-PDFs are typically done with a flavor diagonal operator where τ commutes with τ 3 . It has become standard to compute unpolarized quasi-PDFs using Γ ¼ γ ν with ν satisfying n ν ¼ 0 and helicity quasi-PDFs using Γ ¼ = nγ 5 . In both of these cases, = nΓ= n ¼ −Γ and all of the (anti)commutators vanish, so that the improved operator becomes i.e., only the derivative operator (and no other mixing) contributes at OðaÞ and the only nontrivial effect of chiral symmetry breaking is the dependence on am q . Transversity quasi-PDFs are computed using an operator that has Γ ¼ = nγ ν in the physical basis. If μ q ¼ 0, then Eq. (39) again applies. On the other hand, when using twisted mass fermions at maximal twist we must consider the corresponding twisted-basis operator. For the case of a flavor diagonal transversity operator, the twisted-basis operator has Γ ¼ i= nγ ν γ 5 and suffers from equal-dimensional mixing because fΓ; = ng ¼ 2iγ ν γ 5 . Using a flavor-changing transversity operator would be advantageous: it appears the same in the physical and twisted bases and the form of the improved operator is given by Eq. (39).

Maximal twist
In calculations done using Wilson twisted mass fermions, many correlators benefit from automatic OðaÞ improvement [48,52]. This means that when tuned to maximal twist (m q ¼ 0), the OðaÞ contributions to those correlation functions vanish and there is no need to explicitly tune the improvement coefficients. For the case considered here, the arguments behind automatic improvement do not eliminate all of the OðaÞ contributions, but they do eliminate the contributions that vanish for chiral fermion actions.
We start by working at Oða 0 Þ and examining the equaldimensional mixing. To simplify the expressions, assume Using these two equations we can eliminate O fΓ;= ngτ and obtain We now consider the transformations  [8,9], where equal-dimensional mixing for the transversity operator was not explicitly treated. Now, we move on to the OðaÞ contributions. Specifically, we take the correlation function of the product of O Γ inserted at zero momentum with a renormalized OðaÞ-improved multilocal field Φ. The Symanzik expansion of the lattice correlator is given by where hÁ Á Ái 0 is evaluated in the continuum theory and O Γτ ðξÞ contains the OðaÞ terms that appear in Eq. (38). Assuming ΦO Γτ is invariant under R 1;2 5 and m q ¼ 0, many of the terms vanish and we get for some prefactors c 1;AE and c 2;AE . If the fermion action is OðaÞ improved by including a clover term, 3 then L ð1Þ χ vanishes and one can effectively obtain OðaÞ improvement by using where six parameters remain:Z ϕ;þ ,b 0 ϕ;þ ,c ϕ;þ ,Z ϕ;− ,b 0 ϕ;− , andc ϕ;− . Usually the operator will have a definite G Γ , so that only half of the parameters can contribute, and for many operators the commutator vanishes so that the term proportional to aμ q does not contribute. On the other hand, if the fermion action is not OðaÞ improved, then the second term in Eq. (45) can be eliminated either by explicitly treating the Oða 0 Þ mixing or by choosing an operator such that fΓ; = ng ¼ 0.
The term involving the OðaÞ part of the fermion action survives in Eq. (45) because the connection between dimensional counting and breaking of chiral symmetry, which underlies the usual arguments for automatic improvement, is broken by the equal-dimensional mixing with O fΓ;= ngτ . This leads to another potential worry: if Φ is not OðaÞ improved, as is typically the case for interpolating operators, then an additional nonvanishing term of the form ahΦ ð1Þ O R fΓ;= ngτ i 0 can appear in Eq. (45). However, although this is an OðaÞ contribution in the correlation function, it can be argued that it will not contribute to the ground-state hadronic matrix element determined at large Euclidean time separations, since the latter is independent of the interpolating operator. To see this, consider a correlation function using local hadron interpolators that have been explicitly OðaÞ improved. In this case, the Symanzik expansion is given by Eq. (45) and the OðaÞ terms can be eliminated by including a clover term in the fermion action and tuning the improvement parameters in Eq. (46). Once the correlation function is free of OðaÞ effects, then, following the arguments in Sec. 3.2 of Ref. [52], we find that the matrix element obtained at large time separations will be free of OðaÞ effects. As the same ground-state matrix element can be obtained using any interpolating operator, the OðaÞimproved ground-state matrix element can thus be obtained from the large-time-separation limit of a correlation function even if the interpolator is not OðaÞ improved.

C. Determining improvement coefficients
In principle, parameters associated with breaking of chiral symmetry (r mix ,b ϕ ,b 0 ϕ , b ϕ − b 0 ϕ , andc ϕ ) can be determined using improvement conditions derived from chiral Ward identities [46]; this approach was used for nonperturbative determinations of closely related parameters in the static quark theory in Refs. [56][57][58]. For the remaining improvement coefficients b ϕ and c ϕ , the situation is more difficult as there is no simple continuum physics condition to match onto: one would have to numerically study the continuum extrapolation of suitably chosen observables and tune the parameters to eliminate the linear dependence on a. Alternatively, the parameters could be computed in lattice perturbation theory, which has been done for the static quark theory using a few different lattice actions [50,[58][59][60][61][62][63]. 4

IV. GENERAL LINK PATHS
In this section we generalize our use of the auxiliary field approach on the lattice to describe paths that have corners and paths that are not along a lattice axis. In part, this will serve to understand under which circumstances the assumption made in Refs. [64][65][66][67] (and studied empirically in Ref. [68]), that the continuum renormalization pattern [31] applies to lattice calculations, is valid.

A. Piecewise straight paths
Nonlocal operators with link paths that are not straight are also used for hadron structure; in particular, stapleshaped gauge connections have been used for studying transverse momentum-dependent (TMD) PDFs. If the path is made from a finite number of segments, each of which is a Wilson line propagating along a lattice axis, then it is straightforward to accommodate the nonlocal operator in the lattice auxiliary field framework. An auxiliary field ζ n must be introduced for each segment, where n is the corresponding direction. In addition to suitable bilinears ϕ n i ≡ζ n i ψ andφ n f for the end points, one must also 3 Note that the more recent twisted mass lattice ensembles generated by ETMC, including those at the physical pion mass, do include a clover term [53][54][55]. 4 The parameter r mix can also be determined from the perturbative study of O Γ ðξÞ in Ref. [37]. For the gauge actions common to both calculations, that reference agrees with Ref. [60].
To be concrete, consider the TMD operator with quark fields separated by the four-vector b and the staple of height η in the direction of the unit vector v that is orthogonal to b (see Fig. 3): We introduce the auxiliary fields ζ v , ζ −v , and ζ −b , and obtain In addition to renormalizing the action for the auxiliary fields and the bilinears ϕ, the cusp operators must also be renormalized. The latter are not expected to mix. However, since the operators ϕ v andφ −v connect to auxiliary fields propagating in opposite directions, the mixing pattern allowed by chiral symmetry breaking is different from the straight-line operators used for quasi-PDFs: . This was pointed out by one of us in Ref. [70] and later confirmed at one-loop order in lattice perturbation theory [71]. The latter calculation also found that more generally, the pattern of mixing depends only on the direction with which the Wilson lines hit the quark fields, which is consistent with the prediction of the auxiliary field approach. We note that a numerical study of nonperturbative renormalization of staple-shaped operators was recently done in Ref. [72].
A RI-xMOM renormalization condition for cusp operators is straightforward to formulate. This requires the position-space Green's function for the cusp operator, G C ðξ 0 ; ξÞ ≡ hζ n 0 ðξ 0 n 0 ÞC n 0 ;n ð0Þζð−ξnÞi QCDþζ Performing a position-space amputation, a possible condition is

B. Off-axis paths
It is possible to somewhat relax the constraint that n points along a lattice axis. For example, let n ¼ 1 ffiffi 2 p ðx þŷÞ.
On the lattice, one definition of the straight-line operator is [65] O lat Γ ð0; ξ; nÞ ≡ψðξnÞΓ where W xy is formed from the zigzag product of gauge links alternating between the x and y directions and vice versa for W yx (see Fig. 4). This average is necessary so that the operator has a simple transformation under n-parity. These zigzag Wilson lines can be obtained as propagators of auxiliary fields ζ xy and ζ yx that are defined using twolink covariant derivatives along the n direction, e.g., Because of rotational symmetry breaking on the lattice, the mass counterterm m xy ¼ m yx will in general be different from the on-axis case. 5 The quark bilinear takes the form O lat Γ ð0; ξ; nÞ ¼ where ϕ xy ¼ζ xy ψ, etc. Defining φ AE ≡ 1 ffiffi 2 p ðϕ xy AE ϕ yx Þ, this can be rewritten as where the additional cross terms involving e.g., hζ xyζyx i vanish. Under n-parity, φ AE → AE= nφ AE .
FIG. 3. Staple-shaped operator. 5 In the continuum, the case of an auxiliary field propagating along a multicusp curve that approximates a smooth one was considered in Ref. [30]. In that case, when the number of cusps goes to infinity the action is equal to that of a field propagating along the smooth curve with an added mass term that accounts for the effect of the cusps. This is analogous to the situation here, where a straight line is approximated by zigzags that become infinitely many as the lattice spacing goes to zero.
We find that additional mixing can occur. Defining v ¼ 1 ffiffi 2 p ðx −ŷÞ, φ þ mixes with = nφ þ , = vφ − , and = n= vφ − , the last of which is a chiral-even mixing. This leads to correlation functions of the form hφ þ φ − i, which contain a difference of Wilson lines W xy − W yx . One expects that this difference will vanish in the continuum limit. This expectation can be justified in the auxiliary field approach: in the continuum, there is an SUð2Þ flavor symmetry relating ζ xy and ζ yx . It is broken in the next order of the Symanzik expansion, where a term of the form aζ xy G μν n μ v ν ζ xy − aζ yx G μν n μ v ν ζ yx can appear in the Lagrangian; this implies that hφ þ φ − i is OðaÞ. Therefore, even though the mixing between φ þ and operators containing φ − is equal dimensional, the most serious new mixing effect in correlation functions is suppressed by at least OðaÞ. The remaining effect is that even when using a chiral fermion action, renormalization will be different for operators where = v= nΓ= n= v is equal to þΓ and those where it equals −Γ.
An alternative definition is to average the paths locally. Taking the average of the two local link paths, the lattice action for ζ avg is defined using the covariant derivative, Again, the mass counterterm will in general be different from the previous case. Then, using the bilinear ϕ avg ¼ζ avg ψ, the bilocal operator hφ avg ðξnÞΓϕ avg ð0Þi ζ effectively averages over 2 N link paths, where ξ ¼ ffiffi ffi 2 p Na. In this case, the pattern of equal-dimensional mixing is the same as the on-axis case.

V. CONCLUSIONS
The auxiliary field approach is an invaluable tool for understanding the renormalization and improvement of Wilson-line operators on the lattice. Generically, these nonlocal operators can be represented using local operators in an extended theory involving auxiliary fields. We have shown how this can be done for a variety of operators, beyond the simplest ones involving straight on-axis link paths used for quasi-PDF studies.
Using the Symanzik expansion, we were able to study discretization effects and the form of the improved operators. As foreseen in Ref. [28], we find that the leading effects are linear in a even if chiral symmetry is preserved on the lattice. Likewise, we find that working at maximal twist does not automatically eliminate the effects linear in a, but it can remove some of the contributions and it can produce some simplification by reducing the number of improvement coefficients. Our analysis provides a general framework for the OðaÞ improvement of nonlocal operators used in the quasi-PDF approach for computing PDFs on the lattice.
In order to apply this improvement to a lattice calculation of quasi-PDFs, the relevant coefficients for the choice of lattice action must be determined. For some actions, many of these are already known from lattice perturbation theory calculations done for the static quark theory. There are plans to determine the coefficients for actions used by ETMC [73], and we intend to study the effect of improvement on the approach to the continuum limit. It will be important to establish that discretization effects in quasi-PDF calculations can be controlled. The static quark theory [35,36] is defined using a spinor Q that satisfies Q ¼ P þ Q, where P AE ¼ ð1 AE γ 0 Þ=2. Its 2 spin degrees of freedom do not couple in its action, and its propagator is the same as the ζ-field propagator with n ¼t, multiplied by P þ . This means that, after accounting for the projector, the static quark and auxiliary field theories can be identified with each other. It will be convenient to define the projected bilinears ϕ AE ¼ 1 2 ð1 AE = nÞϕ, where the projectors become P AE when n ¼t. In our approach, neglecting a twisted mass term, the renormalized OðaÞ-improved operators take the form where Z AE ϕ ¼ Z ϕ ð1 AE r mix Þ, b AE ϕ ¼ ðb ϕ AEb ϕ Þ=ð1 AE r mix Þ, and c AE ϕ ¼ ðc ϕ AEc ϕ Þ=ð1 AE r mix Þ. In Ref. [50], the time components of the static-light axial and vector currents are defined as By identifying which projector is contracted withψ, this lets us make the identifications and the corresponding renormalization factors can be equated: For the static-light axial current, Ref. [50] identifies four possible OðaÞ improvement operators: At leading order, the equations of motion for Q give ðδA stat 0 Þ 2 ¼ 0 and those for ψ give ðδA stat 0 Þ 1 − ðδA stat 0 Þ 3 þ ðδA stat 0 Þ 4 ¼ 0; these are used to eliminate ðδA stat 0 Þ 2;3 . The improved, renormalized operator is then given by We can make the identification To relate Eq. (A1) with Eq. (A9), we need to use the equations of motion to eliminate ðδA stat 0 Þ 1 and insert ðδA stat 0 Þ 2 . Our identification allows us to equate the improvement coefficients at leading order: Likewise, Ref. [50] gives the improved, renormalized vector current as and by again using the equations of motion we obtain at leading order b þ ϕ ¼ c stat V þ b stat V and c þ ϕ ¼ c stat V . For a chiral action, our identification leads to Z stat A ¼ Z stat V , c stat A ¼ c stat V , and b stat A ¼ −b stat V , consistent with Ref. [60].

APPENDIX B: COMPARISON WITH WHOLE-OPERATOR APPROACH
In Ref. [49], the symmetry properties of dimensionthree operators of type O Γ ðξÞ and similar dimension-four nonlocal operators were studied in order to understand mixing and OðaÞ effects. We find the same pattern of mixing with three-dimensional operators and four-dimensional operators proportional to m q . However, in general it is much less constraining to consider the operator as a whole rather than using the auxiliary field approach to represent it using two local operators. This leads to the following differences: (1) The auxiliary field approach implies that for a chiral fermion action, the renormalization of O Γ ðξÞ is independent of Γ and depends only on two parameters Z ϕ and m. When chiral symmetry is broken, the splitting between different Γ and the mixing are controlled by a single parameter, r mix . In contrast, the whole-operator approach implies a generic ξdependent and Γ-dependent renormalization. 6 (2) In the auxiliary field approach, the four-dimensional operators with derivatives can only have those derivatives inserted in either local operator, i.e., effectively at either end of the Wilson line. 7 By using the equations of motion for the fermion field and the auxiliary field, the number of improvement terms with derivatives can be significantly reduced. On the other hand, Ref. [49] found operators with = D ⊥ or = nðn · DÞ inserted at any point ξ 0 ∈ ½0; ξ along the Wilson line. This large number of operators makes it appear impractical to attempt OðaÞ improvement. (3) The local operator ϕ ¼ζψ transforms under the fundamental irrep of the fermion flavor symmetry group. This means that there is no mixing of O Γ with nonlocal gluonic operators of the type discussed in Ref. [34]. On the other hand, the whole-operator approach would predict that those gluonic operators could have a divergent Oða −1 Þ contribution from a flavor singlet nonlocal quark operator O Γ . 6 Note that whole-operator nonperturbative renormalization done using RI-MOM type schemes [6,9,15,72] has generically found a dependence on Γ. However, this is by construction in the definition of observables used for imposing renormalization conditions. Once converted to a minimal scheme such as MS, the pattern predicted by the auxiliary field approach is recovered [37], up to the precision of the scheme conversion. 7 Insertions in the Wilson line can occur as a result of higherdimensional corrections to the action of the auxiliary field.