A QCD analysis of near-threshold quarkonium leptoproduction at large photon virtualities

We propose a novel approach to compute the cross section of near-threshold $J/\psi$ and $\Upsilon$ production in electron-proton scattering at large photon virtualities $Q^2$ based on an operator product expansion. We show that the process can be used to extract the gluon part of the D-term gravitational form factor of the proton. At the subleading level, it is also sensitive to the trace anomaly effect of QCD.

numerically evaluate the cross section and study the impact of the D-term as well as the gluon condensate. Section VI is devoted to conclusions.

II. KINEMATICS
We shall be interested in the near-threshold production of a heavy quarkonium vector meson H with mass M in electron-proton scattering ep → e γ * p → e Hp . We have in mind H = J/ψ and Υ. The center-of-mass energy of the virtual photon-proton subsystem at the threshold is where m N = 0.94 GeV is the proton mass. Numerically, W th ≈ 4.04 GeV for J/ψ and W th ≈ 10.4 GeV for Υ. q and p are the virtual photon and proton momenta, respectively, with Q 2 = −q 2 being the photon virtuality. Near the threshold, the Bjorken variable takes the form Eq. (2) shows that, unlike the usual situation in DIS, Q 2 and x B are not independent variables. x B approaches unity as Q 2 goes to infinity. We also see that, somewhat counterintuitively, there is no kinematical restriction in Q 2 . Threshold production can occur even when Q 2 is arbitrarily large. Using the standard variables in DIS, S ep = (p+ ) 2 , y = p · q/p · where µ is the incoming electron momentum, we can write Small-W does not necessarily imply small-S ep when Q 2 is large. In particular, the process can be studied at the future EIC [30]. Let p cm and k cm be the 3-momentum of the incoming proton and outgoing quarkonium, respectively, in the center of mass frame of the γ * p subsystem.
The momentum transfer is At the threshold, k cm = 0 so that We see that |t th | is minimal in photoproduction Q 2 = 0 and monotonously increases with increasing Q 2 . In the heavy-quark mass limit, Q 2 |t th |. Away from the threshold, t takes a value in the range |t max | > |t| > |t min | depending on the angle between p cm and k cm . The differential cross section is given by where J µ em = f e fqf γ µ q f is the electromagnetic current (e f being the charge in units of |e|) and T /L refers to the transversely (T ) or longitudinally (L) polarized virtual photon. The factor 1/2 is for averaging over proton helicities. The nontrivial dynamics of QCD is contained in the hadronic matrix element d 4 ye −iq·y p k|J ν em (y)|p = (2π) 4 δ(k + p − q − p) p k|J ν em (0)|p .
Computing (8) from first principles in QCD is a challenging task. Most of the previous theoretical works have focused on the photoproduction limit. In contrast, in this paper we shall investigate leptoproduction in the large Q 2 M 2 region.

III. OPE AT LARGE Q 2
In this section, we formulate our strategy to calculate the hadronic matrix element p k|J em |p near threshold at large Q 2 M 2 . We have chosen to work in the high Q 2 region for reasons to become clear shortly. For definiteness, we consider J/ψ production, but the case with Υ is completely analogous. In fact, our approach is better justified when M m N . Thus, Υ production is more preferred from a theoretical point of view, though of course experimentally it is more challenging.
Let us first mention that, if the center of mass energy is sufficiently high s = W 2 M 2 , |t|, the process is commonly called Deeply Virtual Meson Production (DVMP). The cross section is known to factorize in perturbative QCD in terms of the generalized parton distribution (GPD) and the meson distribution amplitude (DA) [31]. Near the threshold, s = O(M 2 ), and |t| is comparable to, or even exceeds s depending on the value of Q 2 , see (6). Note, however, that s is small because of the cancellation s = 2p · q − Q 2 + · · · and x B = Q 2 /2p · q stays close to unity. Moreover, 2p · q ∼ Q 2 |t| at least parametrically when M m N (see (6)). This gives us some hope that a perturbative approach is possible.
Our basic argument is that near the threshold, the amplitude (8) is related to the following current-current correlator * where µ (k) is the J/ψ polarization vector. This matrix element is similar to the one that appears in doubly virtual Compton scattering (DDVCS) γ * (q)p → γ * (k)p , or timelike Compton scattering (TCS) in the special case q 2 = 0 (see, e.g., [32]). However, there is a crucial difference. The DDVCS amplitude is given by the correlator J em J em , and is dominated by the light quark degrees of freedom (light quark GPDs) except in the very small-x B region where it is dominated by gluons. In (9), on the other hand, one of the electromagnetic currents has been replaced by the charm quark current operator (bottom quark, in the case of Υ production). As a result, only the charm component of the other J em is relevant, and the matrix element becomes primarily sensitive to the gluonic content of the proton. That the J/ψ production amplitude is related to a DDVCS-like (photon production) amplitude is intuitively reasonable, in view of the fact that in actual experiments, a J/ψ and a timelike photon with virtuality exactly at the J/ψ mass are practically indistinguishable as they are probed via leptonic final states (e + e − pairs). However, in DDVCS or TCS, the resonance region k 2 ≈ M 2 is usually avoided because the nonperturbative final state effect to produce the vector meson comes into play. As a function of k 2 , the right hand side of (9) has a sharp resonance peak near the J/ψ mass shell. Using the LSZ reduction formula, we can write where e f = 2/3 for the charm quark and Γ is the total width of J/ψ. The decay constant is related to the electromagnetic width as Away from the very narrow peak (note that Γ = 93 keV M = 3.1 GeV), the current correlator is expected to behave smoothly. We thus arrive at the relation where β is a c-number of order unity which is not under control. It is understood that the right hand side is evaluated close to, but not too close to the J/ψ mass shell |k 2 − M 2 | M Γ. Our key observation is that in this off-mass-shell region, one can perform an operator product expansion (OPE) when Q 2 is large.
Before doing so, a few additional remarks are in order. (i) On general grounds, one expects corrections to (10) from higher resonances which the operatorcγ µ c can excite. However, this effect will be suppressed near threshold because, at fixed values of s and t, only the resonances with mass smaller than √ s − m N can be produced. There may also be contributions from the deep Euclidean region k 2 < 0 if one considers a dispersion relation for the current correlator in k 2 similarly to [33,34]. 1 Such an analysis may lead to a more precise evaluation of the quarkonium matrix element and help to determine the value of β. We leave this to future work. (ii) Our argument is similar in spirit to the vector meson dominance (VMD) hypothesis. Note that this is different from the VMD assumption used in many literature works on J/ψ photoproduction mentioned in the introduction [3,4,21,22]. In these works, VMD has been applied to the incoming massless photon γ → J/ψ. In photoproduction, this results in a significant mismatch between the initial and final virtualities 0 → M 2 . Here, in a sense, we apply VMD in a reverse way to the outgoing J/ψ → γ * (cf. [8,33,35]). While the difference in virtualities partly remains, this has little impact on the overall kinematics of the reaction because |k 2 − M 2 | Q 2 , |t|. (iii) On the other hand, our approach is different from the nonrelativistic (NR)QCD framework [36] which is commonly used for quarkonium production in hadronic collisions. In NRQCD, the charm and anticharm quarks in the perturbative amplitude which couple to an external J/ψ are both on-shell to leading order in the velocity expansion. However, in (12), c andc are far off-shell with virtuality of the order of Q 2 (see below). Moreover, (9) assumes that J/ψ is produced only in a color-singlet state. This is reasonable because near the threshold, all the energy has to be used to create a J/ψ, and there is little phase space for extra gluon emissions.
Let us now discuss the OPE. The current correlator on the right hand side of (12) can be written as where r = x − y. The product of currents can be expanded if the relative distance |r µ | is small, which is the case when the momentum k+q 2 is deeply spacelike. From (6), near the threshold, Therefore, This can be made arbitrarily negative by choosing Q 2 M 2 . We need to also make sure that the large momentum Q does not 'leak' into the proton vertex which in practice means Q 2 |t|. Very close to the threshold, this is satisfied if M m N . As one goes away (but not too far away) from the threshold, the condition Q 2 |t| is well satisfied when t ∼ t min . 2 As we shall see in Section V, t ∼ t min is the most interesting region.
However, for technical reasons the 'symmetric' form (13) is not very convenient. Being a nonforward matrix element, it can be expressed in several different 'frames' The meaning of the OPE is different in different frames. The final result must be the same, but this equivalence is often difficult to see. We shall return to this issue later. For the moment we find it most convenient to start with the middle expression of (16). We evaluate it as where S is the charm quark propagator in the presence of background gluon fields. Since we work in the regime Q 2 M 2 , |r µ | is typically much smaller than 1/M and the heavy quark mass m c ≈ M/2 can be neglected to first approximation. An important point of our approach is that we shall expand (17) in terms of local operators [39,40], instead of nonlocal light-cone correlators as is usually done in high energy scattering. Near the threshold, the role of light-cone directions appears to be less conspicuous. More importantly, the OPE with local operators is well suited for our purpose of establishing a connection to the D-term which is the matrix element of the (local) energy momentum tensor operator.
Consider the first line on the right hand side of (17). The lowest contribution comes from the operatorcγ µ γ 5 c, followed by higher dimensional operators such ascγ µ D ν c,cγ µ F αβ c andcF αβ γ β γ 5 c. The (nonforward) matrix elements of these operators measure the intrinsic charm contents of the proton which are in general believed to be tiny (see however, [5]). In this paper we simply neglect all of them, although they can be straightforwardly restored if need arises.
We thus focus on the second line of (17). Basically, we only keep dimension-4 purely gluonic operators. This in particular includes the gluon part of the QCD energy momentum tensor whose proton matrix element p |T αβ g |p is what we are ultimately interested in. However, certain higher dimension operators are a priori not suppressed. As in usual DIS or DVCS, the contribution of the leading twist operator with Lorentz spin-j is proportional to (2q · p/Q 2 ) j ∼ (1/x B ) j and x B ∼ 1 for our problem. The difficulty to sum over these higher spin operators with j > 2 is the reason why the local version of the OPE is not commonly used in DVCS. Here, however, we do not attempt to perform this summation. Among the twist-two operators, the energy momentum tensor T g αβ with j = 2 dominates in the sum when Q 2 is sufficiently large. The contributions from the other twist-2 operators with spin j > 2 are relatively suppressed because their anomalous dimensions are nonvanishing. Admittedly, the rate of this suppression is slow, only logarithmic in Q 2 , so a large leverage in Q 2 is needed to isolate the spin-2 component. While this may seem a difficult task, we point out that a very similar problem exists in the current strategy to extract the quark D-term from the DVCS data [27][28][29]. The subtraction constant in the dispersion relation between the real and imaginary parts of the Compton form factor, commonly denoted by ∆(t) [26], is given by the sum of infinitely many Gegenbauer coefficients ∆(t, In order to isolate the quark D-term ∝ d 1 (t) which has the same anomalous dimension as the energy momentum tensor, one needs a large leverage in Q 2 to disentangle different moments. Assuming that such an analysis is feasible at the future EIC, we expect that the same can be done for the gluon D-term.
We shall work in Fock-Schwinger gauge r µ A µ (r) = 0 for actual calculations. In this gauge, in the small-r limit, the massless quark propagator in the background gluon field is given by, in d = 4 − 2ε dimensions (see for example, [37,38]) where F αβ = F αβ a t a with Tr(t a t b ) = δ ab /2 and our convention for the covariant derivative is D µ = ∂ µ + igA µ . In the denominators, r 2 is short for r 2 − i . In (19), we have kept only the terms which contribute to dimension-4 gluonic operators F F . At first sight, the dimension-3 operators of the form D α F βγ are irrelevant because they are matrices in color space so when inserted in (17), they either vanish after tracing over color indices or lead to operators with dimension-5 or larger. However, for the present problem, it turns out that they cannot be neglected. We shall discuss this later. Note that, since the Fock-Schwinger gauge breaks translational invariance, in general S(r, 0) = S(0, −r). However, for the terms listed in (19), the relation S(r, 0) = S(0, −r) actually holds.
In the second line of (17), the unit operator can be neglected because we are computing the nonforward amplitude p |1|p = 0. Consider then the O(g 2 F F ) terms in (19) which lead to a logarithmically enhanced contribution as implied by the prefactor Γ(−ε). Taking the trace of the g 2 F F terms in (19) in color space, we find whereT αβ g is the traceless part of the gluon part of the QCD energy momentum tensor. (20) explicitly shows that the logarithmic part is insensitive to the trace anomaly. Inserting the first term of (20) into the second line of (17), we find Note that in the last step we have dropped the divergent piece 1/ε. It can be absorbed into the renormalization of the twist-two operatorT µν c ∼cγ (µ D ν) c contained in the first line of (17). The coefficient αs 3π can be identified with the anomalous dimension γ c←g of this operator. As we already mentioned, we neglect the matrix element of (renormalized)T µν c so in practice the 1/ε simply disappear. The non-logarithmic terms in (22) combine with those from the second term of (20) and the square of the O(gF ) term in (19). After a tedious but straightforward calculation, we arrive at the total α s F F contribution where the operators are defined at the scale µ R . This is manifestly transverse with respect to q, i.e. q µ A µν = A µν q ν = 0, as a consequence of the Ward-Takahashi identity. In Appendix A, we show that the forward matrix element of (23) reproduces the 1-loop coefficient functions of the DIS structure functions. However, (23) has an obvious problem. The tensor A µν is transverse with respect to q µ and q ν , but this is because we have started with the middle expression in (13). In the present problem, gauge invariance rather implies k µ A µν = A µν q ν = 0. Actually, problems of this kind typically arise in off-forward kinematics. It is known that ensuring the electromagnetic gauge invariance of DVCS amplitudes is a highly nontrivial issue [48,49]. The leading order (leading twist) result does not fully satisfy the WT identity, and one has to include higher twist corrections to restore it. In the context of OPE, this amounts to including operators with total derivatives [50]. In Appendix B, we demonstrate that the dimension-3 operators D α F βγ which were neglected in (17) indeed give rise to total derivative operators. This calculation suggests that a complete treatment of the problem requires the inclusion of dimension-5 and even dimension-6 operators in the expansion (17), which is beyond the scope of this work. Here instead, we suggest an ad hoc solution of the problem. In (23), we set q 2 = −µ 2 R to eliminate the logarithmic terms. In the remainder terms, we implement the following minimal modifications 3 to make A µν transverse with respect to k µ and q ν , and symmetric in q and k where the coupling and the operators are evaluated at the scale O(Q 2 ). In the 'leading-twist' approximation, one can further simplify (see (A2)) −F µα F νβ ≈ 1 2 (g µνT αβ g − g µβT αν g − g ανT µβ g + g αβT µν g ).
Actually, since we are neglecting the twist-2 operators with spin j > 2, it is not entirely consistent to include anything beyond (25) as it corresponds to twist-4 effects. Still, for phenomenological purpose it may be interesting to include at least the trace part of T αβ g in order to assess the impact of the trace anomaly. 3 There is an ambiguity when replacing q 2 with q · k = q 2 − q · ∆, since q 2 → q+k 2 2 = q 2 − q · ∆ + ∆ 2 /4 seems to be an equally good choice (cf. (13)) . However, the difference is subleading because q 2 ≈ 2q · k ∆ 2 in the present kinematics, see (14). This ambiguity can only be resolved by including the dimension-6 operator ∂ 2 T αβ .

IV. TWO-GLUON FORM FACTORS
In order to compute the actual cross section, we need to parametrize the non-forward matrix element of two-gluon operators in (24) in terms of form factors. First we have the gravitational form factors at our disposal [42] p |T µν g |p =ū(p ) A g γ (µ P ν) + B g where ∆ µ = p µ − p µ , P µ ≡ p µ +p µ 2 and A (µ B ν) ≡ (A µ B ν + A ν B µ )/2. All four form factors are functions of t = ∆ 2 and the renormalization scale µ R in the MS scheme. D g is the gluon part of the D-term form factor which we are mainly interested in. (In the literature often the notation C g = D g /4 is often used.) TheC g form factor is related to the trace anomaly [24]. The traceless part reads Next consider the two gluon operator with four open indices Its most general parametrization consistent with parity, hermiticity and time-reversal symmetry is 4 The seven form factors can be partly constrained by requiring consistency with (26). Contracting the indices αβ in (30), we get 4 Terms which contain the antisymmetric tensor µαρλ are not independent. For example, the following identity holds = A gū (p )γ (µ P ν) u(p) where in the second equality we used the following relations which can easily be obtained by term-by-term comparison: We see that only two linear combinations of X, Y, Z enter these relations. By comparing the coefficients of g µν , one should be able to obtain another relation between W andC g . However, this is nontrivial due to the presence of the QCD trace anomaly. If one naively contracts the indices µν in (31) and computes the matrix element of T µν g by forming the linear combination (18), one ends up with a wrong relationC g = −A g /4 (in the forward limit) and W is undetermined. The problem is intimately tied to operator renormalization. In dimensional regularization, the following innocent-looking relation does not hold Namely, operator renormalization and trace operation do not commute. The correct way to proceed is to write on the left hand side of (31) and sum over the indices µν using (26) and (31). This gives On the other hand, the matrix element p |F 2 |p has to be carefully evaluated in a chosen regularization scheme [24,25] (see also [43]). In dimensional regularization, it is given by a linear combination of the gravitational form factors, see Eq. (13) of Ref. [12] p where the quark gravitational form factors A q , B q , D q ,C q are defined analogously to (26) for the quark part of the energy momentum tensor. The coefficients K q,g are defined in [12] and can be evaluated, in principle, to arbitrary order in perturbation theory. At the moment, the three-loop results are available [24,25]. They depend on the number of flavors and the renormalization scale via the QCD coupling α s (µ R ). (37) and (38) give a complicated relation between X, Y, Z, W and the quark and gluon gravitational form factors. In the forward limit t = 0 it somewhat simplifies and we find where we used A q (0) + A g (0) = 1 andC q (t) +C g (t) = 0. The relation to the parameter b often used in the literature [44] is where γ m is the mass anomalous dimension. b is the partition of the trace anomaly into the quark and gluon condensates. It is scheme and scale dependent.

V. NUMERICAL RESULTS
In this section we show numerical results for the differential cross section based on the formula (24). We do not intend to perform a complete calculation which is anyway not possible at the moment as it requires the detailed knowledge of all seven form factors A(t), B(t), ... On the other hand, some information about the gravitational form factors A g , B g , D g ,C g is already available in the literature. Based on this, we consider two interesting cases which allow us to make a quantitative prediction. Case 1: We use the 'leading-twist' approximation (25) and keep only the traceless part of the energy momentum tensor (27). As explained in Section III, in doing so we assume that the contributions from the twist-2 operators with spin j > 2 can be either neglected or separated out by using a large leverage in Q 2 . Case 2: We evaluate the full two-gluon operators (24) including the trace part of the energy momentum tensor. While this is not a consistent approximation (because we keep the twist-4 effect and neglect the twist-2, spin-j > 2 contributions), it is an instructive exercise to assess the impact of the trace anomaly. In both cases, we set B g = 0 following the suggestion from lattice QCD (see, e.g., [45]) that this form factor is numerically small. In Case 2, we also set X = Y = Z = 0 as we know nothing about these form factors. On the other hand, the W form factor is related to the trace anomaly and will be given full consideration.
We use the following parametrization of the gravitational form factors The tripole form of the D-term is motivated by the quark counting rule [46]. Since the form factors are evaluated at a large scale µ 2 R = Q 2 , to first approximation we can use the asymptotic results 2Nc and n f = 3 represents the number of light flavors in the proton. The value D g (0) is our main interest and should be determined by future experiments. Here, for the sake of demonstration, we use the results of a recent lattice simulation D g (0) = −7.2 (C g (0) = −1.8) with m A = 1.13 GeV and m C = 0.76 GeV at µ R = 2 GeV [47]. (We neglect the scale dependence of these parameters.) On the other handC g at zero momentum transfer is related to the QCD trace anomaly [24]. Asymptotically µ 2 R → ∞, where b is introduced in (40). To one-loop, β 0 = 11N c /3 − 2n f /3 and γ m = 3C F αs 2π . A more precise expression can be found in [24,25].
Under these assumptions, (37) and (38) reduce to a simple formula For simplicity, we use the one-loop result K q,g where (α s = α s (Q 2 ))  See [25] for the three-loop result. Finally, the square of the prefactor in (9) is evaluated as (including the factor e f from (13)) where we set e c = 2/3, M = 3.1 GeV, Γ e + e − = 5.55 keV and α em = e 2 /4π = 1/137. The corresponding value for Υ is where e b = −1/3, M = 9.46 GeV and Γ e + e − = 1.34 keV. The parameter β should be determined by fitting the data (for example the total cross section at some value of W ) for each quarkonium species. In the numerical results below we set |β| = 1. In Fig. 1, we show the total and differential cross sections for J/ψ at Q = 8 GeV, α s (Q) = 0.2. The latter is evaluated at W = 4.4 GeV. In both plots, the upper and lower dashed curves correspond to Case 1 with D g = 0 and D g = −7.2, respectively. We see a dramatic impact of the gluon D-term. 5 A negative (positive) D-term tends to shrink (enhance) the differential cross section. The same tendency has been observed in [9] in the case of photoproduction Q 2 = 0. The upper and lower solid curves correspond to Case 2 with b = 1 (zero gluon condensate) and b = 0 (zero quark condensate), respectively. We see that the dependence on the parameter b is significant. We also see that the gluon condensate tends to reduce the cross section, which is actually opposite to what was found in [9]. It is not clear to us whether this is due to the fact that different processes were considered (photoproduction vs. leptoproduction), or perhaps due to the deficiency of the model used in [9].
Next, in Fig. 2 we show the result for Υ at Q = 18 GeV, α s (Q) = 0.16. Near the threshold (W th = 10.4 GeV), the cross section becomes very small. In the right panel we selected a somewhat large value W = 12.5 GeV considering the realistic luminosity of EIC. Again we see a large effect of the D-term. However, the impact of the trace anomaly and the split between b = 1 and b = 0 are barely visible.
In photoproduction, the J/ψ total cross section is about 1 nb at W = 4.5 GeV [15]. In leptoproduction, we see that the cross section is several orders of magnitude smaller. For Υ production, there is another 2 orders of magnitude suppression. Thus, near-threshold leptoproduction is a luminosity-hungry observable. Moreover, as explained in Section III, one needs a large leverage in Q 2 to extract the D-term. Given these requirements, we think that the best place to test our proposal is J/ψ (and possibly also Υ) production in the high luminosity mode of EIC.  FIG. 2: Υ total and differential cross sections at Q 2 = 18 2 GeV 2 . See the caption of Fig. 1 for the explanation of each curve.

VI. CONCLUSION
In this paper we have proposed a novel strategy to compute the cross section of near-threshold quarkonium production at large momentum transfer. Compared to photoproduction, near-threshold leptoproduction has so far attracted much less attention due to the lack of strong phenomenological motivations. We have demonstrated that the process is useful for probing the gluon D-term, quite complementary to the ongoing effort to extract the quark D-term in DVCS. The possible impact of the D-term on the differential cross section dσ/dt has been already pointed out in the case of photoproduction using holography [9,13]. In leptoproduction at large Q 2 M 2 , the problem can be studied within the perturbative framework. Moreover, at the subleading level the cross section is also sensitive to the value of the parameter b defined in Eq. (40) which characterizes the structure of the QCD trace anomaly. The proposed measurements require high luminosity and a large leverage in Q 2 . The only machine that can deliver these requirements is the future EIC.
Our analysis in this paper is only the first step and there are a number of directions for future research. In particular, it is interesting to see if a similar approach can be applied to photoproduction using the heavy quark mass as a hard scale. On the phenomenological side, the contribution from the Bethe-Heitler process e − → e − γ * needs to be investigated as J/ψ and Υ are reconstructed from lepton pairs in actual experiments. The seven form factors introduced in (30) can be calculated in lattice QCD. The values of t considered in this paper are rather large, and the extrapolation to the forward limit is a serious challenge. Lattice calculations of these form factors will be very valuable in this respect.
where A g is the fraction of the proton momentum carried by gluons. We then decompose the operator F µα F νβ as 6 and extract the energy momentum tensor component T αβ g ∼ −F αλ F β λ . The remainder tensor C µανβ has the same symmetry as F µα F νβ except that it is traceless with respect to any pair of indices. Its forward matrix element vanishes. We thus have p| − F µα F νβ |p ≈ A g (g µν p α p β − g µβ p α p ν − g αν p µ p β + g αβ p µ p ν ). (A3) This gives i d 4 re ir·q p|cγ µ c(0)cγ ν c(−r)|p q α q β (g µν p α p β − g µβ p α p ν − g να p µ p β + g αβ p µ p ν ) .
From this one can read off the known one-loop coefficient functions for the DIS structure functions where we used the identity ∂ ν T µν g = F µ ν D α F αν . Note that this is antisymmetric in µ and ν. The 1/ε divergence in (B3) is absorbed into the renormalization of the operator cγ µ gF α λ γ λ γ 5 γ ν c − (µ ↔ ν) ∼ µλνρc γ ρ gF α λ c, which comes from the first line of (17). The remaining terms contain total derivative operators. In the nonforward matrix element, the derivative operator is replaced by the momentum transfer p |∂ ρ T ρν g |p = i∆ ρ p |T ρν g |p where ∆ ρ = q ρ − k ρ . This is how, in principle, total derivative operators from higher dimensional terms can restore the WT identity through the addition of ∆ corrections. However, (B3) is not sufficient to make the logarithmic terms in (23) transverse with respect to k µ . For that, we would need operators like ∂ ν T µβ g and g µν ∂ α T αβ g . We presume that the missing terms come from the dimension-5 and dimension-6 operators in the expansion of S(r, 0). We leave this to a future work.