Gauge invariance and kaon production in deep inelastic scattering at low scales

This paper focuses on hadron mass effects in calculations of semi-inclusive kaon production in lepton-Deuteron deeply inelastic scattering at HERMES and COMPASS kinematics. In the collinear factorization framework, the corresponding cross section is shown to factorize, at leading order and leading twist, into products of parton distributions and fragmentation functions evaluated in terms of kaon- and nucleon-mass-dependent scaling variables, and to respect gauge invariance. It is found that hadron mass corrections for integrated kaon multiplicities sizeably reduce the apparent large discrepancy between measurements of $K^+ + K^-$ multiplicities performed by the two collaborations, and fully reconcile their $K^+/K^-$ ratios.


II. LEADING ORDER MULTIPLICITIES AT FINITE Q 2
The z-integrated hadron multiplicities measured by the HERMES collaboration [14,19] are defined as a ratio of the semi-inclusive to inclusive cross sections, where x B = Q 2 2p·q is the Bjorken scaling variable, Q 2 = −q 2 the virtuality of the exchanged photon, z h = p·p h p·q is the fragmentation invariant, and the remaining kinematic variables are defined in Fig. 1 left. The COMPASS collaboration has measured integrated multiplicities as averages over y of the differential ones, dz h M h (x b , y, z h ) y [20,21]; however, since y < 0.7 within the COMPASS kinematic cuts, the two definitions are approximately equivalent, and in this paper we will use Eq. (1) for both experiments. The integration over the initial state DIS invariants is performed over the experimental kinematic acceptance of each measurement [20,26], denoted in short by "exp". In more detail, the integral over dx B is performed over the bin of nominal value x exp B , and the integration over dQ 2 is performed within x B -dependent limits defined by the kinematic cuts of each experiment. These cuts, as well as plots of the (x B , Q 2 ) phase space with the experimental x B bins are shown in Fig. 2. The integration limits on z h are explicitly indicated in Eq. (1), with the COMPASS value in parentheses when different from that used at HERMES. As noted in Ref. [26], it is important to perform the full integration in Eq. (1), rather than evaluating the cross section at an average value for, in particular, Q 2 .
In order to study Hadron Mass Corrections in SIDIS, we will consider Nachtmann-type scaling variables defined by light-cone fractions of the photon's momentum q and, respectively, the proton and hadron momentum. In the so-called "(p, q) frame" [22], in which p and q are collinear in 3-dimensional space and oriented along the z−direction, i.e., have zero transverse component (p T = q T = 0), one finds where M is the nucleon target mass and m h is the mass of the detected hadron [23]. In the case of ζ h , note the interplay between initial and final state masses and Lorentz invariants that complicates the analysis. Had the data been binned in z e = −p · p h /q 2 , there would have been no mixing [23]. One can easily verify that in the Bjorken limit, where M 2 /Q 2 → 0 and m 2 h /Q 2 → 0, one recovers the usual massless scaling variables x B and z h .

A. Collinear factorization with non-zero average parton virtualities
At leading order (LO) in the strong coupling constant, one needs to calculate the diagram in Fig. 1 left. The resulting hadronic tensor reads where Φ q and ∆ h q are quark-quark correlators associated with the quark distribution and fragmentation functions, respectively [27][28][29][30]. For clarity of presentation, we are here considering only one quark flavor. In general, the right hand side would sport a charge-weighted sum over quark and antiquark flavors; details can be found in Ref. [23].
Obtaining a factorized expression for the hadronic tensor (4) requires two steps. The first one is an expansion of the quark-quark correlators in inverse powers of the leading components of the parton momenta entering and exiting the hard scattering H in Fig. 1, namely, k + and k − . In this paper, we limit our attention to the first order terms in this expansion, and write Φ q = k + φ 2 (k)n / + O(1/k + ) and ∆ q = k − δ 2 (k )n / + O(1/k − ) . The lowercase Greek letters indicate scalar functions of the momenta, and the unit light-cone plus-vector and minus-vector are denoted, respectively, by n µ andn µ . One then obtains In this expression, there is a clear separation between the partonic "twist-2" correlation functions φ 2 and δ 2 on the one hand, and on the other hand the dynamics and overall kinematics of the hard scattering process (namely the trace and δ-function). Note that the neglected pieces are not forgotten, but in fact contribute to restore gauge invariance in "higher-twist" (HT) diagrams that include additional parton exchanges between the hard scattering and the non-perturbative blobs in Fig. 1 [28,31,32]. The second step consists in approximating the incoming and outgoing partonic momenta appearing in the fourmomentum conservation δ-function, namely k ≈ k and k ≈ k . It is important to remark that this is the only place where we perform an approximation rather than a controlled expansion. Contrastingly, the trace part is expanded in inverse powers of the plus and minus light-cone momenta; one could then improve on this approximation by retaining higher order terms, and considering in addition diagrams involving multi-parton correlators. After contraction with the leptonic tensor, the result would be a "twist" expansion of the SIDIS cross section in powers of Λ/Q, where Λ is a scale quantifying the strength of parton-parton correlations inside the proton target and the detected hadron. In this paper, however, we only consider terms of order (Λ/Q) 0 .
In collinear factorization, one chooses k and k to be collinear in 3D space to the momentum of the target nucleon and the detected hadron, respectively. In the (p, q) frame, these read As a consequence, the hadronic tensor turns out to depend only on the 1-dimensional, "collinear" Parton Distribution Function q(x) = dk − d 2 k T φ 2 (k), and Fragmentation Function D q (z) = (z/2) dk + d 2 k T δ 2 (k), where x ≡ k + /p + and z ≡ p − h /k − . At variance with the conventional treatment, we consider generic "average virtualities" v 2 ≈ k µ k µ and v 2 ≈ k µ k µ for the incoming and outgoing partons, whose values will be fixed later rather than put to 0. It is indeed clear from Fig. 1 that k µ k µ ≥ m 2 h , and that this bound cannot be a priori neglected at the kinematics we are interested in. In general, the average virtualities of the quarks entering the diagram in Fig. 1 are determined by the dynamics of the scattering and hadronization processes, see for example Refs. [33,34], and can be in principle different for the scattering and scattered quarks. It is also important to to keep the quark's current mass m q and virtuality v or v conceptually separated. It is only when a quark line is cut, i.e., when the scattered quark appears in the final state, that k µ k µ = m 2 q ; this is clearly not the case for either quark in the handbag diagram of Fig. 1. Furthermore, it is m q , rather than the virtualities v or v (as claimed in Ref. [24]), that appears, via Dirac's equation, in the so-called "equations of motion relations" essential to the treatment of HT terms [28].
Finally, the SIDIS hadronic tensor factorizes into a convolution of a quark PDF, a quark FF, and a hard scattering tensor H µν as [23] 2M where For ease of interpretation and discussion, in this formula we explicitly separated the virtualities v and v from the "massless" partonic momenta k 0 and k 0 defined as In Eq. (8), we reinstated the sum over quark flavors and neglected terms of twist higher than 2 in the twist expansion; a detailed discussion of factorization at twist 3 and twist 4 in inclusive and semi-inclusive DIS can be found in Refs. [28,31,32,35]. Let us remark that, as a result of our calculation and factorization scheme, the trace term appearing in Eq. (9) can be interpreted as the matrix element squared for the scattering of a virtual photon with a parton of momentum k µ 0 and current mass m q = 0; however, when the average virtualities v and v have non-zero values, the δ-functions impose different values of x and z than if this was an actual physical process.
It is interesting to note that, if one does choose v = v = 0, the hard scattering tensor can be rewritten as Then the whole SIDIS hadronic tensor can be interpreted in terms of the parton model scattering and fragmentation of a fictitious free quark of zero mass collinear with the target nucleon and the produced hadron. This model was proposed by Feynman as a heuristically well motivated approximation to the full QCD process in the "infinite momentum" p + ∼ Q frame [36]; however, at sub-asymptotic values of Q, such as those investigated here, the resulting approximation may not be optimal. In the original parton model, the masses of the target and of the detected hadron are neglected. Our Eq. (12), instead, supplements the parton model with mass corrections in a way that was already proposed by Albino et al. in Ref. [37] and shown to provide improved fits of experimental data. If, however, one chooses v = 0 or v = 0, the δ-functions in Eq. (9) cannot be interpreted as expressing four momentum conservation of the fictitious free quark as it scatters on the virtual photon; therefore, the hard scattering tensor cannot be given a parton model interpretation. This lack of intuitive interpretability is not to be considered a show stopper: on the contrary, the hard scattering tensor (9) satisfies by inspection the Ward identity q µ H µν = 0, and therefore the hadronic tensor (8) is a legitimate, gauge invariant approximation of the full scattering diagram in Fig. 1. In fact, as argued in Refs [22,23] and discussed next, v = 0 is a necessary condition to respect 4-momentum conservation in the SIDIS process. Thus, our Equations (8)- (9) provide the means to go beyond the parton model, and to implement this kinematic requirement in a gauge invariant way in the collinear factorization framework.
More in general, the 2-steps procedure discussed above defines an "approximator" of the hard scattering which is analogous to that introduced by Collins, Rogers and Stasto in Ref. [33]. Compared to that one, our approximator takes into account kinematical hadron mass effects, and furthermore allows one to define fully integrated collinear PDFs and FFs instead of the totally unintegrated ones considered in the mentioned reference. For a full proof of factorization, one would furthermore need to verify that this scheme extends at least to NLO, and that our approximator allows one to use Ward identities to factor out and resum longitudinal gluons into the Wilson lines needed to ensure gauge invariance in the operator definition of the PDFs and FFs. The successful phenomenological approach of our LO scheme, to be discussed in detail in Section III, justifies future efforts in this direction.

B. Choice of virtualities
Upon integration over the delta functions, one obtains and, clearly, in the Bjorken limit, one recovers x ≈ x B and z ≈ z h . To proceed further, it is necessary to specify a choice for the average virtualities v and v . For this purpose, we minimally require that the approximated, internal k and k momenta respect four-momentum and baryon number conservation in the factorized diagram, or in other words that the "internal" (approximated) collinear kinematics at parton level matches the "external", hadron-level kinematics. The limits that this requirement places on the possible values of v and v have been derived in detail in Refs. [22,23]. Here we only recall the main results, and defer to Appendix A a subtler point regarding the treatment of baryon number conservation that was not sufficiently explained in those papers. Due to the interpretation of the trace term in the hard scattering tensor H µν as due to a "massless" quark of momentumk 0 scattering on a virtual photon, it is desirable to choose More importantly, if we applied this formalism to semi-inclusive hadron production in electron-positron scattering events, a value v 2 = 0 would be imposed by the cut on the non-fragmenting (light) quark line in the leading order diagram. Thus, Eq. (14) is in fact a necessary condition for the proposed HMC formalism to be universal. Fortunately, a zero value for v is also compatible with the external kinematics [22,23]. As a result, Eq. (13b) requires z = ζ h . The fragmentation of the scattered parton into a massive hadron, instead, requires a non vanishing virtuality v 2 . More precisely, by requiring four momentum conservation in the right hand side diagram of Fig. 1 (i.e., by matching the internal approximated partonic kinematic with the external hadronic kinematics of the fragmentation process), one finds that a parton of light-cone momentum fraction z needs an average virtuality v 2 ≥ m 2 h /z to fragment into a hadron of mass m 2 h . Then, compatibly with the LO constraint on z just discussed, we choose Inserting these virtualities into Eqs. (13a)-(13b), one finds with the Nachtmann-type scaling variables ξ and ζ h defined in Eqs.
(2)-(3). Note that ξ h is reminiscent of the χ = x B (1 + 4m 2 Q /Q 2 ) scaling variable in the ACOT-χ treatment of heavy quarks in DIS [38,39], where an unobserved open heavy flavor of mass 2m Q is produced in the final state, much like the detected hadron of mass m h and momentum fraction z h discussed in this paper. Further exploration of this similarity is left for future work.

C. Factorized hadron multiplicities
Collecting the above results, contracting the hadronic tensor with the leptonic tensor, and accounting for mass corrections also in the inclusive cross section appearing in the denominator [40], one can see that, at finite Q 2 values, the LO hadron multiplicity integrand in Eq. (1) factorizes in terms of products of quark PDFs and FFs, D h q , but evaluated at the scaling variables ξ h and ζ h just derived, and that where J h = dζ h /dz h is a scale-dependent Jacobian factor [23]. This expression is gauge invariant and incorporates the kinematic requirement for the scattered quark to have a non-zero virtuality in order to fragment into a hadron of non-zero mass m h and invariant momentum fraction z h . Note that in the Bjorken limit all masses become negligible (m 2 /Q 2 → 0) and one recovers the usual "massless" M h multiplicity, with its usual parton model interpretation.
The arguments leading to the factorized formula (18) and the proof of its gauge invariance were already laid out in Refs. [22,23], although in a way that seems to have originated some misunderstanding and confusion in the literature [24]. It is our hope that the present discussion will dispel the doubts raised there on the validity of our treatment of hadron mass corrections.

III. NUMERICAL RESULTS FOR KAON MULTIPLICITIES
The HERMES and COMPASS data on integrated kaon multiplicities [19][20][21] appear to be incompatible with each other, a well known and hotly debated fact [26,41,42]. While most of the discussion has centered on kinematic and binning issues, here we argue that HMCs may also play an essential role due to the relatively low average values of Q 2 dominating the HERMES and COMPASS measurements.

A. Data over theory ratios
One way to compare HERMES multiplicities to COMPASS multiplicities is through the ratio between experimental data and theory calculations, in which both the difference in kinematic cuts and Q 2 evolution between the two experiments are canceled, having been included in the theory calculation Eq. (1): with a perfect theory (and in the absence of unaccounted for experimental systematics) the ratio should be equal to 1 within statistical fluctuations.
In Fig. 3, we can observe the data over theory D/T ratio for different sets of PDFs [5,43,44] and FFs [45,46]. The effect of HMCs can be observed comparing the ratio calculated using the massless theory (left panel) and the theory with HMCs (right panel). There clearly is a large FF systematics due to the poor knowledge we have of kaon fragmentation functions, but this amounts largely to an overall multiplicative factor; the PDF systematics is definitely smaller.
In the massless ratios, even looking at only one given FF set, one can notice a difference in size, as well as shape, of the HERMES and COMPASS D/T ratios. Using HMCs the size discrepancy between the two experiments is reduced and the ratio is flatter for both sets of data. In particular the COMPASS ratio is rather flat over the whole x B , while HERMES still persists having a downward slope and a concave shape.

B. Multiplicities in a massless world
A more direct data-to-data comparison of HERMES and COMPASS results, that also reduces the effect of the FF and PDF systematics, can be obtained by defining "theoretical correction ratios" that produce (approximate) massless parton multiplicities at a common beam energy. This method also allows one to interpret the corrected data at face value using parton model formulas such as our M  The first step in the calculation consists in removing the mass effects from the original data using the "HMC ratio" where M h(0) is the massless hadron multiplicity calculated theoretically using Eq. (19) and M h is the finite Q 2 multiplicity from Eq. (18). Using this, the product M h exp × R h HM C can be interpreted as a "massless" experimental multiplicity. In other words, this is the multiplicity that one would expect to measure in a world where nucleons and kaons are massless.
Next, we address evolution effects, i.e., the difference in the Q 2 reach of each x B bin of HERMES and COMPASS. For this, we choose the COMPASS kinematics to be the one at which we want to compare the data. Then, we bring HERMES data to COMPASS energies through an evolution ratio that we define as: Here, the numerator is the massless multiplicity calculated integrating over each one of the HERMES x B bins, but using the kinematic acceptance of the COMPASS experiment; namely, we integrate over the black hatched vertical stripes in the (x B , Q 2 ) phase space shown in the right panel of Fig. 2. The denominator is the massless multiplicity integrated using the original HERMES kinematic cuts (blue vertical stripes in the right panel of Fig. 2). As a result, multiplying the massless HERMES multiplicity found in the previous step by this ratio, we are effectively "evolving" HERMES results to COMPASS energy and spectrometer. Finally, the massless and evolved experimental multiplicities can be defined by multiplying the original data by the appropriate correction ratios; in our case, The correction ratios were evaluated numerically and plotted in Fig. 4 and we find, as expected, that these are relatively stable with respect to the choice of FFs, because the FF systematics shown in Fig. 3 is canceled in the ratios defined in Eqs. (20) and (21). The PDF systematics is also small. Furthermore, hadron mass effects are dominant compared to evolution effects, that are rather small. For COMPASS, the HMC corrections are smaller than at HERMES because the Q 2 accessed at COMPASS is higher at a given x B than at HERMES due to the higher beam energy. The PDF and FF systematic uncertainties are calculated by varying these among the fits listed in Fig. 3, and are typically smaller at COMPASS due to the higher Q 2 reach. The green FF systematic band for the COMPASS R K HM C is very small compared to the HERMES case, and almost invisible in the plot. The purple PDF systematic band for R H→C evo is very small compared to the FF green band.  Left: Parton level multiplicities after applying the theoretical correction ratios given by Eq. (22) to the data shown on the right. The green FF systematics band for COMPASS is very small compared to the HERMES case, and almost invisible in the plot. This is due in part to the larger Q 2 value in each bin, and in part to the additional systematics in the HERMES band introduced by the evolution ratio in Eq. (22b).
In Fig. 5, we plot the experimental K + + K − multiplicity data M K exp on the left and the "massless" multiplicities M K(0) exp on the right using Eqs. (22a)-(22b). In the D/T ratios, HMCs were included in the theoretical calculations; here, instead, HMCs are "removed" from data. After furthermore evolving HERMES data to COMPASS energy (which was automatically achieved in the D/T ratios), the discrepancy in size between HERMES and COMPASS is also largely reduced. Moreover, the corrected data, which can be interpreted directly in terms of parton model formulas, now show for both experiments a negative slope in x B that agrees much better with the (1 − x) β power law behavior of any PDF, including the s-quark. Clearly the slopes and shapes in x B of the HERMES and COMPASS data do not match yet, which indicates that corrections other than HMCs, or unquantified systematic uncertainties, are at play.

IV. NUMERICAL RESULTS FOR KAON MULTIPLICITY RATIOS
Another interesting observable is the K + /K − multiplicity ratio, because one can expect the systematic and theoretical uncertainties in each experiment, as well as Q 2 evolution effects, to largely cancel between numerator and denominator. However, one may still expect some residual mass effect because of the different slopes in z of the K + and K − FFs.
The theoretical correction ratios for the K + /K − multiplicity ratio are plotted in Fig. 6. As expected, the corrections are smaller than in the K + + K − multiplicity sum. The HMCs are non negligible (up to −15% for HERMES and −10% for COMPASS) and of the same order of the HERMES to COMPASS evolution effects. (The FF systematics has not been evaluated because the HKNS fit cannot extract reliable charge separated fragmentation functions.) The original and "massless" data for both HERMES and COMPASS experiments are then plotted in Fig. 7. In this case, the slopes are compatible already in the original data, which shows that much of the systematics difference between the two experiments is not irreducible, but affects only the charged K + + K − multiplicities. However, the discrepancy in size persists. After removing the mass effects and compensating for evolution, the "massless" kaon ratios become fully compatible between the two experiments. A possible exception is the last HERMES x B bin, that shows a sharp change in slope as it also happens for the case of the summed K + + K − , but lies just outside the COMPASS range. In the charged multiplicity case, this could be partly attributed to nuclear binding and Fermi motion effects in the Deuteron target. However, nuclear effects should largely cancel out in the K + /K − ratio, and the origin of the slope change (which is however marginally compatible with the rest of the data within systematic and statistical uncertainties) remains to be understood.

V. OTHER Q 2 -DEPENDENT CORRECTIONS
The HMC calculations presented in the previous sections have been performed at leading twist and leading order in α s accuracy. As these do not seem to exhaust the sources of difference between HERMES and COMPASS integrated multiplicity ratios (although they can potentially explain just by themselves the difference in the kaon ratio measurements) it is worthwhile commenting on other Q 2 -dependent corrections.
Higher-Twist contributions in unpolarized scattering scale as Λ 2 /Q 2 , where Λ is a dynamical non-perturbative scale that quantifies quark-gluon correlations inside the nucleons. Since this is the same kind of scaling exhibited by HMCs, that are however kinematic in origin and scale as m 2 K /Q 2 , one may wonder if HT corrections might explain the residual difference in kaon multiplicities. This is certainly possible, although quantifying those corrections is outside the scope of this article. Here, we just note that the HT phenomenology in inclusive DIS is well developed [44,[47][48][49], while we are not aware of similar studies for SIDIS.
Likewise, one may want to consider NLO corrections, that, however, depend only logarithmically on Q 2 . These may therefore slightly tilt the data/theory ratios for massless multiplicities, but not necessarily close the remaining gap between the HERMES and COMPASS data, as also suggested by the calculations presented in Ref. [50]. This will be explored in a forthcoming paper [51].

VI. DISCUSSION AND CONCLUDING REMARKS
In this paper, we have argued that Hadron Mass Corrections of order O(m 2 /Q 2 ) in SIDIS are non negligible for kaon production at HERMES and COMPASS, where integrated multiplicities have an average Q ranging from 1 to 4 GeV, quite comparable to the kaon mass. These corrections can be captured in a gauge-invariant way at leading twist by new massive scaling variables that incorporate the need for the struck quark to be sufficiently off the free-particle mass shell in order to fragment into a massive hadron. At leading order in the coupling constant, the leading-twist cross section still factorizes into a product of PDFs and FFs, but is evaluated at the Nachtmann variable ξ h of Eq. (16) and the fragmentation scaling variable ζ h of Eq. (3), respectively.
After accounting for HMCs in this way, we found that the discrepancy between the integrated kaon multiplicities measured by the HERMES and COMPASS collaborations is reduced.
For the charge-summed K + + K − multiplicity there are still some differences in slope and shape that need to be investigated. From a theoretical side, one would certainly need to evaluate the effects of Higher-Twist contributions, while NLO corrections, that scale logarithmically in Q 2 , do not seem likely to close the gap remaining between the two experimental measurements.
In the case of the K + /K − multiplicity ratio, where much of the theoretical and experimental systematics can be expected to cancel, the slopes were already similar in the published data, and HMCs can fully reconcile the remaining discrepancy in size. The only possible exception are the last two x B bins of the HERMES measurement, that however lie just outside the reach of the COMPASS experiment. It would therefore be interesting to repeat these measurements at the 12 GeV Jefferson Lab upgrade (JLab 12), where a higher x B range could be covered at Q 2 values comparable to the average HERMES Q 2 , but retaining nonetheless a considerable overlap in x B with both HERMES and COMPASS. Likewise, measuring pion multiplicities at JLab 12 would allow one to investigate the large difference in that overlap region between existing measurements at Jefferson Lab and HERMES noted in Ref. [52], but with an intermediate energy beam.
The nearly perfect agreement in the overlap region of the kaon multiplicity ratios after HMCs are taken into account is a strong indication that the remaining differences in the charged kaon multiplicities are of systematic origin -whether theoretical or experimental remains to be ascertained. This conclusion is strengthened by observing that in the case of the much lighter (and essentially HMC-free) pion, the ratios measured by the two experiments also agree despite displaying strong differences in the charged multiplicity data.
As an outlook, we would like to include deuteron nuclear corrections in our analysis to see if this may explain the large x B behavior of the HERMES data. More importantly, however, we need to prove that factorization extends to NLO in perturbation theory when including a non vanishing average virtuality v 2 = 0 for the fragmenting quark; it will also be necessary to verify that the hard scattering approximator defined in Section II allows one to resum the longitudinal gluons into a Wilson line as it happens in "asymptotic" factorization theorems [27,33,53]. The analogies of our scaling variable ξ h with the χ variable of the ACOT-χ heavy quark scheme also deserve further investigation.
Finally, the results presented in this paper point at the necessity of using hadron mass corrected theoretical calculations in QCD fits of fragmentation function that include HERMES and COMPASS data [54][55][56], in order to avoid deforming the kaon FFs to compensate for the neglected mass effects. Likewise, other power-suppressed corrections such as Higher-Twist terms in SIDIS should also be included, but this is still, to our knowledge, a largely unexplored topic from a phenomenological point of view.
Associates, LLC operates Jefferson Lab and DOE contract No. DE-SC0008791.