Matching of Fracture Functions for SIDIS in Target Fragmentation Region

In the target fragmentation region of Semi-Inclusive Deep Inelastic Scattering, the diffractively produced hadron has small transverse momentum. If it is at order of $\Lambda_{QCD}$, it prevents to make predictions with the standard collinear factorization. However, in this case, differential cross-sections can be predicted by the factorization with fracture functions, diffractive parton distributions. If the transverse momentum is much larger than $\Lambda_{QCD}$ but much smaller than $Q$ which is the virtuality of the virtual photon, both factorizations apply. In this case, fracture functions can be factorized with collinear parton distributions and fragmentation functions. We study the factorization up to twist-3 level and obtain gauge invariant results. They will be helpful for modeling fracture functions and useful for resummation of large logarithm of the transverse momentum appearing in collinear factorization.


Introduction
Semi-Inclusive DIS(SIDIS) have three different kinematical regions in which different QCD factorizations are used. These regions are classified according to the momentum of the observed hadron in the final state. In the region where the transverse momentum is large, the production rate can be predicted with standard QCD collinear factorization. In the region where the observed hadron moves in the forward region of the virtual photon with a small transverse momentum, one can use Transverse-Momentum-Dependent QCD factorization [1,2], which involves TMD parton distributions and TMD parton fragmentation functions. The third region is called as target fragmentation region, where the observed hadron moves in the forward region of the initial hadron with a small transverse momentum k ⊥ . In this region, the standard collinear factorization fails because of that the perturbative part in the factorization becomes divergent as powers of ln k ⊥ .
Experimentally, the production in the target fragmentation region has been observed in HERA experiment [3]. It has been suggested in [4] that the production rate can be factorized with fracture functions. An one-loop study of SIDIS in [5] shows that the factorization holds at the order. One can prove the factorization with fracture functions at all orders, as discussed in [6,7]. For the production of one hadron combined with a lepton-pair in hadron collisions, factorizations with fracture functions have been shown to hold at one-loop level in [8,9], where the produced hadron is in the forward-or backward regions. However, the factorization for hadron collisions can be failed as discussed in [10]. Fracture functions, also called as diffractive parton distributions, are parton distributions of an initial hadron with one diffractively produced particle observed in the final state. Currently, most information about fracture functions comes from analysis of HERA data, e.g., in [11,12]. The production in target fragmentation region can be studied with current experiments at JLab and planned experiments at Eic [13] and EicC [14]. Hence, more information of fracture functions will be available.
It is noted that the factorization in the target fragmentation region holds for Q ≫ k ⊥ where Q is the square root of the momentum transfer of the virtual photon of SIDIS. In the case of Q ≫ k ⊥ ≫ Λ QCD , one can also use the collinear factorization to factorize the differential cross-sections with collinear parton distributions and fragmentation functions. Therefore, for k ⊥ ≫ Λ QCD , fracture functions in SIDIS can be factorized with collinear parton distributions and fragmentation functions. In this paper, we study the factorization up to twist-3. The matching of twist-2 fracture functions is straightforward, while the matching of twist-3 fracture functions is nontrivial in the sense that one has to keep the gauge invariance. Here we will mainly study the factorization or matching at the level of twist-3.
Fracture functions at twist-2 relevant to SIDIS have been defined in [6], where their asymptotic behavior has been derived for the region where the momentum fraction of the struck parton approaches its maximal value. In [15] TMD quark fracture functions have been classified for a polarized spin 1/2 hadron in the initial state. An one-loop matching of various twist-2 fracture functions for the process where the production of a lepton pair with large invariant mass in hadron collisions associated with a diffractively produced photon instead of a hadron in different kinematical regions, has been studied in [16].
In the target fragmentation region with Q ≫ k ⊥ SIDIS can be conveniently described with fracture functions. If one uses collinear factorization in the region of Q ≫ k ⊥ ≫ Λ QCD , then there will be terms with large log's like ln Q/k ⊥ in perturbative coefficient functions. Using the factorization with fracture functions and the matching studied here, these terms with large log's can be re-summed with helps of evolution equations. In collinear factorization, twist-3 effects involve in general twist-3 parton distributions and twist-3 parton fragmentation functions. It is interesting to notice that the latter is absent in the matching of twist-3 studied here. The reason for this is helicity conservation in QCD as we will show.
The content of our paper is: In Sect.2. we introduce our notations and definitions of relevant parton distributions, fragmentation functions and fracture functions. The results of matching of twist-2 fracture functions are also given in this section. In Sect.3. we derive the matching of the twist-3 fracture function which is relevant to single transverse-spin asymmetry. In Sect.4. the matching of the twist-3 fracture function responsible to a double-spin asymmetry is given. Sect.5 is our summary.

Definitions of Fracture Functions and SIDIS in Target Fragmentation Region
We consider an initial hadron h A with the momentum P µ = (P + , P − , 0, 0) and another hadron in the final state with the momentum k µ = (k + , k − , k 1 ⊥ , k 2 ⊥ ) with k + = ξP + . The transverse momentum k ⊥ is small in comparison with k + . The quark fracture functions or quark diffractive parton distributions can be defined with the density matrix: with ij stand for Dirac-and color indices. Here the transverse momentum of the struck quark is integrated over. To make the density matrix gauge-invariant, the gauge links along the n-direction are implemented: The decomposition of the density matrix has been done in [15]. We consider the case that the initial hadron is a spin-1/2 particle and the polarization of the final hadron is summed if it has nonzero spin. The spin vector of the initial hadron can be decomposed as: with s L as the helicity. Neglecting terms beyond twist-3, the density matrix can be decomposed as: where · · · denote those distributions defined with chirality-odd operators. These distributions will not contribute to SIDIS in target fragmentation region. The function F q and ∆F q are of twist-2, while F qT and ∆F qT are of twist-3. We neglect here the contributions beyond twist-3. P P k k We consider the SIDIS process: where the initial hadron is of spin-1/2 with the spin vector s. The initial electron can be polarized with the helicity λ e . We assume that the polarization of the hadron in the final state is not observed. At leading order of QED, there is an exchange of one virtual photon with the momentum q = k e −k ′ e between the electron and the initial hadron. The relevant hadronic tensor is: The standard variables for SIDIS are: The masses of hadrons and leptons are neglected. We take a frame in which the momenta of the initial hadron and the virtual photon are: i.e., the hadron moves in the z-direction, and the virtual photon moves in the −z-direction. We consider the case in which the transverse momentum k ⊥ is small with k 2 ⊥ ≪ Q 2 . In the region z h ∼ O(1), the produced hadron moves almost in the forward region of the virtual photon. This region is called as current fragmentation region. Transverse-Momentum-Dependent(TMD) factorization applies in this region. In the region of z h ≪ 1, the produced hadron moves almost in the forward region of the initial hadron. This region is called as target fragmentation region. In this region the hadronic tensor can be factorized with the above introduced fracture functions.
In the target fragmentation region, we write the momentum of the produced hadron as: We consider the experimental situation in which the initial hadron is polarized and the transverse part s µ ⊥ of the spin vector s µ is nonzero. The incoming-and outgoing lepton span the so-called lepton plane. In our frame specified in the above the azimuthal angle between the spin vector s ⊥ and the lepton plane is denoted φ s . Similarly, one defines the azimuthal angle φ h for the produced hadron. The azimuthal angle of the outgoing lepton around the lepton beam with respect to the spin vector is denoted ψ. In the kinematical region of SIDIS with large Q 2 , one has ψ ≈ φ s [17]. With this specification the differential cross-section is given by [17,18]: where α is the fine structure constant. At the leading order, the hadronic tensor receives the contribution from Fig.1 where the gray box represents the density matrix in Eq.(1). From Fig.1 the tensor is given by: where · · · are power-suppressed terms. With the fracture functions in Eq.(4) we can derive the hadronic tensor at the order. With the leptonic tensor at the leading order of QED, we have the differential cross-section in the target fragmentation region expressed with fracture functions as [15]: with We notice that the factorization in Eq. (12) holds in fact for the case with k ⊥ ≪ Q. In the case of Q ≫ k ⊥ ≫ Λ QCD , collinear factorization also holds. Therefore, fracture functions in this case can be matched to collinear parton distributions. The function F q and ∆F q are matched to twist-2 parton distributions and fragmentation functions. The function F qT and ∆F qT are matched to twist-3 parton distributions and twist-2 parton fragmentation functions. It is interesting to note that there is no contribution involving twist-3 parton fragmentation functions because of helicity conservation. This will be explained later in detail.
and fragmentation functions. The function qT and ∆ qT are matched to twist-3 parton distributions and twist-2 parton fragmentation functions. It is interesting to note that there is no contribution involving twist-3 parton fragmentation functions because of helicity conservation. This will explained later in detail.

Matching of and
The fracture function and ∆ can be matched to twist-2 parton distributions and twist-2 parton fragmentation functions, in the case of QCD . These twist-2 quark-and gluon distributions can be defined as: In the above, ) is the spin-averaged quark distribution function, ∆ is the longitudinally polarized quark distribution, is the quark transversity distribution. · · · stand for twist-3 or higher twists. At twist-2 there are two gluon distributions.
) is the spin-averaged gluon distribution function, ∆ is the longitudinally polarized gluon distribution. At twist-2 the relevant parton fragmentation functions of an unpolarized hadron are defined as:

Matching of F q and ∆F q
The fracture function F q and ∆F q can be matched to twist-2 parton distributions and twist-2 parton fragmentation functions, in the case of k ⊥ ≫ Λ QCD . These twist-2 quark-and gluon distributions can be defined as: In the above, q(x) is the spin-averaged quark distribution function, ∆q is the longitudinally polarized quark distribution, h 1 is the quark transversity distribution. · · · stand for twist-3 or higher twists. At twist-2 there are two gluon distributions. g(x) is the spin-averaged gluon distribution function, ∆g is the longitudinally polarized gluon distribution. At twist-2 the relevant parton fragmentation functions of an unpolarized hadron are defined as: where dq is the anti-quark fragmentation function, d g is the gluon fragmentation function. The observed hadron in the final state is with the momentum k µ = (k + , 0, 0, 0), where its mass is neglected.
With the twist-2 parton distributions and fragmentation functions, one can easily calculate the fracture functions F q and ∆F q at tree-level, where nonperturbative effects are represented by these distributions and functions. The contributions to F q and ∆F q are given by Fig.2. Calculations of these diagrams are rather standard. Therefore, we give our twist-2 results directly. More technical details will be explained in subsequent sections. From Fig.2 we obtain the fracture functions F q and ∆F q : with We notice that for k ⊥ ≫ Λ QCD , the twist-2 fracture functions behave like 1/k 2 ⊥ . The other two fracture functions F qT and ∆F qT are for the initial hadron which is transversely polarized. At first look, the twist-2 transversity distribution h 1 , which is chirality-odd, can give contributions to these functions from Fig.2, where twist-3 parton fragmentation functions are involved. At the order we consider, helicity or chirality is conserved in the amplitude of the upper part in Fig.1. Therefore, contributions involve chirality-odd parton distributions are absent. Hence, the two fracture functions only receive contributions where twist-3 parton distributions and twist-2 parton fragmentation functions are involved. Before we turn to the matching of these two functions, we give definitions of twist-3 parton distributions in the next subsection.

Definitions of Twist-3 Parton Distributions
From the quark density matrix one can define a set of three twist-3 parton distributions with one transverse derivative [19]: withs µ = ǫ µν ⊥ s ⊥ν . The three distributions are real. Another set of twist-3 distributions can be defined with a pair of quark fields with one gluon field strength operator. They are given by the matrix: where we have suppressed the gauge links for short notations. The matrix can be decomposed into Properties of these twist-3 parton distributions have been studied in [20,21,22,23,24]. The first two have been introduced in the study of single transverse-spin asymmetries. The last two are chirality-odd. They will not give contributions in the matching because of helicity conservation. The chirality-even parton distributions satisfy: Corresponding to the two distributions in Eq.(20) one can define additionally two twist-3 distributions by replacing the field strength tensor These two functions will not appear in our calculation. In fact they can be expressed with the distributions given in Eqs. (18,20) as shown in [25].
The three twist-3 distributions in Eq.(18) and the two twist-3 chirality-even distributions in Eq.(20) are not independent. One can show: where P stands for the principle-value prescription. The first relation has been derived in [26]. The second relation is derived in [19]. It should be emphasized that the second relation is for SIDIS, for q ′ ∂ defined with future pointing gauge links. The distribution q ′ ∂ (x) in Drell-Yan processes is defined with gauge links pointing to the past. With the symmetries of time-reversal and parity one can show that there is a sign-difference between the two distributions.
There are twist-3 gluon distributions defined only with gluon fields. One can first define three density matrices: where the field strength tensor G +ρ (λ 2 n) in the last line and the covariant derivative D ρ (λ 2 n) are in the adjoint representation. The three matrices are not independent. They are related through the relation: Besides the three matrices at twist-3 one can also find an additional twist-3 matrices conveniently in the light-cone gauge G + = 0. In this gauge, one can introduce where all indices α, β and γ are transverse. One can identify that: Therefore, there are two twist-3 matrices built with three gluon field strength tensor operators. From Bose-symmetry and covariance the two tensors take the form [27,28,29]: with the properties of the functions O and N The matrix M µνρ ∂ (x) can be parametrized with constraints of symmetries as: with The two functions g ∂ and g ′ ∂ are real. From the relation in Eq.(24) one has: Therefore, all twist-3 gluon distributions are determined by N (x 1 , x 2 ) and O(x 1 , x 2 ). In our notations, all twist-3 parton distributions have the dimension 1 in mass and are proportional to Λ QCD .

Matching of F qT
F qT describes the single transverse-spin asymmetry of SIDIS in target fragmentation region. It is well-known that such an asymmetry is a T-odd effect which requires nonzero absorptive parts in the scattering amplitude. In the case of k ⊥ ≫ Λ QCD the asymmetry is a twist-3 effect in the collinear factorization [20,21]. In this case, F qT can be matched with twist-3 parton distributions and twist-2 parton fragmentation functions. Because of helicity conservation, the contribution involving the twist-2 quark distribution h 1 with twist-3 parton fragmentation functions is absent as discussed before. At the order considered, there is no contribution from Fig.2 because of the absence of the required absorptive part. The contributions are given by Fig.3, Fig.4 and Fig.5, where the diagrams are given in terms of amplitude. The contributions are obtained from the interference of the amplitudes of Fig.3a, Fig.4 and Fig.5 with those of Fig.3b. There is one propagator in Fig.3a, Fig.4 and Fig.5 with a short bar. This implies that only the absorptive part of the propagator is taken into account. This gives the required absorptive parts in the amplitudes. The contributions are classified according to the parton momenta. The hard-pole contribution is from Fig.3 where the gluon parton carries nonzero momentum. The soft-gluon contribution is from Fig.4 where the gluon parton carries zero momentum. The soft-quark contribution is from Fig.5 where the quark parton is with zero momentum.
All contributions from these diagrams involve only the twist-2 gluon fragmentation function. The general structure of each type of the contributions from these diagrams can be written in the form with the quark-gluon correlator: d g is the gluon fragmentation function. H a,ρ L is the sum of upper parts of all diagrams in which the gluon line from the correlator in the left side of the cut, or in the left part of diagrams. The quark in the right part of the diagrams is with the momentum k A , the gluon is with k 1 and the quark in the left part is with the momentum k A − k 1 .
The quark-gluon density matrix can be written in the form: where · · · denote the terms beyond twist-3. The quark-gluon matrix elements are The power counting for the parton momenta and the quark-gluon matrix elements are: In the case of F qT , the matrix elements M ρµ ⊥ and M ρµ A⊥ contribute at twist-3 only with ρ = +, where one can neglect all transverse-and −-components of parton momenta. For ρ = +, the contributions from these two matrix elements are beyond twist-3.
We take the contributions from M ρ (k A , k 1 ) in Eq.(35) to discuss the collinear expansion and the gauge invariance of the results. The discussion in the case of M ρ A (k A , k 1 ) is similar. To obtain twist-3 contributions we need to expand the upper parts around the parton momenta: The upper part of the contributions involving M ρ (k A , k) is given by: The collinear expansion is: where terms represented by · · · give the contributions beyond twist-3. With the expansion, the contribution at the leading order of λ can be written in the form: where contributions beyond the leading power of λ are neglected. It is clear that the first term in Eq.(40) is gauge-invariant at the considered order of g s . The secondand third terms in Eq.(40) are not gauge-invariant. However, because of the additional cut, i.e., that the propagator with a short bar represents an on-shell particle, one can find the following Ward identity and the identity with the momentak 1,A : With the first identity, one can show that in the hard-pole-and soft-quark contributions the second-and third gauge-variant terms are zero. In the case of the soft-gluon pole contribution, because of that k + 1 = 0 one can not use the argument of the identity. But, we find by adding the contribution from complex conjugated diagrams that the final result is gauge invariant. The situation here is similar to the analysis of twist-3 contribution of SIDIS in [30]. For the contributions from M ρ A (k A , k 1 ) the results are similar except the soft-gluon-pole contribution which is nonzero and gauge-variant.
For the contributions from the matrix elements M ρµ ⊥ and M ρµ A⊥ with ρ = + there are similar identities to those in Eq.(41). In the hard-pole-and soft-quark contributions only the momentum component k + 1 and k + A are not zero. Therefore, there are no hard-pole-and soft-quark contributions from these two matrix elements because of an identity similar to the second one in Eq.(41). For soft-gluon-contributions, the argument can not be used because of that k + 1 = 0. With explicit calculations we find that the softgluon-pole contribution from M +µ ⊥ is zero if we add the contribution from complex conjugated diagrams. The soft-gluon-pole contribution from M +µ A⊥ and that from M ρ A are not zero. Both contributions can not be written in a gauge-invariant form. We notice that in the matrix element M +µ A⊥ one of the quark fields must be the −-component with the decomposition discussed in the next section. With equation of motion one can show that the sum of soft-gluon-pole contributions from M +µ A⊥ and M ρ A can be written in a gauge invariant form, i.e., in terms of T F . Therefore, the soft-gluon-pole contributions come from all matrix elements in Eq.(35) except M ρµ ⊥ with ρ = +. The calculation is tedious but straightforward. Hence, we will give our results of this section directly without giving the details of our calculations.
The hard-pole contribution is from Fig.3. The result is: The soft-gluon-pole contribution is from Fig.4. The result is There is no soft-gluon contribution from T ∆ because of that T ∆ (y, y) = 0. The soft-quark contribution is from Fig.5. The result is In these results y is given by y = x + ξ/z. Besides the twist-3 contributions from quark-gluon correlators, there are contributions from purely gluonic distributions at twist-3. Since an additional cut or an absorptive part is required, there is no contribution from Fig.2b. The contributions are only from Fig.6, where one quark propagator or eikonal propagator is with a short bar.
After completing the collinear expansion related to the produced hadron, the contribution from Fig.6 can be written in the form: This contribution is associated with the fragmentation function of an antiquark into the observed hadron. H abc Lµ 1 µ 2 µ 3 (k 1 , k 2 ) is the perturbative part which is the sum of the upper parts of Fig.6. In Fig.6a the left gluon line carries the momentum k 1 and the right one carries k 2 . In Fig.6b, the gluon line carries the momentum k 3 . Because of the cut and that the propagator with a short bar represents an on-shell particle, one can find the following identities: The contributions from the last two diagrams in Fig.6a are proportional to δ(k + 2 ) because of the cut on the eikonal propagator. With this fact we can find the identity: However, this identity is useless in the case of soft-gluon-pole contributions here. This is similar to the case discussed after Eq.(41). We find the leading contribution from Fig.6 and their complex conjugated diagrams can be written in a gauge-invariant form. Other contributions are beyond twist-3. We have the result with the twist-3 gluon distributions: with y = x + ξ/z. The final result for F qT is the sum: where terms in the sum can be found in Eqs.(42,43, 44,48).

Matching of ∆F qT
∆F qT represents the double spin asymmetry. Unlike F qT , the asymmetry is not zero in the absence of absorptive parts in the scattering amplitude. In the matching of ∆F qT , there are contributions from Fig.2. We first discuss the contribution from Fig.2a. It involves only gluon fragmentation function and can be written in the form: where k g is the momentum of the gluon which is fixed as zk g = k. p is the momentum of the right quark line. H kl (p) is the sum of the perturbative parts represented by the upper parts of the diagrams in Fig.2a.
With the power counting in Eq.(36) for parton momenta we can expand H kl (p) around the momentum p µ = (p + , 0, 0, 0): where · · · stand for contributions at higher orders. Taking the leading term and the twist-2 part of the quark density matrix in Eq. (14), one obtains the twist-2 contribution to F q and ∆F q given as the first terms in Eq. (16).
The twist-3 contribution is obtained by taking the second term in the expansion in Eq.(51) or the first term combined with the twist-3 part of the quark density matrix. The result is: with These contributions seems to be with the twist-3 distributions q T and q ∂ , respectively. But, they are not exactly those distributions. The contributions of gauge links are not included. When we consider the contributions in Fig.7 where an additional gluon exchanges, parts of them will be the contributions of gauge links at the considered order.
In calculating the contributions with one-gluon exchange given by Fig.7, we have the same quarkgluon correlator given in Eq.(33) and its decomposition in Eq.(35). By calculating the contribution with M ρ A , we find that a part of the contribution can be added to the second term in Eq.(52) so that the sum is obtained by inserting gauge links in the matrix element at the considered order of g s . With the gauge links the matrix element is that used to defined q ∂ . Excluding this part, the contribution from Fig.7 with In Eq.(54), the first two terms are gauge invariant, while the third-and fourth terms are not. The contribution from M ρ in the quark-gluon density matrix can be written: where the first two terms are gauge invariant. The third-and fourth terms are not gauge invariant. It is interesting to note that they have the same perturbative coefficient function F 0 as that in the gauge variant contribution in Eq.(54). Now we consider the contributions with M ρµ ⊥ and M ρµ A⊥ in the quark-gluon density matrix. The twist-3 contributions are given by taking the index ρ = +, i.e, only the field component G + is involved. We find that the contribution from M ρµ ⊥ is proportional to C A . In the contribution from M ρµ A⊥ , there is a part proportional to C F . This part can be summed with the first term in Eq.(52). The sum is obtained by inserting gauge links in the matrix element in the first term so that the matrix element is that used to defined q T . The sum is then gauge invariant at the considered order of g s . The sum of the remaining contribution from M ρµ A⊥ and that from M ρµ ⊥ can be written in the form: This sum is not gauge invariant and it involves the same perturbative function F 0 . We notice that the quark field can be written as the sum of a +-and −-component: The −-component is not independent. With the equation of motion one has: In Eq.(57) one quark field has to be the −-component. Using the solution in Eq.(59) one has the sum: Comparing the sum with the gauge variant contribution in Eq.(54, 56), we find that all gauge variant contributions are cancelled each other so that the remaining contribution is gauge invariant. We have then the result which is gauge invariant: with H 2p,∂ and H 2p,T are given in Eq.(52). H(x, ξ, x 2 ) and H A (x, ξ, x 2 ) are given as: ∆F qT receives contributions from twist-3 gluon distributions, where the antiquark fragmentation function is involved. These contributions are given by diagrams of two-gluon exchanges as given in Fig.2b and by three-gluon exchanges in Fig.6 without the cut on the quark propagator and the gauge link. It is noted that the contribution of three-gluon exchanges from the second-and third diagram in Fig.6a is the same. Because of Bose-symmetry, one should either take any one of the two diagrams, or the half of the sum into account.
The contribution of two-gluon exchanges from Fig.2b can be written in the form: where H µν (k 1 ) is the sum of the perturbative parts represented by the upper parts of the diagrams in Fig.2b. k 1 is the momentum carried by the gluons. It is easy to check the following Ward identities: Before doing the collinear expansion in k 1 aroundk µ 1 = (k + 1 , 0, 0, 0), we can use these identities to manipulate the expression: In the first term, the leading contribution in the collinear expansion of H µ ⊥ α ⊥ (k 1 ) gives the twist-2 contribution, i.e., to ∆F q . The twist-3 contribution is obtained by expanding the next-to-leading contribution. The leading contribution of the second-and third term is at twist-3. The last term is a twist-4 contribution which can be neglected. In the collinear expansion of H µν , we find that the next-to-leading contribution H µ ⊥ α ⊥ (k 1 ) is antisymmetric in the indices µ ⊥ and α ⊥ . Because of symmetries, the related matrix element is symmetric in the indices µ ⊥ and α ⊥ , as shown in Eq. (29). Hence, the first term gives no contribution. Therefore, only the second-and third gives nonzero contribution at twist-3. The contribution can be expressed by the twist-3 gluon distribution: g 3T (x)s µ = i x dλ 2π e −ixλP + h A |G +− (λn)L † n (λn)L n (0)G +µ (0)|h A .
We obtain the contribution from Fig.2b: It is noted that from Fig.2b we only obtain the contribution with g 3T defined only withĜ µν and without gauge links. The contribution alone is not gauge invariant. There is a difference of a two-gluon term betweenĜ µν and G µν . In the calculation of three-gluon exchanges given by diagrams in Fig.6, we find that a part of the contribution from Fig.6 gives the needed two-gluon term in G µν and another part forms the gauge links in g 3T at the considered order of g s . The contribution from three-gluon exchanges is given by diagrams in Fig.6. It can be written where H abc µ 1 µ 2 µ 3 (k 1 , k 2 , k 3 ) is the sum of the perturbative parts represented by the upper parts of the diagrams in Fig.6. The two gluons in Fig.6a carry the momentum k 1 and k 2 , respectively. k 3 is the momentum carried by the gluon in Fig.6b. To obtain the gauge invariant result, we first note that there The result of ∆F qT from Fig.2b and Fig.6 is: −N (y − x 2 , y) + (2ξy + yx 2 z − y 2 z − ξx 2 ) N (y, x 2 ) + O(y, x 2 ) The complete matching result is the sum of that in Eq.(61) and that in Eq.(75): 5. Summary SIDIS in target fragmentation region can be conveniently described with fracture functions, i.e., one can use QCD factorization with fracture functions to make predictions in this region. If the transverse momentum of the produced hadron is in the region Q ≫ k ⊥ ≫ Λ QCD , the standard collinear factorization can also be used. Therefore, fracture functions can be factorized or matched with parton distributions and fragmentation functions. We have studied the matching up to twist-3 level. At the order of α s considered in this work, fracture functions are only matched to twist-2 parton fragmentation functions with twist-2or twist-3 parton distributions. There is no contribution from chirality-odd parton distributions. We have derived perturbative coefficient functions in the factorization or matching of twist-2-and twist-3 fracture functions. Especially in the derivation of twist-3 fracture functions, we find that the results can be written in a gauge invariant form. These results will be useful for modeling of fracture functions and resummation of large logarithms of k ⊥ in collinear factorization.