1-jettiness with jet axis at $O(\alpha_s)$ in jeep inelastic scattering

We present $O(\alpha_s)$ analytic predictions for event shape 1-jettiness $\tau_1$ distribution aiming measurements in deep inelastic scattering process at future Electron Ion Colliders. The result depends on conventional variables $x$ and $Q$ as well as on $\tau_1$ and is relatively compact and easy to implement for numerical calculation. Three different choices of axis, with respect to which $\tau_1$ is measured are considered in the Breit frame. The first is the one optimally adjusted to minimize $\tau_1$ and the second and third are taken from anti-$k_T$ and Centauro jet algorithms defined with a jet radius parameter $R$, respectively. We find that the first and second give the same result at this order and are independent of $R$, while the third depends on the radius. This fixed-order result provides a nonsingular contribution to be combined with a singular log-resummed contribution to give the full spectrum in $\tau_1$ space and also shows how fixed-order and resummation regions change as a function of $x$ and $Q$.


Introduction
Jets, energetic hadron bunches produced in high-energy collisions, are useful tool to study the strong dynamics induced by Quantum Chromodynamics (QCD) and to probe new physics beyond the Standard Model. Its application to future QCD machines that probe internal structure of ions with electron beam is extensively investigated in Electron-Ion-Collider in the US (EIC) [1,2] and in Electron-ion-collider in China (EicC) [3]. The jets can be defined and studied exclusively using jet finding algorithms [4][5][6][7][8][9], which conventionally require parameters like a jet radius and jet veto. On the other hand, they can also be studied inclusively with classic observables called event shapes [10], which are theoretically simpler with a small number of parameters, hence easier to achieve higher accuracy compared to exclusive jet study. An example of event shapes is a thrust which has been predicted up to N 3 LL+O(α 3 s ) accuracy [11][12][13][14][15][16] in e + e − annihilation. From the event shape, the strong coupling constant was determined at 1% precision [15][16][17], which is one among the most precise determinations listed in Particle Data Book [18]. The thrust as well as other event shapes in deep-inelastic scattering (DIS) was studied [10] and measured in HERA experiment [19][20][21][22][23][24]. Due to limited detector coverage they were defined from particles in a current hemisphere, to which products of hard scattering usually belong while initial-state radiations and beam remnants do not. On the other hand, the definition in e + e − is done with both hemispheres. More recent developments including improved accuracy and/or new event shape predictions [25][26][27][28][29][30][31] assume that future machines can cover both regions. With the assumption, DIS thrust, which we call a 1-jettiness [25,26] can be defined as 1 where p i is the momentum of particle i in the final state X and Q is a hard momentum transferred by a virtual photon. q B and q J are beam and jet axes, onto which momentum p i is projected. The operator 'min' takes smaller one among two scalar products and this makes particles grouped into one of two regions, beam or jet. The value of τ 1 is small when the final state X contains two collimated bunches along each of beam and jet axes. Otherwise, for multi-jet final state the value is not small. In [25], a version called τ b 1 was defined by z axis in the Breit frame for q J and known to be same as a version of DIS thrust called τ Q measured in HERA [33]. It was computed analytically at the first order in α s [34].
In this paper, we study another version called τ a 1 [25,35] for which the jet axis q J is determined by axis finding algorithms, while q B defined next section is held fixed to the beam axis. Its distribution was numerically computed in the Lab frame [26] by using the axis obtained from anti-k T algorithms. In this paper we consider three different algorithms in the Breit frame in Sec. 2 and show analytic expressions of our fixed-order predictions at the first order in α s and separately next-to-leading power (NLP) terms in small τ 1 limit in Sec. 3 as well as numerical results in Sec. 4. Finally, we summarize in Sec. 5.

1-jettiness with three jet axes
In this section we define the jet axis q J and write down the expression of 1-jettiness defined with the axis for 2-body final state. We consider three different jet axes in the Breit frame. The first axis is the axis that minimizes the value of 1-jettiness as in e + e − thrust and we call it a 1-jettiness axis. Next two axes are obtained from jet momentum defined by anti-k T algorithm [9] and by Centauro algorithm [36].
We first make a quick review of kinematics in the Breit frame. In the frame, the virtual photon with momentum q and the proton with momentum P are aligned along the z axis. One of advantages of this frame is that the initial state radiation moving along the proton direction is well separated from the other particles produced by a hard scattering. We take 2 . (2.4) The scalar products of the beam axis q B and p i contributing to 1-jettiness can be expressed as which is a constant in v. We used the hadronic variable x in advance instead of the partonic variable z, which is replaced by x at the end. The maximum occurs on the heavy lines in Fig. 1.

Anti-k T algorithm τ kt
In the k T -type algorithms [9] designed to be invariant under longitudinal boosts, the distance between particles is defined in term of transverse momentum p T and invariant angular distances where d ij is the distance between outgoing particles i and j, and d iB is the distance of a particle i from the beam. p T i is the magnitude of the transverse momentum of particle Figure 1. Two-body phase space for 1-jettiness, τ jt,kt with jettiness/anti-k T axis (left) and τ ct with Centauro axis (right). 1-jettiness takes the same expression in first three regions while the fourth region is only for τ ct . The bound of fourth region for jet radius R = 1 and 2 is shown in red dashed and red solid, respectively.
i, R is a jet radius parameter of O(1) and ∆R ij is an angular distance defined in terms of distances in rapidity y and in azimuth ϕ. The power p on the momentum is a real valued parameter. When p = 1, 0, −1, it is called k T , Cambridge/Aachen and anti-k T , respectively. In this paper we take the anti-k T algorithm p = −1, which is a conventional choice by many experiments. With above distances, the jet algorithm finds minimum d min among all the d ij and d iB . If d min is a d ij , particle i and j are merged into a single particle and if d min is a d iB , the particle is declared to be a final jet and removed from the list of final particles. We repeat this with new list of particles until all the particles in the list are gone.
For a two-body final state, we need to compare d 12 , d 1B and d 2B to classify jets. Recall that in the Breit frame the transverse momenta of two particles in Eq. (2.3) are back-toback so that The distance in azimuthal angle between p 1,2 is (ϕ 1 −ϕ 2 ) 2 = π 2 and the distance in rapidity is (y 1 − y 2 ) 2 = ln 2 1−v v . Then, the distance d ij is Typical size of jet radius R is smaller than π and ln 2 1−v v > 0 hence, d ij is always greater than d 1B . This means that there is no merging in the two-body final state. The anti-k T algorithm gives two jets and each of them is p 1 and p 2 and the jet momentum of each jet is p i . The jet axis q J is defined from the jet momentum as where in the second equality, we used the fact that p i is massless. In higher multiplicity, the jet algorithm returns more jets. Among those jets, what we need is collinear jets with a large momentum that are candidates for the jet axis q J and by implementing a veto condition we can reject unnecessary soft jets with a small momentum. Among the candidates, the one that gives the smallest value of 1-jettiness will be selected. Then, 1-jettiness with the candidates p 1 and p 2 is given by where in the second line we set q J = p 1 and q J = p 2 after first and second summation symbols, respectively. Then, as was done in Eq. (2.7), all possible combinations are taken into account in the third line. To get the fifth line, we used Eq. (2.5) then, deleted trivial terms that cannot be taken as a minimum. Interestingly, we come up with a fact that τ kt is equal to τ jt in the two-body case. This is because of simplification with a small number of final particles. In general, the radius R dependence enters in the case of anti-k T algorithm and we do not expect that τ jt and τ kt are the same with higher multiplicity. We do not consider higher multiplicity in this paper and from now on, we do not separately treat τ kt except for the case when we need to distinguish τ kt from τ jt .

Centauro algorithm τ ct
For the third jet axis we would like to adopt the Centauro algorithm [36], which is more recently introduced algorithm for DIS taking into account the target-current asymmetry in the Breit frame. It is still longitudinal-boost invariant like the anti-k T algorithm but it allows to capture jets close to beam axis while the anti-k T algorithm cannot form a jet in the region since the rapidity distances between particles across the beam line become very large (y i − y j ) → ±∞. The Centauro algorithm defines the following distance measure where ∆f ij and ∆φ ij are differences between f i and f j and φ i and φ j , respectively. The function f i is a function of an angular variableη i of particle i whereη diverges in forward region asn z ·p i → 0 and thus prevents jets from enclosing the proton beam direction, while it decreases as particles become closer in backward region For the function f , we take the simplest form in [36] f (η) =η . (2.18) Then, f 1 and f 2 and the distance d 12 are given by , (2.20) where R is a radius parameter of order O(1). If d min is a d iB , we have the same jet as the anti-k T jet. If d min is a d 12 , we have a jet with momentum p 1 + p 2 , which is not allowed in the anti-k T algorithm. The condition for two particles merging into a single jet is d 12 < d iB , which indicates the following region where v ± are upper and lower bounds of the variable v and v ± = 1/2 at the value (2. 22) In this region, the jet momentum is p 1 + p 2 and q J always is in the opposite to the proton direction wheren 12 = (p 1 +p 2 )/|p 1 +p 2 |. Otherwise q ct J = q kt J is one of p 1 and p 2 as in Eq. (2.14). The right panel in Fig. 1 indicates the phase space for 1-jettiness with the Centauro algorithm. The bounds of region IV for R = 1, 2 are shown in the figure in red curves. The region reduces to zero as R → 0 and invades the regions I and II for R > 2. In many experiments, R is taken to be less than or, equal to 2 and the phase space in this range is also relatively simple. So, we constrain the region of our interest to be R ≤ 2 . (2.24) Then, in the region IV the value of τ ct can be computed with Centauro-jet axis in Eq. (2.23) and in the regions I, II, and III τ ct = τ jt . Now τ ct in any region can be expressed as Note that τ ct is always positive in the region IV as well as the other regions because the minimum value of z is z c (R) ≥ 4/5 in the range given by Eq. (2.24).

Analytic cross section at order α s
In this section we summarize the analytic result for 1-jettiness cross section in DIS computed at the first order in α s . Details of the calculations are given in App. A. Here, our results are given in the form of cumulative cross section obtained by integrating the differential cross section The DIS cross section is conventionally decomposed into the structure functions F i This can be expressed in terms of F 1 , F 2 by using the relation between structure functions F 2 = F L + 2xF 1 . Note that the structure functions are the functions of τ as well as x and Q. At the point τ = τ max , the cumulative functions σ c and F i reduce to the inclusive cross section and inclusive structure functions, respectively. The F i are coefficients of basis tensors that decompose a current-current correlator called the hadronic tensor W µν defined in Eq. (A.3) and one can read off the individual coefficients by projecting with orthogonal tensors like P µ P ν W µν and g µν W µν as shown in Eq. (A.6).
The structure functions are written in terms of these projected correlators where A i and B i are respectively equivalent to P µ P ν W µν and g µν W µν up to multiplicative factors that can be found in Eqs. (A.54) and (A.57). They contain logarithmic terms (singular) that were obtained by the factorization formula in [25,35]. Non-logarithmic terms (nonsingular) are obtained in this study by fixed-order QCD calculations. Each of A, B can be written like A = A sing + A ns and B = B sing + B ns . For the completeness we copy and paste the singular parts where P qq (z) and P qg (z) are the splitting functions and L n (1 − z) is a plus distribution Note that the anti-quark contributions Aq and Bq are the same as A q and B q except for the quark PDF replaced by the anti-quark PDF.
The nonsingular parts are computed in App. A. The final expressions are given by where Θ 0 represents physical region of τ for a given value of x In comparison to τ b , the singular and nonsingular parts of τ jt have many terms in common with those of τ b in [34]. Their differences are summarized in App. C. As shown in Eq. (2.25), τ ct is different from τ jt in the region IV. We take their differences in the region and denote them by δA i and δB i . Then, the structure functions for τ ct are obtained by replacing A by A + δA and B by B + δB in Eq. (3.3). Their final expressions are given by and the parameters z jt,ct and r(z, R) are given by The differential distributions can be obtained by differentiating Eqs. In addition, by expanding A ns q , A ns g , B ns g and B ns q in the τ → 0 limit we can obtain the NLP term which is the power correction to the singular terms in Eq. (3.4). At the leading power the corrections contain terms like τ and τ ln τ , which is suppressed by τ compared to the singular terms. The NLP obtained from Eq. (3.6) can be expressed as At O(τ ) the power corrections are zero except for δB ns q . So, expanding up to first nonzero correction we have It is worth to note that PDFs derivatives appear in the NLP results which should be identified in a factorization at the subleading power and O( √ τ 1 ) corrections are absent for any of the jet algorithms presented since power corrections of 1-jettiness scales like O(τ n 1 ) with integer n. Another feature to note is that the appearance of a jet radius dependence in Eq. (3.14) is expected from the effective theory analysis [37] that implies that power corrections are sensitive to the jet algorithms and the logarithmic structure due to the phase-space cutoff typically appear at the leading power in jet substructure observables [38] and in transverse vetos [39].
In the context of high-order calculations in perturbative QCD, the observable Njettiness is used to control and to subtract the infrared singularities in numerical computation [40][41][42], the results in Eqs. (3.13) and (3.14) are essential ingredients required to improve the subtraction accuracy to subleading order [43,44], while the singular part in Eq. (3.4) is for subtraction at the leading order. The results also serve as an important crosscheck to be reproduced by effective field theory approach using the factorization at subleading power on the way to the higher-order in α s . It is worth to point out an interesting observation in our results. The algorithms, 1-jettiness and anti-k t , are the same at the order α s hence, their subtractions to arbitrary accuracy are identical at this order. On the hand, the jet algorithms, anti-k t and Centauro, are different in the subtraction at the subleading power because their difference in Eq. (3.14) contains non-vanishing τ ln τ terms.

Numerical results
In this section, we show our numerical results obtained by using analytic 1-jettiness cross sections at order α s in Sec. 3. We compare singular and nonsingular parts at different values of x and Q. Numerical difference between τ jt and τ ct are also presented for several values of R. We also compare our result to τ b result given in [34]. The crossing point between singular and nonsingular parts implies a boundary between fixed-order and resummation regions and the point as a function of x at several values of Q is presented. In choosing the values for the kinematic parameters x and Q in the section, we considered the region that can be studied in future experiments such as EIC or EicC.
Main achievement of this paper is the analytic expressions of nonsingular part of τ jt and τ ct . In Fig. 2   About the other variant τ ct , its singular part is the same as that of τ jt and the nonsingular is similar to τ jt . Instead of showing the similar style of plots for τ ct , we compare difference between τ jt and τ ct in Fig. 3 by using the expression in Eq. (3.8). The size of difference increases with R increasing and it becomes as large as 25% when R = 1.8. With smaller values of R the difference quickly reduces to zero and this is consistent with the R dependence of the region IV in Fig. 1 left panel, where the region shrinks to zero as R → 0. A feature of sharp turning upwards in the plot is associated with the contributions from τ ct vanishing in τ ≥ (16 − R 2 )R 2 /256. Then, when the other contribution from τ jt vanishes in τ ≥ R 2 /16 the difference becomes zero as shown in the figure.
In Fig. 4 we compare our result to another version of 1-jettiness called τ b in [34]. Because of difference in factorized formula between them, their singular parts are different [25]. However the difference in singular part is proportional to δ(τ ) and is not visible in a differential distribution like Fig. 4. Therefore, the difference from τ b shown in the plot is from nonsingular parts and the size of difference is larger at lower value of Q. Note that as shown in the plot τ b does not vanish beyond τ = 1/2 because the maximum of τ b is 1, while it is 1/2 for τ jt and τ ct . For the uncertainty figure, the uncertainties are obtained by varying the renormalization scale from µ = Q by a factor of two up and down (dσ(Q) − dσ(2Q))/dσ(Q) and (dσ(Q) − dσ(Q/2))/dσ(Q) as upper and lower boundary, respectively. We can find that τ jt and τ ct share almost same uncertainty. They only have slightly difference with each other within τ = 1/2. This uncertainty increases with the decrease of the values of Q. We also compare uncertainties for τ ct with different R value and find that the effect is around 0.01% changes.
The crossing point where singular and nonsingular parts meet each other can be understood as a boundary between fixed-order and resummation regions. In the fixed-order region, the singular part is smaller than the nonsingular part hence, an ordinary fixedorder QCD result is valid, while in the resummation region, the singular part is larger than nonsingular part due to increasing logarithms and resummation of those logarithms is necessary. In Fig. 5 the crossing points in τ as a function of x are shown and the results for τ jt (solid) and τ b (dashed) are similar as implied from Fig. 4.
The color density plot in Fig. 5 represents a relative size of the nonsingluar part to the full differential cross section at Q = 15 GeV for τ jt . An absolute values are taken for simplicity. In the blue left and lower corner the nonsingular is small while the singular dominates the cross section. In the light-colored region the singular and nonsingular are comparable to each other. Finally, in deep-red right and upper corner unphysical singular and nonsingular are largely cancelled to give the full cross section. In the region their absolute magnitudes easily become greater than full cross section and the deep-red region should be understood as 100% or, greater relative to the full cross section.
An important feature in the plot is that the resummation region increases with decreasing value of x and the region gets close to the maximum of τ jt near x = 0.01 while it does near x = 10 −5 for τ b . These crossing points imply when the resummation should be turnoff in τ spectrum and specifically in the scale profile function [25,34] the points can be taken to be the value of a parameter t 2 which is the point where the resummation begins to be turned off. In order to understand the behavior shown in Fig. 5, we take a following form where the left side is from the singular part and the coefficients α(x) and β(x) are obtained by differentiating Eq. (3.4). where δ(x) = (c ns /α(x)) exp[−β(x)/α(x)]. We can find the values of c ns that makes the approximate solution fitted to the curves in Fig. 5 obtained using the known nonsingular part. We find that c ns between 0 and −1 times β(x), represented by the grey dotted lines in Fig. 5, shows a reasonable description to the curves. Therefore, the singular part is mainly responsible to the behavior of crossing point as a function of x and in the absence of nonsingular part, the approximation in Eq. (4.4) can be used as a first time estimate for the crossing points and to fix a corresponding parameter of scale profile function.
Our numerical results given in this section are purely perturbative results. Including nonperturbative corrections and hadronization effects are important in precision predictions. The hadronic effects are power suppressed by Λ QCD /(τ Q) for τ Λ QCD /Q in small τ region and the correction at the leading power can be parameterized using a nonperturbative parameter Ω 1 [16,45,46]. The parameter is universal for any version of 1-jettiness and the dependence on the jet algorithm would remain small in this region. One also can take a shape function method that takes into account nonperturbative behavior as well as hadronic power correction [25,47]. However more recent studies imply that more careful considerations are required in a scheme related to renormalon subtraction [48] and in hadronization effect away from dijet region [49]. Furthermore, recent analysis using Monte Carlo simulations implies its fine-tuning is required to explain HERA measurements [50]. Therefore, a quantitative analysis on these effect in 1-jettiness would be beyond the scope of this work and could be done in future project.

Conclusions
In this paper we study the event shape 1-jettiness in DIS at the first order in α s that can be measured in future experiments. We considered three different jet axes, onto which a particle momentum is projected to compute the value of 1-jettiness. τ jt is a version with the jettiness axis that is optimally adjusted to minimize the value of 1-jettiness. The other two versions are τ kt and τ ct that take their axes from exclusive jet algorithms such as anti-k T and the Centauro algorithms in the Breit frame, respectively. We find that τ jt and τ kt are equivalent for the two-body final state, i.e., at the order α s .
Our main result is the predictions for τ jt and τ ct distributions at the first order in α s analytically expressed in Eqs. (3.6) and (3.8). They are expressed in the form of the cumulative distribution and in the Appendix the differential distribution is also given. The results are expressed such that one can easily write the structure functions such as F 1 (x, Q 2 , τ ) and F L (x, Q 2 , τ ) as well as the cross section. The results of τ ct share many terms in common with those of τ jt and their difference depending on R is given so that τ ct is obtained by adding the difference on the top of τ jt result. Comparison to the analytic result for τ b using Breit frame axis are also given in App. C.
Numerical results of our predictions are presented at different values of x and Q. Singular and nonsingular parts are compared in τ jt distributions. For τ ct , we studied the difference from τ jt and found that it is sensitive to the value of radius R and the magnitude is in the range of 5 ∼ 25% of τ jt distribution when the value of R is from 1 to 1.8. We also studied the singular-nonsingular crossing point τ cross as a function of x and Q. The crossing implies a boundary between resummation and fixed-order regions and the value can be used as a reference point for turning off the resummation. The value of τ cross hence, the resummation region increases with the decreasing value of x, and the region gets close to τ max = 1/2 at x = 0.01, while the value is less sensitive to Q. We found that the x dependence is well explained by the singular part and obtained an approximate expression for τ cross , which can be used in the absence of nonsingular part.
Our results provide an important piece of information toward precision predictions of event shapes in DIS. Our prediction for nonsingular part combined with resummed singular part can be measured in the future EIC and EicC and can be used to determine the strong coupling constant and a universal hadronic nonpertubative parameter.

APPENDIX A Details of calculations
In this section we calculate 1-jettiness QCD cross section. In derivations of the cross section, we follow the steps done in [25]. The cross section is expressed in the product of lepton and hadronic tensors.
where L II µν is the lepton tensor given in [25] and the index I = V, A implies vector and axial currents. Although, we consider II = V V in this paper but we keep these index to maintain generality for a moment. For vector currents, we have where −g T µν = −g µν + 2/Q 2 (k µ k ν + k µ k ν ), with k µ and k µ being the momenta of incoming and outgoing leptons and the photon momentum is given by q µ = k µ − k µ . One can show that the lepton tensor is transverse to virtual photon as implied by the Ward identity q µ L V V µν ∝ q µ g T µν = 0. The hadronic tensor W µν that measures the 1-jettiness from final hadronic state X can be expressed as where there is an average over incoming spins implicit. In the second line we specify the number of final particles and the corresponding phase-space integral. The hadronic tensor is decomposed into several tensors where F i are called structure functions, which is differential in τ here and integrating F i over τ gives ordinary structure function. The tensor T µν i are given by Multiplying g µν and P µ P ν by Eq. (A.3), one obtains linear combinations of F 1 and F 2 . Solving for F 1 and F 2 , they are expressed in terms of the hadronic tensor. Similarly, F 3 can be obtained by multiplying T 3 µν by Eq. (A.3).
In this paper we consider II = V V , there is no F 3 contribution. From hereafter we drop the index for the hadronic tensor W µν = W µν V V and one needs to calculate g µν W µν , P µ P ν W µν . The cross section in Eq. (A.1) in terms of the structure functions F i is given by Eq. (3.2).

A.1 Phase space integral
The phase-space integral can be explicitly written as where d = 4 − 2 . For n = 1, For n = 2, Fig. 1 shows three and four regions in v-z space divided for τ jt,kt and for τ ct , respectively. Then the phase-space integral in Eq. (A.9) splits into three pieces as where Θ 0 and v ± are given in Eqs. (3.7) and (2.22), respectively. Note that the region III in Eq. (A.10c) is that of τ jt shown in Fig. 1. In τ ct computation we use the same expression, which incorrectly include the region IV then, in Eq. (A.10d) the incorrect contribution from τ jt is subtracted. It makes Eq. (A.10d) finite so that we can set to zero and also makes the expression of τ ct cross section simple with an additional term added to τ jt cross section.

A.2 Hadronic tensor for incoming quark
The hadronic tensor defined in Eq. (A.3), the matrix element with incoming proton can be factorized into convolutions of proton PDF and Wilson coefficients. The Wilson coefficients can be determined by matching the factorization formula with incoming parton to the hadronic tensor for the parton state. The hadronic tensors W q µν and W g µν for initial quark and gluon, respectively should be calculated. To O(α s ) W q µν involves one-loop calculation and W g µν is simply tree level calculation. In this section, we calculate W q µν . W q µν receives contributions from tree, virtual, and real diagrams so can be written as Here, we suppressed the superscript q on right side and we will suppress it in the middle of calculations.

A.2.1 Tree-level and virtual contributions
The tree-level amplitude is M (0) µ = Q fū (p 1 )γ µ u(P ). The tree-level hadronic tensor is given by where we performed the spin average explicitly over quark spins σ then used the 1-body phase space in Eq. (A.8). The transverse metric is g T µν = g µν − (n zµnzν + n zνnzµ )/2. This gives the projected hadronic tensors where the second line vanishes by the Dirac equation of massless particle P / u(P ) = 0. The virtual contribution can be found from many literature for instance (14.19) [51] and the contribution is the same as that of τ b in [34] −g µν W vir where L = ln µ 2 Q 2 . The factor (1 − ) is not expanded because it is to cancel 1/(1 − ) in Eq. (A.6). In the second line of Eq. (A.14) we have used the MS scheme by re-scaling the scale µ 2 by e γ E /(4π) such that we use following replacement Note that the finite part in Eq. (A.14) is the α s term in the hard function in [25], derived in SCET in [52,53]. The P µ P ν W vir µν is zero because P / u(P ) = 0.

A.2.2 Real contribution
The definition of τ in Eqs. (2.7) and (2.25) divide the phase space into three or four regions as shown in Fig. 1 hence the 2-body phase-space is accordingly divided as follows where in the last line the factor 1/2 takes account the average over quark spins σ. In this subsection, we discuss the results of regions I,II, and III, which is relevant to τ jt and for the region IV, which is IR finite and simpler, we just give the final result in App. B.
The real emission contribution to the two projections of the color-averaged squared amplitudes that we need are given by where the first line can be found from (14.21) in [51].
µν , there is no singular term and we can safely set = 0. Including the spin average, we obtain where the Θ 0 (τ, z) is given in Eq. (3.7). The sum of Eq. (A.20) is given by where the superscript q representing the incoming quark is made explicit, again. Note that the tree-level and virtual contributions are zero and the real contribution is the total.
The regions I and II share the same Θ 0 (τ, z), so we combine them and expand in where Θ 0 (τ, z) was divided into two parts, and the first part multiplied by terms 1/τ 1+ and 1/(1−z) 1+ gives 1/ pole, while with the other part the only term like 1/τ 1+ gives the pole. We now expand above equation in powers of . For the term like θ(−τ + z 1−z )/τ 1+ /(1−z) 1+ we use a plus distribution identity in Appendix A of [34].
where P qq is the splitting function in Eq. (3.5). The functions I s fin and I ns fin contribute to singular and nonsingular parts, respectively.
where I ns is identical to corresponding part in τ b and is given in Eq. (A.8) of [34]. We use following relations to compress Eq. (A.24) where standard plus distributions L n (z) are defined by [34] L n (z) ≡ lim ε→0 d dz With a test function g(z) well behaving near z = 0, the integration against the function gives A variation of L n with a variable lower bound z 0 and its integral is given by For the contribution from the region III, it contains singular terms when τ or, 1 − z approaches zero. We isolate the singular part carefully by subtracting and adding singular terms where we used the delta function in the first line and replace z by τ . Finally we expand Eq. (A.30a) in It is a non-trivial cross check to show that sum of −g µν W (i) µν (τ ) and P µ P ν W where L = ln µ 2 Q 2 . First, the two logarithmic terms above are cancelled by the same terms in virtual part Eq. (A.14). The 1/ term above will be replaced by 1-loop correction of the proton PDF during matching procedure. The logarithmic scale dependence proportional to P qq in the last term should cancel the same scale dependence from RG evolution of the proton PDF. Therefore, all µ dependence at O(α s ) are cancelled.

A.5 Convolution with PDF
The hadronic tensor with a hadron state h is expressed in the factorized form as where f i/h is the PDF for the initial hadron h into a parton i and the superscript i on the coefficient w i represents the contribution from the parton i. Note that we obtained the tensors W q,g µν for incoming quark and gluon and by the perturbative matching performed in [34] one can translate them into the coefficients w i . We do not describe the matching procedure, and the expressions for w i are essentially the same as W q,g µν . Instead of writing w i , we give expressions for the hadronic tensor Eq. (A.53) convolved with the PDF. The (C.2)