Power corrections in semi-inclusive deep inelastic scatterings at fixed target energies

The COMPASS collaboration published precise data on production cross section of charged hadrons in lepton-hadron semi-inclusive deep inelastic scattering, showing almost an order of magnitude larger than next-to-leading order QCD calculations when PhT and zh are sufficiently large. We explore the role of power corrections to the theoretical calculations, and quantitatively demonstrate that the power corrections are extremely important for these data when the final-state multiplicity is low and the production kinematics is near the edge of phase space. Our finding motivates more detailed studies on power corrections for upcoming experiments at Jefferson Lab, as well as the future Electron-Ion Collider.

The COMPASS collaboration published precise data on production cross section of charged hadrons in lepton-hadron semi-inclusive deep inelastic scattering, showing almost an order of magnitude larger than next-to-leading order QCD calculations when P h T and z h are sufficiently large. We explore the role of power corrections to the theoretical calculations, and quantitatively demonstrate that the power corrections are extremely important for these data when the final-state multiplicity is low and the production kinematics is near the edge of phase space. Our finding motivates more detailed studies on power corrections for upcoming experiments at Jefferson Lab, as well as the future Electron-Ion Collider.

I. INTRODUCTION
Unveiling the structure of nucleons (or hadrons, in general) in terms of quarks and gluons of quantum chromodynamics (QCD) is one of the central goals that have been actively pursued by the science community since the first lepton-proton deep inelastic scattering (DIS) experiment taken place at SLAC about 50 years ago [1]. However, owing to the fact that no isolated quarks and gluons have ever been seen in a modern detector, nearly all analyses of high energy scattering events with identified hadron(s) rely on the QCD factorization theorem [2] that provides the link between the observed hadrons and the quarks and gluons, or collectively, partons, that participated in hard scatterings. The inclusive lepton-proton DIS at a large momentum transfer, Q Λ QCD ∼ 1/fm, is dominated by the scattering of the lepton off one active quark/parton inside the colliding proton [3]. With one hard momentum scale Q, the inclusive DIS cross section is not very sensitive to the dynamics at a typical hadronic scale ∼ 1/fm Q, and can be factorized into the lepton-quark scattering at a short-distance (∼ 1/Q) multiplied by corresponding quark parton distribution functions (PDFs), φ i/P (x, µ 2 ), interpreted as the probability distribution to find this active quark of flavor i inside the colliding proton at the factorization or probing scale µ ∼ Q, carrying the proton's momentum fraction x. This factorization is known as the QCD collinear factorization, with active quark's transverse momentum k T integrated into the PDFs and overall corrections suppressed by inverse powers of Q. The measurement of the inclusive DIS cross section has provided good information on the proton's partonic structure, encoded in these factorized PDFs.
It is the QCD factorization that provides the "probe" -the short-distance partonic scattering to enable us to "see" quark, gluon and their dynamics indirectly. The predictive power of such factorization approach relies on both the precision of the probe, which we could achieve and improve by calculating the partonic interactions at the scale Q order-by-order in QCD perturbation theory, and the universality of these PDFs, so that we are able to extract them from data of some experiments and use them to predict and to be tested in other measurements. In terms of this factorization formalism, QCD has been extremely successful in interpreting almost all available data from high energy scatterings with probing distance less than 0.1 fm (or equivalently, with the momentum transfer greater than 2 GeV) [4,5]. It is this success that has provided us the confidence and the tools to discover the Higgs particle and to explore new physics beyond the Standard Model of particle physics in high energy hadronic collisions [6][7][8].
Instead of summing over all hadronic final states, the semi-inclusive DIS (SIDIS) identifies one final-state hadron of momentum P h , as illustrated in Fig. 1, and covers a part of the inclusive lepton-hadron DIS cross section. By detecting one hadron in the final-state, SIDIS enables us to explore the emergence of color neutral hadrons from colored quarks and gluons, in addition to the information on finding a quark or gluon inside the colliding hadron. By selecting different type of observed hadrons (pion, kaon, ...), SIDIS provides opportunities to study the flavor dependence of QCD dynamics. With a large momentum transfer carried by the virtual gauge boson of momentum q (Q 2 ≡ −q 2 1/fm 2 ), as shown in Fig. 1, SIDIS could provide a short-distance probe with an additional and adjustable momentum scale by measuring the hadron at different momentum P h . In the frame where the exchange virtual gauge boson of momentum q and the colliding hadron of momentum P are headed on, the leading contribution to the SIDIS is naturally from the region where the transverse momentum of the observed hadron, P h T Q, and the scattering provides a short-distance probe with two very different momentum scales, from which the harder scale Q localizes the hard collision to "see" the particle nature of quarks and gluons, while the soft scale P h T is sensitive to the confined motion of quarks and gluons in the direction perpendicular to the direction of the colliding proton. In this kinematic regime where Q 2 P 2 h T 1/fm 2 , similar to the inclusive DIS, SIDIS cross section can be factorized into a product of perturbatively calculable lepton-parton scattering at the hard scale Q, corresponding transverse momentum dependent (TMD) parton distribution functions (or simply, TMDs), φ i/P (x, k T , µ 2 ) with k T being the active parton's transverse momentum perpendicular to the direction of the colliding hadron of momentum P , and TMD fragmentation functions (FFs), D j→h (z, p T , µ 2 ) with the emergent hadron-type h carrying momentum fraction between z and z + dz of the fragmenting parton of momentum p and p T being the parton's transverse momentum off the direction of the observed final-state hadron of momentum P h , where i, j = {q,q, g} represent the active parton flavors [9,10]. In terms of this TMD factorization formalism, with the corrections of O(P h T /Q), lepton-hadron SIDIS is an excellent process to probe three-dimensional (3D) confined motion of quarks and gluons inside a bound proton, and has been actively pursued by experimental programs at all lepton-hadron scattering facilities, such as COMPASS at CERN [11] and various experiments at Jefferson Lab [12], as well as the future Electron-Ion Collider (EIC) [13].
On the other hand, when P h T ∼ Q Λ QCD , the SIDIS is dominated by a single hard scale, no longer sensitive to the dynamics at the scale of active parton's transverse momentum k T (or p T ), which is typically much less than Q, and the cross section should be better described by QCD collinear factorization approach. At the leading power (LP) of momentum transfer, the collinearly factorized SIDIS cross section is given by where Bjorken variables are and dσ l+i→l +j+X is the perturbatively calculable partonic hard part for the lepton to scatter off a collinear on-shell parton of momentum xP and flavor i to produce an active parton of momentum p = P h /z and flavor j, which fragments into the observed hadron of momentum P h . The φ i/P (x) and D j→h (z) in Eq. (1) are collinear PDFs and FFs, respectively, and their factorization scale dependence is suppressed. The corrections to Eq. (1) are suppressed by 1/P 2 h T (or 1/Q 2 ) for spin-averaged cross sections. A smooth transition and matching between the TMD and collinear factorization approaches has been developed and tested for various two-scale observables, e.g., heavy gauge boson W/Z production in hadron-hadron collisions [14][15][16]. However, the precise SIDIS data, recently published by COMPASS Collaboration [11], show almost one order of magnitude discrepancy between the data and the leading order (LO) theoretical prediction for the region where P h T ∼ Q and z h are large [17], and the next-to-leading order (NLO) corrections provide very little help to reduce this discrepancy [17,18]. Soft gluon resummation does help enhance the prediction of theoretical calculations for the high P h T region, but, not enough to make up the order of magnitude difference [19,20]. Although the P h T -integrated cross section, which is dominated by small-P h T regime, has been well described within the collinear factorization, a successful description of the full P h T -distribution is critically important for understanding QCD dynamics and extracting the multidimensional partonic structure of the nucleon.
The QCD factorization for high energy scattering cross sections with identified hadron(s), both TMD and collinear, is an approximation that neglects corrections suppressed by powers of the large momentum transfer of the hard collisions. When z h → 1 in SIDIS, as shown in Eq. (1), the phase space for the fragmenting parton to shower (to radiate softer partons) is vanishing. Consequently, the fragmentation function, D j→h (z) ∝ (1 − z) n → 0 at large z with n ≥ 1 for meson production (or a larger power for baryon production), strongly suppresses the probability for the produced active parton to become the observed hadron, and gives a power-like suppression to the LP contribution to the SIDIS cross section in Eq. (1). If next-to-leading power (NLP) contributions to the cross sections are factorizable, corresponding perturbatively calculable hard parts will be suppressed by inverse powers of the large momentum transfer, such as 1/Q α or 1/P α h T with α positive, in comparison with the hard parts of LP contributions. However, the overall factorizable power suppressed contributions to the measured cross section could be sizable or even more important if the active partonic states, produced through the power suppressed short-distance partonic subprocesses, are much more likely to become the observed hadrons than the partonic states produced at the LP. One example is the heavy quarkonium production at large transverse momentum at collider energies [22][23][24], where producing a heavy quark pair at large P h T is power suppressed comparing to the production of a single gluon at the same P h T , but, the probability for the produced heavy quark pair to become a bound quarkonium could be much greater than the probability for the produced gluon to become a heavy quarkonium in some kinematic regions.
For the production of a positively charged hadron in SIDIS, such as π + , producing a quark-antiquark pair at high P h T , as sketched in Fig. 2 (right), is certainly suppressed in comparison with the production of a single quark or antiquark at the same P h T , as shown in Fig. 2 (left). However, a quark-antiquark pair with the right quantum number, e.g., ud for π + , could be much more likely to become the measured meson than a single quark or antiquark through the fragmentation process when the phase space for the radiation is vanishing and the multiplicity for the events is low. For example, when z h → 1, the factorized LP contribution to π + production in Eq. (1) is suppressed by powers of (1−z) from the FFs with z ∼ z h , and the factorizable NLP contribution, as we show in this paper, could have the leading transition from ud to π + to be proportional to δ(1 − z) without the power suppression in (1 − z). While the production of the ud pair is suppressed by inverse powers of P h T , as we will demonstrate below, the trade off between the 1/P 2 h T suppressed hard parts at the NLP and the power suppressed FFs at the LP could make the formally power suppressed contributions to the SIDIS cross section very important for low multiplicity events.
In this paper, we investigate the NLP corrections to SIDIS production of charged mesons near the threshold, where P h T ∼ Q and z h → 1, in terms of QCD collinear factorization approach. Instead of calculating all possible corrections at the NLP, we focus on the partonic subprocesses that could have the best chance to compete with the better studied LP contribution, and more specifically, we calculate the LO perturbative contribution in α s to the production of a quark-antiquark pair that have the right flavor combination of the observed mesons. In order to demonstrate the impact of the power corrections, quantitatively, we estimate the leading quark-antiquark pair FFs to a charged meson, which is proportional to δ(1 − z), in terms of the better-known light-cone meson distribution amplitude square, by neglecting contributions suppressed by powers of (1 − z). We find that the NLP corrections to SIDIS are extremely important for the production of charged mesons when the final-state multiplicity is low and the production kinematics is near the edge of phase space. Our finding warrants a much more detailed study of power corrections to the SIDIS cross section near the edge of phase space where the hadron multiplicity is very low, which provides a unique opportunity to explore the mechanism of hadronization and color neutralization in QCD -the emergence of hadrons from produced quarks and gluons in high energy collisions.
The rest of this paper is organized as follows. In Sec. II, we introduce the factorization formalism for SIDIS to the accuracy of NLP, provide a leading order calculation of the short-distance hard parts at the NLP from the channels in which a quark-antiquark pair is produced with the same flavor combination of the valence content of the observed meson, and derive an approximate relation between the quark-antiquark FFs to a meson and the square of distribution amplitudes of the same meson. In Sec. III, we show our numerical estimation of the size of the power corrections in comparison to the size of the LP contribution, quantitatively. Contributions from channels other than the direct production one are discussed in Sec. IV. Finally, conclusions and outlooks are given in Sec. V.

II. NEXT-TO-LEADING POWER CONTRIBUTION TO SIDIS
With effectively one large momentum scale, P h T ∼ Q Λ QCD observed, SIDIS cross section with a large P h T hadron could be factorized in terms of the QCD collinear factorization approach. With two identified hadrons, the initial-state hadron and the observed finalstate hadron, only the leading and the first subleading power contributions to the SIDIS cross sections, in terms of the inverse power expansion of the observed large momentum transfer, are perturbatively factorizable, similar to the collinear factorization for inclusive Drell-Yan cross sections [21]. In general, the sub-leading power contributions are much smaller comparing to the LP contributions because of the power suppression of the observed large momentum scale for the short-distance hard part, unless they can get enhancement from the hadronization [22][23][24], sufficient corrections to a steeply falling spectrum near the edge of available phase space [25,26], or the multiple scattering in a large size and/or dense medium [27][28][29]. Here, we consider possible enhancements from both the hadronization and the steep falling spectrum near the edge of phase space.

A. The factorization formalism
With the approximation of one-photon exchange, as shown in Fig. 1, the leptonic contributions to the SIDIS cross section is well-understood and well-defined. In the rest of this paper, we present our calculations in terms of scattering of a virtual photon γ * of momentum q with Q 2 = −q 2 > 0 on a hadron A of momentum P to produce a charged hadron h of momentum P h : γ * (q) + A(P ) → h(P h ) + X. The corresponding formalisms can also cover the situation of photoproduction with a real photon Q 2 = 0. When the photon is either deep virtual with Q 2 = −q 2 ∼ P 2 h T or real with Q 2 = 0, the scattering cross section can be expressed in terms of a collinearly factorized formalism [21,23], where a, f (and f ) run over all parton flavors: q,q, and g, κ runs over all spin and color states of [f f ], φ a/P (x) and D f →h (z) are defined above as PDFs and FFs, respectively, with l = xP and P h = zp, and Fig. 3, with the momentum fraction of the pair carried by the observed hadron, z = P h /p, and ξ and ζ being relative momentum fractions of the two partons in the amplitude and its complex conjugate, respectively [22][23][24]. 3. Sketch of the fragmentation function for a quarkantiquark pair of momentum p to fragment to a hadron of momentum P h .
, respectively, are perturbatively calculable short-distance hard parts for producing a single parton of flavor f and a pair of partons of flavor combination [f f ], while the hard part for producing a pair partons is powerly suppressed by 1/P 2 h T in comparison with the hard part for producing a single parton at the same momentum p. In Eq. (3), we neglected other NLP terms proportional to twist-4 parton distributions and the same single parton FFs, φ [ab]/P (x, x , x ) ⊗ D f →h (z), since corresponding partonic hard part is powerly suppressed comparing the LP hard part, while no potential enhancement from the hadronization since the same FFs are used. There could be additional contributions to the factorized formalism in Eq.(3), and they are typically suppressed by even higher power in 1/P 2 h T and/or 1/Q 2 . These contributions are not expected to be factorizable. The partonic short-distance hard parts in Eq. (3) at the leading order in power of strong coupling constant α s are given by (5) where parton level Mandelstam variables are defined aŝ with the constraint imposed by the phase space δ-function of the massless two-particle final-state and the momentum conservation of the 2 → 2 partonic subprocess. In Eqs. (4) and (5), the denominator 2(ŝ + Q 2 ) = 2E γ * 2E a |v γ * − v a | is the flux factor of the partonic scattering process, and |M γ * +a(l)→f (p)+X | 2 and |M γ * +a(l)→[f f (κ)](p)+X | 2 , respectively, are squared LP and NLP partonic scattering amplitudes with initial-state spin and color averaged and final-state spin and color summed.

B. Power suppressed partonic hard parts
The theoretical prediction (or an estimation) for the size of power corrections in Eq. (3) relies on our ability to calculate the short-distance hard parts to produce a pair of partons and our knowledge of the nonperturbative double-parton FFs. For getting the most enhancement from the pair's hadronization to compensate the power suppression of the partonic hard parts to produce the pair at NLP, we calculate the partonic hard parts from subprocesses that can produce a quark-antiquark pair with the same valence quark flavor combination as the produced charged meson.
At the lowest order in power of α s , we need to consider the subprocess: γ * + q → [qq (κ)] + q , with corresponding Feynman diagrams given in Fig. 4 where the quark and antiquark flavors q andq should match the valence flavors of the measured meson. The color state of the produced quark-antiquark pair can be either a singlet " [1]" or an octet " [8]" state, and corresponding color projection operators are respectively proportional to δ ab and t A ab , where a, b = 1, 2, ..., N c are color indices for the quark and the antiquark, and t A ab is the generator of the SU (N c ) color group in the fundamental representation with A = 1, 2, ..., N 2 c − 1 and N c = 3 in QCD. In this paper, we use the same color projection operators for the hard part as done in Refs. [23,24], where d, c are color indices for the quark-antiquark pair in the complex-conjugate of the scattering amplitudes. The corresponding color projection operators for the quarkantiquark FFs in Fig. 3 are [23,24] C [1] ab,cd = C [8] ab,cd = 1 They satisfy the normalization condition, where I, J = [1], [8].
The spin projection operators for the four spin states of the quark-antiquark pair can be given by (γ · p) ij , (γ · p γ 5 ) ij , and (γ · p γ α ⊥ ) ij with α = 1, 2, or their linear combinations. They could be referred to as the vector (v), axial-vector (a), and tensor (t) projections. Following Refs. [23,24], we choose the spin projection operators for the hard part as and corresponding spin projection operators for the quark-antiquark FFs to be where n is a null vector with n 2 = 0, defined to be conjugated to the momentum of the quark-antiquark pair p, such that p · n is the only nonvanishing component of p µ if p 2 = 0. As required, the spin projection operators satisfy the normalization condition, where s, s = v, a, t.
The spin state of a virtual photon of momentum q can be either transversely polarized with the polarization vectors µ ± or longitudinally polarized with the polarization vector µ L . The transverse spin polarization tensor is defined as and the longitudinal spin polarization tensor is given by where v andv are two null vectors introduced to pick the "+" and "−" light-cone components of the photon momentum q with v 2 =v 2 = 0, v ·v = 1, and The color factors for all the squares of diagrams in Fig. 4 are the same for each color projection: where the factor 1/N c is from the average of initial state quark colors. Thus we can factor it out from amplitude squares of the diagrams in Fig. 4. Now, we calculate the invariant scattering amplitude squares from the diagrams in Fig. 4 for the production of a quark-antiquark pair in the axial-vector spin state. For example, the invariant amplitude square of the diagram (a) in Fig. 4 with transversely polarized photon is given by where the factor 1/2 is from the spin average of the initial state quark, C = C [1] (or C [8] ) is the color factor, d µν (q)/2 projects out one transversely polarized state of the colliding photon of momentum q, and Feynman gauge was used for the gluon propogators. For calculating contributions from the scattering of a longitudinally polarized virtual photon, one only needs to replace d µν (q)/2 above by K µν (q) defined in Eq. (21). We obtain the invariant amplitude squares of all diagrams in Fig. 4 with a transversely polarized photon and the produced quark-antiquark pair in an axial-vector spin state, where the parameters are defined as A η =t + η(ŝ +û), with η = ξ, ζ,ξ,ζ andξ = 1 − ξ,ζ = 1 − ζ, respectively. Similarly, for a longitudinally polarized photon, we have where all parameters are the same as those defined in Eqs. (29) -(31).
By replacing the axial-vector spin projection for the produced quark-antiquark pair by a vector spin projection, we find that the invariant amplitude squares of all diagrams in Fig. 4 have the same results as those from Eqs. (25) to (35). This is easy to understand at this order of calculation since we neglect the quark mass and the two γ 5 in the spin trace, as those in Eq. (24), can be combined into γ 2 5 = 1.

C. The quark-antiquark fragmentation function
In order to quantitatively estimate the size of the power corrections to the semi-inclusive production of large P h T hadrons, in terms of the factorization approach in Eq. (3), we need the knowledge of quark-antiquark FFs in addition to the calculation of perturbative partonic hard parts. For getting the most enhancement from the fragmentation of the quark-antiquark pair, we focus on the production channels in which the flavor of the produced quark-antiquark pair [f f ] matches the flavor combination of the valence components, [qq ], of the produced meson h. As introduced in Ref. [23], the quark-antiquark FFs D [qq (κ)]→h (z, ξ, ζ) are defined as where C and P are respectively the color and spin projection operators defined in Sec. II B with color and spin indices suppressed. The Φ n (y − ) in Eq. (36) is the gauge link in the fundamental representation of QCD color group, given by where P and G A represent the path ordering and the gluon field, respectively, and t A is the generator of SU(3) color with color index A, as introduced in Sec. II B.
Like the single-parton FFs, the quark-antiquark FFs are nonperturbative and cannot be calculated within the QCD perturbation theory, while their factorization scale dependence could be calculated if the variation is within the perturbative regime [23]. In order to estimate the size of the power corrections, quantitatively, we make the following approximation. We assume that at an input scale µ 0 , the quark-antiquark pair FFs are dominated by the final-state in which there is no additional hadron produced other than the observed meson, which mimics the physical condition when the observed meson is produced near the edge of phase space with very low multiplicity. That is, the hadronic final-state in Eq. (36) is approximated as |h(P h )X ≈ |h(P h ) . Under this approximation, as shown below, we can relate the quark-antiquark FFs to the square of meson distribution amplitude, to help us to estimate the size (possibly a lower limit) of the power corrections quantitatively by using the better known knowledge on the meson distribution amplitudes.
For pseudoscalar mesons, e.g., pions and kaons, the distribution amplitude φ h (x, µ) has the following matrix element definition, where the contraction of color and spin indices are explicitly shown, In deriving the above result, the relation in Eq. (38) was used. Under this extreme approximation, |h(P h )X ≈ |h(P h ) , the nonperturbative quark-antiquark pair FFs can be expressed in terms of a product of two nonperturbative meson distribution amplitudes φ h (ζ, µ 0 ) and φ h (ξ, µ 0 ), which have been better studied than the quark-antiquark pair FFs. We will come back to discuss the corrections to our extreme approximation in Sec. IV.

III. COMPARISON BETWEEN LP AND NLP CONTRIBUTIONS
In order to understand the relevance and potential impact of the NLP corrections, we evaluate and compare the size of the LP and NLP contributions to the SIDIS production of a charged meson at large transverse momentum, numerically, in this section.
The differential multiplicities for charged hadrons in lepton DIS off a deuteron target were recently measured by the COMPASS Collaboration [11]. The differential multiplicity is defined as the ratio between the SIDIS and the inclusive DIS differential cross sections: where x B and z h are defined in Eq. (2), and P h T is defined in the photon-target frame. Since the produced hadrons are dominated by pions, we only calculate SIDIS cross sections for charged pions, π ± in this section, instead of the sum of all long-lived charged hadrons, h ± , as included in the COMPASS data.
Since the purpose of this paper is to show the relevance and potential impact of NLP contributions to warrant a urgent and more detailed study of the NLP power corrections to the low multiplicity observables in SIDIS, instead of a precise fitting to the data, we perform straightforward leading order calculations in power of α s for both LP and NLP contributions to the differential multiplicity defined in Eq. (41). Since the inclusive DIS cross section in the denominator is dominated by the low P h T region and consistent with the LP contribution alone, we include both LP and NLP contributions to the SIDIS cross section in the numerator while having the LP contribution to the inclusive DIS cross section in the denominator in our numerical evaluation of the differential multiplicity below. For PDFs, we use CT14 PDF set in our numerical evaluation [30]. For single-parton FFs in the LP contribution, we use the NNFF1.0 FF sets [31]. For the NLP contribution, we only consider the channels with produced quark-antiquark pair matching the valence flavors in the color singlet and axial-vector spin projection. In addition, we approximate, |h(P h )X ≈ |h(P h ) , in the definition of quark-antiquark FFs, as discussed in Sec. II C, and express the quark-antiquark FFs in terms of two distribution amplitudes, as in Eq. (40). For the distribution amplitude, we adopt those from Ref. [32]. Our numerical results are shown in Fig. 5 along with COMPASS data [11]. Consistent with what was found in Refs. [17,18], the LP contribution alone (the dotted curves) is about one order of magnitude smaller than the data. While adding next-to-leading order corrections to the LP contribution does not help much [18], it is clear from Fig. 5 that the NLP contribution (the difference between the solid and dotted curves) is large, and could be as large as a factor of five of the LP contribution when z h and P h T are large, near the edge of phase space. Therefore, in this regime where the multiplicity is low and there is not much phase space for radiation (into light hadrons), it is very important to include the NLP corrections in the QCD global fitting for extracting PDFs and FFs. It is also an opportunity for studying QCD power corrections and the formation or emergence of hadrons from perturbatively produced quarks and gluons.

IV. DISCUSSIONS AND FUTURE OPPORTUNITIES
As emphasized earlier, it is not our goal of this paper to fit the COMPASS data to extract the NLP contributions, since the precise LP and NLP contributions to one physical observable, or more precisely, to the differential multiplicity in Eq. (41), depend on more than one unknown, nonperturbative function. In principle, we need theoretical calculations for more physical observables, which are also sensitive to the same quarkantiquark FFs, and corresponding data to perform QCD global analyses to extract both PDFs and FFs, as well as these new quark-antiquark FFs, which could provide much more insights to the color neutralization and formation of light hadrons, complimentary to what we have learned from the LP single-parton FFs. The predictive power of this QCD factorization approach beyond the LP is our ability to calculate the short-distance hard parts and the universality of these new multi-parton FFs.
In this section, we will discuss the source of possible contributions to these new quark-antiquark FFs to gain some insights into their potential structure and functional forms, and to identify new physical observables that could also be sensitive to the same quark-antiquark FFs, so that we could test the universality of these new multi-parton FFs and QCD dynamics beyond the LP contributions in the future work.
FIG. 6. Feynman diagram representation of the FFs for a pair of ud to fragment into a π + meson.
The quark-antiquark FFs are non-perturbative functions and cannot be calculated within QCD perturbation theory. However, like PDFs and FFs, we might be able to gain some insights into these nonperturbative functions' asymptotic behavior as the variables of these functions approach to an extreme limit, such as z → 1 (or x → 1 or 0 in the case of PDFs). With the operator definition in Eq. (36), in principle, we could represent these quark-antiquark FFs in terms of Feynman diagrams -a Feynman diagram representation. For example, the FFs for a ud pair to fragment into a π + could be represented by an infinite number of Feynman diagrams, as shown in Fig. 6. The first diagram on the right of the "≈" sign is effectively the lowest order diagram in power of α s in the approximation, |h(P h )X ≈ |h(P h ) , which led to the approximated expression of D [qq (1a)] (z, ξ, ζ, µ 0 ) in Eq. (40). With additional radiation of gluons, other diagrams in Fig. 6 could also contribute to D [qq (1a)] (z, ξ, ζ, µ 0 ), but, cannot be proportional to δ(1 − z), instead, proportional to (1 − z) n as z → 1. Although the power of n is a nonperturbative number and depends on the scale at which the FFs are measured, the power n should be positive that leads to a power-like suppression to the NLP contribution from these diagrams, similar to the suppression from LP single parton FFs when z → 1 as discussed earlier in this paper. In addition, a quark-gluon pair could also fragment into a meson as illustrated by Feynman diagrams in Fig. 7, which is suppressed by the power of 1 − z as z → 1. In general, in a confining theory, like QCD, the neutralization of color of these radiated gluons (and/or quarks) requires them to turn into physical hadrons in the final-state, such as the pions with the lightest mass, which strongly suppresses their contributions to the physical cross sections near the edge of phase space (or those with a very small multiplicity).
Feynman diagram representation of the FFs for a ug pair to fragment into a π + meson.
Although the factorized NLP contribution to the SIDIS cross sections in Eq. (3) is formally suppressed by the hard scale P h T , the impact of the NLP contribution to the physical cross sections does not vanish as the power of 1/P h T [23,33,34]. With the factorization formula in Eq. (3), which is a factorization of perturbative collinear singularities of partonic scattering, we must modify the DGLAP evolution of LP single-parton FFs to be consistent to the collinear factorization accuracy at the NLP. Following the discussion in Ref. [23], we can derive the evolution equations for both the single-parton and double-parton FFs from the factorization formalism in Eq. (3). Since a physical observable should be independent of the choice of the factorization scale, we have where ⊗ represents the convolution of momentum fractions as defined in Eq. (3). From Eq. (42), we obtain a closed set of evolution equations for the FFs [23] where is the evolution kernel for resumming logarithmic collinear contribution to the doubleparton FFs, γ f →f is the normal LP DGLAP-type evolution kernel for resumming logarithmic collinear contribution to the one-parton FFs, and γ f →[f f (κ )] is a newtype of evolution kernel for resumming the collinear contributions from the diagrams such as those in Fig. 8 to the one-parton FFs. Although the second terms on the right-hand-side of Eq. (44) is power suppressed by 1/µ 2 , its contribution to the physical observables, such as the SIDIS cross section in Eq. (3), does not vanish as 1/P h T since it contributes to the slope of D f →h , not the D f →h itself [23,33,34]. That is, in order to understand the true impact of the NLP contribution to the SIDIS, we need to do a simultaneous QCD global fitting of PDFs and FFs [35], together with double-parton FFs if one wants to include the COMPASS data or other data near the edge of phase space.
To avoid the difficulty of having too many doubleparton FFs, as the leading approximation, it might be practical for now to make the "valence quark" approximation to keep the quark-antiquark flavors -[f f (κ)] in both the cross section calculations and evolution equations to be the same as the valence quark flavors of the observed meson. That is, for the double-parton fragmentation functions, we ignore the contributions from the double-parton FFs in Fig. 7, while keeping the contributions from the diagrams in Fig. 6.
To close this section, we estimate the NLP contribution to charged pion and kaon productions in upcoming SIDIS experiments at Jefferson Lab (JLab), which have much lower collision energies than what COMPASS had, and thus, should have much less high multiplicity events. Therefore, the impact of the NLP contribution at JLab kinematics could be more significant than that at COMPASS kinematics. In Fig. 9, we present our calcu- lated differential multiplicities, defined in Eq. (41), for both pion and kaon production in SIDIS experiments with a typical kinematics at JLab. We show the LP (dashed), NLP (dot-dashed) and LP+NLP (solid) contri-butions, respectively. As expected, the NLP term dominates large-P h T region, due to the strong (1 − z)-power suppression from the single-parton FFs to the LP contribution, even though the NLP contribution is formally suppressed by extra power of 1/P h T . To make this point even more quantitative and transparent, we plot the fractional contribution to the differential multiplicity from the LP (dashed) and NLP (dot-dashed) in the lower panels, respectively. In addition to the P h T dependence, the NLP contribution has different rapidity distributions from that of the LP contribution. In Fig. 10, we show the LP, NLP and LP+NLP contributions separately as a function of the rapidity y h of the measured charged meson. The rapidity y h is defined in the photon-target collinear frame with virtual photon momentum q = (−Q/ √ 2, Q/ √ 2, 0 ⊥ ). As expected, the NLP contribution favors more negative rapidity, which corresponds to larger z h regime where the LP contribution is suppressed by powers of (1 − z) from its single-parton FFs. Because the phase space for the produced partons to radiate to light hadrons is smaller at JLab energies, the NLP contribution is sizable, about 20%, even in the mid-rapidity region comparing with the LP contribution. Therefore, more detailed studies of the NLP corrections are urgent for upcoming SIDIS experiments at JLab.

V. SUMMARY AND CONCLUSIONS
We have presented the first calculations of power corrections to charged meson productions in SIDIS at large transverse momentum, P h T ∼ Q Λ QCD , in terms of the QCD collinear factorization formalism. We found that the power corrections are very important for the events near the edge of phase space where the hadron multiplicity is low, and z h → 1 and P h T is large.
By expanding the SIDIS cross section in terms of the inverse power of the observed large momentum scale, 1/P h T or 1/Q, we perturbatively factorized the leading and the first subleading power contributions into shortdistance hard parts and corresponding nonperturbative PDFs and FFs. The first subleading power (or the NLP) contributions are formally suppressed by powers of the large momentum scale in the hard part comparing to the LP contributions, but, as we have shown in this paper, the net size of its contribution to the cross section is not necessarily smaller than the formal LP contribution, since the hadronization for the produced quark-antiquark pair might be more effective when the pair has the right quantum number to match to the measured meson. We demonstrated that contributions from such more direct production channels are very important for cross sections near the edge of phase space where events are typically with a low multiplicity. When the phase space for parton shower is closing out at the edge of the kinematic limit, it makes much harder for the fragmenting quark of the LP contribution to neutralize its color, which is consistent with the (1 − z)-type power suppression from the single-parton FFs when z → 1. On the other hand, for the NLP contribution, the leading transition for the produced quark-antiquark pair to a measured meson could be proportional to δ(1 − z) without the phase space suppression, which could help make up the 1/P h T -type suppression to produce the quark-antiquark pair at the NLP instead of a single fragmenting parton at the LP.
To quantitatively estimate the size of the NLP contributions, possibly its lower limit, we performed a calculation considering only the partonic subprocesses that produce a quark-antiquark pair with the same flavor combination of the valence constituents of the detected meson, because these subprocesses could have the best chance to compete with LP channels. From the operator definition of the quark-antiquark FFs, and an approximation for the hadronic final-state of the FFs, |h(P h )X ≈ |h(P h ) , we were able to express the unknown quark-antiquark FFs to a charged meson in terms of this meson's distribution amplitudes, which are better studied. With the perturbatively calculated LO hard parts and the approximated quark-antiquark FFs, we were able to present numerically the approximate size of the NLP contributions, and in particular, to demonstrate quantitatively how important the NLP contributions are in comparison with the LP contributions for the COMPASS kinematics, as well as the energy regime at JLab. At large z h and P h T , we found that the NLP contribution could be as large as five times of the LP contribution. Although the purpose of this paper is not to fit COMPASS data, our finding ensures the importance to take into account NLP corrections in the QCD global analysis if the COMPASS data or other data near the edge of phase space are going to be included in the analyses. Therefore, it is urgent to have more detailed studies on the NLP corrections.
The new double-parton FFs could provide more insights to the color neutralization and the formation of light hadrons in complimentary to the knowledge we have learned from the single-parton FFs. As nonperturbative functions, they cannot be calculated within the QCD perturbation theory, and a simultaneous QCD global fitting of PDFs, single-parton FFs, and double-parton FFs, together with more physical observables that are sensitive to double-parton FFs is needed. The DGLAP evolution equation of the single-parton FFs should also be modified consistently to the accuracy at the NLP. Though the correction term from the NLP is suppressed by 1/µ 2 in the evolution equation, its effect on physical observables does not vanish even at high scales, because it contributes to the slope of FFs [23,33,34]. The predictive power of such factorization approach is the universality of the new FFs, together with our ability to calculate the short-distance hard parts. Upcoming experiments at JLab, as well as those at the future EIC, will provide ample opportunities to study the multi-parton FFs to investigate a new domain of QCD dynamics sensitive to multi-parton correlations.