TMD factorization for di-jet and heavy meson pair in DIS

We study a transverse momentum dependent (TMD) factorization framework for the processes of di-jet and heavy meson pair production in deep-inelastic-scattering in an electron-proton collider, considering the measurement of the transverse momentum imbalance of the two hard probes in the Breit frame. For the factorization theorem we employ soft-collinear and boosted-heavy-quark effective field theories. The factorized cross-section for both processes is sensitive to gluon unpolarized and linearly polarized TMD distributions and requires the introduction of a new soft function. We calculate the new soft function here at one loop, regulating rapidity divergences with the $\delta$-regulator. In addition, using a factorization consistency relation and a universality argument regarding the heavy-quark jet function, we obtain the anomalous dimension of the new soft function at two loops.


Introduction
It is well known that gluons are an essential constituent of nuclei and that the gluon parton distribution functions (PDFs) can numerically be much bigger than the corresponding quark distributions, especially when the parton energy fraction is small.Gluon transverse momentum dependent distributions (TMDs) are also expected to be similarly enhanced, however they result to be difficult to access due to the lack of clean processes where the factorization of the cross-section holds and incoming gluons constitute the dominant effect.An example of such a process is the Higgs production in hadronic dijet LO process: heavy meson pair at LO: < l a t e x i t s h a 1 _ b a s e 6 4 = " J 8 c o P N 5 h k T 5 1 P A 5 E v k X I y t Y y S M A = " > A A A B 8 3 i c b V B N S w M x E M 3 W r 1 q / q h 6 9 B I t S P Z T d K u i x 4 M V j B f s B 3 b X M p t l t a J J d k q x Q S v + G F w + K e P X P e P P f m L Z 7 0 O q D g c d 7 M 8 z M C 1 P O t H H d L 6 e w s r q 2 v l H c L G 1 t 7 + z u l f c P 2 j r J F K E t k v B E d U P Q l D N J W 4 Y Z T r u p o i B C T j v h 6 G b m d x 6 p 0 i y R 9 2 a c 0 k B A L F n E C B g r + V U / B i H g 4 R z H Z / 1 y x a 2 5 c + C / x M t J B e V o 9 s u f / i A h m a D S E A 5 a 9 z w 3 N c E E l G G E 0 2 n J z z R N g Y w g p j 1 L J Q i q g 8 n 8 5 i k + s c o A R 4 m y J Q 2 e q z 8 n J i C 0 H o v Q d g o w Q 7 3 s z c T / v F 5 m o u t g w m S a G S r J Y l G U c W w S P A s A D 5 i i x P C x J U A U s 7 d i M g Q F x N i Y S j Y E b / n l v 6 R d r 3 k X t f r d Z a V x m s d R R E f o G F W R h 6 5 Q A 9 2 i J m o h g l L 0 h F 7 Q q 5 M 5 z 8 6 b 8 7 5 o L T j 5 z C H 6 B e f j G 4 N c k J k = < / l a t e x i t > ( ⇤ f ) < l a t e x i t s h a 1 _ b a s e 6 4 = " h x h j w w K S L w f A q W d Y u d 8 9 c C 8 O w L k = " > A A A B 8 3 i c b V B N S w M x E M 3 W r 1 q / q h 6 9 B I t S P Z T d K u i x 4 M V j B f s B 3 b X M p t k 2 N M k u S V Y o S / + G F w + K e P X P e P P f m L Z 7 0 O q D g c d 7 M 8 z M C x P O t H H d L 6 e w s r q 2 v l H c L G 1 t 7 + z u l f c P  colliders [1][2][3][4][5].However, extractions of gluon TMDs from the Higgs transverse-momentum spectrum is challenging due to the nature of the scalar boson and its large mass (f.i.see also [6] including jet veto considerations and [7] which do not include all TMD effects).Even for the relatively clean process of Higgs spectrum, both the unpolarized and linearly polarized gluon distributions appear in the leading power factorization of the cross-section in the q T /M H expansion, where q T and M H are the boson transverse momentum and mass, respectively.The absence of a color neutral scalar at low energies has driven the attention to quarkonium production both in semi-inclusive deep inelastic scattering (SIDIS) at Electron Ion Collider (EIC) and LHC [3,[8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27].However, the factorization of these processes is also challenging (see [22,23]) and a series of QCD effects are present because of the color structure of quarkonia and the complexity of the non-relativistic expansion (commonly used in quarkonium production studies).
In this work we consider two alternative processes which are presently attracting increasing attention: the dijet [28] and heavy-meson pair [29][30][31] production in an electron-hadron collider, as generated by the (γ * g) and/or (γ * f ) hard interactions.The processes are where and are the initial and final state leptons, h is the colliding hadron, and J i and H/ H are the jets and heavy mesons, respectively.All undetected particles in (1.1) are represented by X. Dijet production has been the object of several studies at the EIC, as it is sensitive to polarized and unpolarized gluon TMDs [32][33][34][35].The produced jets analyzed in the Breit frame have typically a p T ∈ [5,40] GeV and are found in the central rapidity region.Recent studies (see for example [35]) suggest that the experimental observation of the dijet imbalance is possible at the future EIC.The kinematic constraints we consider here for the dijet process need to be such that do not create hierarchies among the partonic Mandelstam variables, i.e. we demand that ŝ ∼ t ∼ û.If such hierarchies exist, they will induce large logarithms in the hard factor of the cross-section and can potentially ruin the convergence of perturbative expansion, unless further resummation/refactorization of the hard factor is performed.
The heavy-meson pair case is instead experimentally more challenging due to the necessity to reconstruct the momenta of the heavy meson from its decay products.In addition, the large energy required to produce a boosted heavy-meson pair makes the process less likely to be observed compared to the dijet production process.On the other hand, recent investigations using monte-carlo generators suggest that for charmed mesons this observable could be possible.The factorization we construct below requires the transverse momenta of the heavy mesons, p H/ H T , be parametrically larger than their mass, m H , i.e. p

H/
H T m H .Although an alternative factorization can be constructed when this condition is violated, we will not pursue this task here.
At leading order (LO) the two processes are schematically shown in fig. 1.In the case of jets we have that the initial parton can be either a gluon or a quark, while in the heavy-meson case only the gluon initial state is relevant.We consider the differential cross-section where x is the Bjorken variable, and η i , r T and p T are respectively the rapidity, the sum of the transverse momenta (with respect to the beam axis) and the average scalar transverse momenta of the two final jets.In the Breit frame, where the virtual photon and target-hadron directions are backto-back, the factorization holds when |r T | p T .The factorization of the cross-section involves the standard TMDPDFs (we have unpolarized and linearly polarized gluon TMDs, and/or quark TMDs), jet or heavy meson distributions, and a new TMD soft function built with Wilson lines aligned along the directions of the incoming hadron and the two outgoing jets.This three-direction soft function has some similarities with the one found in vector boson + jet processes in hadronic colliders [36,37], however the structure of rapidity divergences is very different compared to the soft function discussed in those studies.The perturbative calculation of this new TMD soft function is performed here at oneloop using the modified δ-regulator introduced in [38,39], and we explicitly check the consistency of the factorization at the same order.For the photon-gluon-fusion channel, higher orders of the anomalous dimension of the new soft function can be deduced from the consistency of the anomalous dimensions of the factorized cross-section, since for the case of heavy-meson pair production all other pieces of the cross-section are also known at higher orders.
The presence of the new TMD soft function raises the question of the universality of TMDs beyond the conventional processes of Drell-Yan, semi-inclusive DIS (SIDIS) and di-hadron production in electrorn-positron annihilation1 .While the non-perturbative evolution is universal (in the processes we are considering here and other processes such as SIDIS), the non-perturbative corrections to the soft matrix element are not yet connected to other processes.It is therefore non-trivial to independently extract the gluon TMDPDFs.However, the size of universality breaking effects can be estimated phenomenologically by comparison with simpler processes2 .For a quantitative analysis of these effects further theoretical advancements are needed.
The paper is organized as follows: in sec. 2 we set the notation to be used in the rest of the paper and we give the dijet process factorization theorem.In sec. 3 we extend the discussion to the heavy-meson pair production and we comment on the universality of the heavy-quark jet functions that appear in our factorization theorem, and the corresponding fragmentation shape functions used to describe heavy-meson fragmentation at threshold.Finally we conclude in sec.4. In appendices A and B we collected known results from the literature which we use.In appendix C we give a pedagogical review of some loop calculations made in this work.

Dijet imbalance
In this section we discuss the factorization of the cross-section for the dijet case in DIS within the softcollinear effective theory (SCET).We do not give a detailed derivation of the factorization theorem, but we rather summarize the final result.We also present the NLO calculation for the new threedirection soft function and perform a consistency check of our results using the invariance of the cross-section under renormalization group evolution.The notation and kinematics that we develop here are also useful for the heavy-meson pair production presented in the subsequent section.

Notation and kinematics
Assuming that the direction of the beam is along the ẑ axis it is useful to define the four-vector We also define a conjugate vector nµ by reversing the sign of the spacial coordinates.Thus, n µ and nµ satisfy, Using the vectors n µ and nµ we can decompose any other four-vector, p µ , into its light-cone components, where For the direction of the two jets we use v 1 and v 2 , normalized as where the conjugate vectors vJ , as above, are defined by reversing the sign of the spacial components.We define the standard Lorentz-invariants, where q µ is the momentum of the virtual photon, P µ is the momentum of the target hadron.In the Breit frame we have q µ = (0, 0, 0, Q) and neglecting mass corrections we can solve for target hadron momentum, The fraction of the longitudinal momenta of the incoming parton and the target hadron we denote with ξ, where k µ the momenta of the parton incoming to the hard process.We can then express the variables Q and ξ in terms of the Born level kinematics using the pseudo-rapidities, η 1 and η 2 , and the transverse momentum, p T , of the two outgoing partons, where In this expressions we have neglected corrections from the target hadron mass.The partonic Mandelstam variables in terms of the same quantities are, where p µ 1 and p µ 2 are the momenta of the outgoing partons.It is easy to check that the partonic Mandelstam variables satisfy, (2.12) We denote the transverse momentum imbalance of the two jets with r T , where the hard transverse momentum p T corresponds, up-to power corrections, to the average transverse momenta of the two jets, where the sub-index 1,2 refers to the final jets.At Born level p 1T = −p 2T and thus r T = 0.However, hadronization of the outgoing partons will form jet-like configurations along similar directions and wide angle radiation could escape the jet clustering algorithm, which will then contribute to the imbalance.

Factorization theorem for dijet production
There are two channels for the dijet process that we need to consider: a) the gluon-photon fusion channel which corresponds to the partonic process, γ * g → f f , and b) the incoming quark or antiquark channel from the partonic process, γ * f → gf .The factorization that we propose holds when |r T | p T and we treat the cross-sections only at leading power in an |r T |/p T or |r T |/Q expansion.Within this constraint, the gluon-photon channel factorized cross-section is where the sum runs over all light quark and antiquark flavours f .Here we consider jets with their momentum reconstructed with the so called E-scheme, that is, the momentum of the jets is given by the sum of all jet-constituents.For these jets and for small jet radius (R 1) the cross-section can be factorized in terms of the collinear-soft function C i (b, R), that describes the soft radiation close to the jet boundary and the exclusive jet functions J i , that describe the collinear and energetic radiation confined within the jet. 3 These functions are calculated up to NLO for generic k T -type and cone jet algorithms in [36,44].The corresponding operator definitions are given in the appendix A. In addition, the factorization theorem contains the dijet soft function S γg , which is discussed and calculated at one-loop in next section.Finally, we have the gluon TMD function F g , whose operator definition can also be found in appendix A.
We notice at this point that the latter two functions, S γg and F g , have an intricate interplay due to the rapidity divergences.In sec.2.4 we explain this issue and the role played by the zero-bin subtractions, which lead to the proper definition of the dijet soft function and the TMDs.
The gluon TMD for an unpolarized proton can be further separated into two pieces, the unpolarized gluon distribution f 1 (ξ, b) and the linearly polarized gluon contribution h ⊥ 1 (ξ, b): with g µν T = g µν −n µ nν − nµ n ν .Both of these functions are known perturbatively up to next-to-next-to leading order (NNLO) [5,39,45].The evolution of the TMDs, which is universal (see e.g.[46,47]), is also known up to N 3 LO [38,48,49].In (2.15) we have included only twist-2 TMDs, neglecting higher twists in the TMD expansion.This is sufficient since we consider higher-twist functions suppressed, which is consistent with SIDIS studies as in [50].However, we have no quantitative estimate of these functions and further investigation is important.The inclusion of the higher-twist contributions is beyond the scope of this work and we leave such considerations for future studies.The hard function is then decomposed in two tensor structures: where we have ignored all terms proportional to the four-vector n µ (since they vanish after Lorentzcontraction with the gluon beam function) and any anti-symmetric combinations (since the crosssection is integrated over angles).The coefficients σ gU (L) 0 are introduced such that the leading order hard functions are normalized to the unity, i.e.H U (L) LO = 1 + O(α s ).With this we can now separate the cross-section into a contribution from the unpolarized gluons and one from the linearly polarized gluons.We write the cross-section as where    and We used the shorthand notation s b and c b for the sine and cosine of the angle between the vectors b and v 1T , respectively.The hard factors are calculated up to NNLO in the unpolarized case in [51,52], while for the linearly polarized gluons they are calculated at LO in [53].All hard coefficients are reported in appendix A. Finally, the incoming quark channel has contributions only from the unpolarized gluon jets and the cross-section is given by the following factorized formula: where F f is the f -flavor quark unpolarized TMDPDF.

The dijet soft function at NLO
The only new matrix element in the previous section is the soft function.Here we give the operator matrix element definition of the soft function and we proceed with the NLO calculation.The details of the calculation are collected in the appendix.We start defining the soft function for the photon-gluon fusion process: The soft function corresponding to the case of incoming quark or antiquark is obtained with the exchange Ŝγf = Ŝγg (n ↔ v 2 ) . (2.22) The Wilson lines are defined as We omit here T -Wilson lines which are needed in singular gauges.We also distinguish Wilson lines in the adjoint and fundamental SU(N c ) representations using S and S respectively.Note that the δ-regulator is introduced only in S n (S n ).We write the soft function as a series in a s = α s /(4π) At tree level Ŝ[0] = 1, and at one-loop only the diagrams with a real gluon give a non-zero result.We can therefore write the one-loop soft function as a sum of diagrams with one real gluon exchange between the three soft Wilson lines, where the suffix indicates the Wilson lines connected by the exchanged gluon, including the mirror diagrams as well.Note that the mirror diagrams for ŜJB will introduce an additional iπ component, as discussed in more detail in the appendix.The coefficients C ij are the color factors and they are different for the two channels γ * g and γ * f , The relevant diagrams are shown in fig. 2 and the corresponding contributions are where J = 1, 2 and the corresponding integrals are with d = 4 − 2 .The results at all orders in are and with the shorthand notation, A b is also dimensionless, longitudinal boost invariant, and it is bounded to negative values less than −1, i.e., A b ≤ −1.The function 2 F 1 is the standard hypergeometric function and the expansion for this function can be written as follows Adding all contributions and expanding in we obtain the bare soft function.For the γ * g-channel we have: for the γ * f -channel: where the finite part of the soft function is To simplify the results we have used the notation (2.37)

The zero-bin subtraction and the universal TMDs
Here we reorganize the factorization theorem such that the cross-section is expressed in terms of rapidity divergence-free TMDs as presented in the factorization theorem in (2.18), (2.19) and (2.20).
In this way we make clear the dependence on the universal TMDPDFs and the unknown TMD soft function.To do this we write the TMD beam function in terms of the zero-bin un-subtracted term divided by the back-to-back soft function, The back-to-back two-direction soft function operator definition can be found in appendix A. Then following [38,54] we can factorize the soft function as where ν is an arbitrary positive number which plays a booking role and will be removed from the final result, introducing this way a constraint on the product of rapidity scales.The bare function where here and in the rest of this manuscript we used the shorthand notation Thus we can now reorganize the beam and soft function matrix elements (denoted by the "hat" notation) into a product of TMDs as they appear in the factorization theorem, where the functions in the r.h.s. of (2.42) are respectively that is, the universal TMDPDF as defined in other observables such as Drell-Yan and semi-inclusive DIS and that is the unknown soft function now incorporated in a rapidity divergent-free ratio.We can thus eliminate the dependence on the arbitrary parameter ν by introducing the following constraint where ŝ, t, and û are the partonic Mandelstam variables and thus the combination ζ 1 ζ 2 is Lorentz invariant and in the Breit frame, or any other frame boosted along the proton direction, we have: Notice that the procedure to obtain (2.45) is totally analogous to the one used in Drell-Yan or SIDIS [55,56], TMD factorization theorem.In that case we have that ζ 1,2 have both a square mass dimension and ζ 1 ζ 2 = Q 4 , while in the present case ζ 2 is dimensionless quantity but ζ 1 , as usual, has dimensions of mass squared.The natural way to choose the values of ζ 1 and ζ 2 is This way we have the standard evolution for the TMDPDF up to the hard scale and the ratio of soft functions in (2.44) has no large rapidity logarithms and thus does not require evolution in rapidity.
The renormalized soft function can then be written in terms of the renormalization kernel Z S and the bare soft function which is simply the ratio of the bare functions that appear in (2.44), In the MS scheme for the (γ * g) channel we have and for the (γ * f ) channel we have where the corresponding renormalization functions are and The soft anomalous dimension can be obtained from the renormalization functions as follows, The one-loop results for the soft function anomalous dimensions are collected in the next section.

Consistency check
Each element of the factorized cross-section has a factorization scale dependence and it satisfies a renormalization group equation, where G runs over all the functions in the factorization theorem and γ G is the corresponding anomalous dimension.On the other hand, the cross-section is renormalization group invariant.Therefore, as required by consistency, the sum of all anomalous dimensions of the terms appearing in the factorized cross-section must vanish.In the impact parameter space, where the cross-section is written as a product of these functions, we have, and We write the perturbative expansion of these anomalous dimension as with a s = α s /(4π).For the two channels in the dijet process the relevant anomalous dimensions up to one-loop are, γ The imaginary component in the collinear-soft anomalous dimension is denoted by κ(v i ) where Note that both in the anomalous dimension and the finite terms the imaginary terms cancel in the sum resulting in a real cross-section as a function of b, which we have checked explicitly up to the NLO contributions.These anomalous dimensions can be found in [36,38,44,[51][52][53].We also used (2.31) to expand A b in the soft function anomalous dimension in terms of ŝ, p T , and c b .It is now easy to confirm the cancelation of the anomalous dimensions at O(α s ) which also serves as confirmation of the factorization theorem at the same order.

Heavy-meson pair imbalance
In this section we consider the process of heavy-meson pair production, + p → + H + H + X, in the back-to-back limit for which the transverse momentum imbalance r T is measured, where we use the notation H for generic heavy meson and H for the corresponding anti-particle.The imbalance is measured in the Breit frame and in the region sensitive to TMDs, i.e., |r T | p H, H T , the two heavy mesons are fragmented near the kinematic end-point and carry most of the energy of the heavy quark coming from the hard process.
In contrast to the dijet process, for the heavy-meson pair production we only need to consider the photon-gluon fusion channel.From this perspective the formalism is simpler, but on the other hand the mass of the heavy meson, m H , introduces a new scale which we need to consider.In the case when the heavy mesons are highly boosted, i.e., p H T m H , the factorization is similar to the dijet production discussed above.The cross-section is then expressed in terms of the same hard, soft, and beam functions, but the production of the final state heavy mesons is described by a heavy-quark jet function, J Q→H [57,58], with η H and η H referring to the pseudo-rapidities of the heavy mesons.Similarly to the dijet factorization, in the hard function we do not consider corrections due to the quark mass and we define, The decomposition into the unpolarized and linearly polarized gluon contributions follows the same steps as in section 2 and thus we do not repeat here.

Refactorization of heavy-quark fragmentation function
The fragmentation of a heavy quark to a heavy meson is described by the heavy-quark fragmentation function which is studied in a plethora of processes (see for example [29][30][31]).The large scale of the process, which is introduced by the mass of the heavy quark, allows for the use of perturbation theory to calculate the fragmentation function up to small and universal non-perturbative corrections.
In our case the heavy-quark jet function, J Q→H , describes the fragmentation of heavy mesons from heavy quarks and is differential in the two-dimensional transverse momentum of the fragments w.r.t. the beam axis.In the limit r T p T there are two parametrically different scales which are involved in the fragmentation process, Logarithms of ratios of these scales will appear in the perturbative calculation of the jet function and could potentially ruin the convergence of perturbative expansion.Thus, resummation of these logarithms is essential to ensure the convergence of the expansion.Note that these logarithms are the same to the logarithms resummed by the TMD evolution.The resummation of the logs generated by the scales in (3.4) can be achieved through the means of factorization.To this end we employ the boosted-heavy-quark effective field theory (bHQET) [59] which will allow us to factorize the jet function into a hard matching coefficient and a transverse momentum dependent matrix element.To demonstrate how such a factorization occurs, we give a brief description of the relevant modes.
First consider the momentum of the heavy quark in the heavy meson rest frame, p µ Q , which can be decomposed into a mass term and the residual soft component, where k µ s is the typical size of soft (light) degrees of freedom in the heavy meson and β µ = (1, 1, 0 ⊥ ) v .Note that here we decomposed four-vectors into the light-cone coordinates along the direction of the boosted heavy meson, v.The momenta k µ s , in the boosted frame sets the size of energy loss of the heavy quark during fragmentation to a heavy meson.To obtain the boosted momenta we simply apply the following transformations, We can obtain Λ by comparing the momentum of the heavy quark at the rest frame of the heavy meson in (3.5) to the momentum of the boosted heavy quark (up to power-corrections ∼ m H /p H T ), and thus Λ = 2E H /m H .The transformations in (3.6) give the momentum scaling of the so called "ultra-collinear" modes which are simply the soft modes of HQET boosted to the frame where we perform the measurement, The contribution to the transverse momentum spectrum w.r.t. the beam axis, comes from the large component v •k uc .Therefore, we can estimate the typical size of the transverse momentum imbalance, T m H , the typical soft scale, r T , is perturbative and that justifies the approach of perturbative matching to the TMD matrix elements to which non-perturbative effects are incorporated as corrections.
To proceed with the factorization we match from the massive-SCET [59,60], which includes collinear degrees of freedom, onto the boosted HQET where the degrees of freedom are the ultracollinear modes.We are interested in the matching of the massive, collinear gauge invariant, quark building block, onto the HQET heavy quark fields, where C + (m Q , µ) is the short distance matching coefficient and β + denotes the heavy quark velocity in the boosted frame.With this matching we can now factorize the jet function into a short distance matching coefficient and a bHQET matrix element that depends on the transverse momentum of the ultra-collinear fragments, where The operator definition of the two-dimensional shape function is where v is a Euclidean, two dimensional, transverse component of light-like four-vector v µ pointing along the direction of the boosted heavy meson.The impact parameter space expression is obtained by simply taking the Fourier transform, The hard matching coefficient is known up to two-loops but the jet function J Q→H (r), as defined above, appears here for the first time.However, as we discuss later in this section, this jet function is related at the operator level to the fragmentation shape function from [57,58] in the near-end-point limit, (z H → 1).The one-loop hard function, H + is, and the corresponding anomalous dimension is In the following section we show the calculation of the bHQET matrix element, J Q→H at NLO and we use this result to derive the one-loop anomalous dimension.We demonstrate the consistency of anomalous dimensions for this process at NLO and we give an all order statement that connects the matrix element J Q→H to the near-end-point fragmentation shape function for heavy mesons.

The bHQET matrix element at NLO
The one-loop contributions to J Q→H (b) in (3.13) are shown in fig. 3. We have non zero contributions only from diagrams (f) and (h) where a real gluon is exchanged.Virtual diagrams (e) and (g) are   scaleless and vanish in our scheme, where where we have made the standard M S scale replacement µ 2 → µ 2 exp(γ E )/(4π).Expanding (3.16) in the limit → 0 and keeping all poles and finite non-vanishing terms we have The renormalized jet function which in the MS-scheme is simply given by the finite terms in (3.18) is defined by the following equation, The corresponding anomalous dimension is It is now trivial to show the consistency of the anomalous dimensions and therefore of factorization at NLO, since it is sufficient to show, The result in (3.18) involves logarithms of the scale m Q /(p T b 0 ), where b 0 = b exp(γ E )/2.On the other hand the hard function H + in (3.14) involves only logarithms of m Q .This suggests that we have successfully separated the two scales and we can resum ratios of those (i.e., logarithms of p T b 0 ) by evaluating each function at its canonical scale and then evolving each up to a common scale by solving with the corresponding renormalization group equation.

Connection to the fragmentation shape function
In this section we show that the two dimensional bHQET jet function in (3.12) is related to the fragmentation shape function.We consider the operator definition of the shape function as in eq.(4.23) of [58].The same same shape function was written with a different notation in [57].
Taking the one dimensional Fourier transform of this expression with respect to ω we have5 Comparing this result with the two dimensional Fourier transform of (3.12) as prescribed in (3.13) we have Using this equation we can confirm our perturbative calculation at NLO in (3.18) by comparing the finite terms of this equation against eq.(4.24) of [58].To do that one needs the Fourier transformations of the regular plus-distributions which can be found in the literature but we also give here for completeness, (3.25) Using these equations we can easily check that indeed (3.24) is satisfied up to NLO, although the perturbative validity of (3.24) is inferred beyond NLO.Therefore, since the anomalous dimension of the fragmentation shape function and the hard function H + is already known up to two-loops, we can use the consistency of factorization to solve for the anomalous dimension of the global soft function S γg up to two-loops, The anomalous dimensions at two and three loops are given in appendix B. We give the result for the soft function by organizing it into a term proportional to the cusp anomalous dimension, γ cusp and a "non-cusp" term, where and we have used the same notation for the perturbative expansion of the cusp anomalous dimension as in (2.56).The lengthy and not so intuitive three-loop non-cusp component, δγ [3] Sγg , is given in the appendix (see (B.13)).
With this result we can push the calculation of the heavy-meson pair production up to NNLL with no additional perturbative calculations.Furthermore with the knowledge of the soft anomalous dimension and using (2.54) we can also solve for the collinear-soft anomalous dimension to the same order.This can now give us the NNLL cross-section of the dijet photon-gluon-fusion (γ * g) process.
For the full NNLL dijet cross-section we are still missing the global-soft or collinear-soft anomalous dimensions from the photon-quark-initiated (γ * f ) process.

Conclusions
In this work we have established a new factorization theorem for dijet and heavy-meson pair production in DIS which can be valuable in the quest of processes with a clear sensitivity to gluon TMDs.The factorization involves a new soft function, which we have calculated at one-loop, and whose anomalous dimension has been deduced at two and three loops from consistency relations.All the calculations have been performed with the δ-regulator, combined with standard dimensional regularization.The factorized cross-section is then written terms of TMD parton distribution functions, the new TMD soft function and two final-state jet functions or heavy hadron distributions.The cross-section is sensitive to both unpolarized and linearly polarized gluon TMDs.
The influence of this new soft function is certainly an element that should be studied in the future.In particular one should understand how large is its non-perturbative contribution to the cross-section and whether it appears in multiple processes, i.e. whether it is a universal quantity.For a recent discussion on the universality of TMDs with multiple collinear directions see [61].
In our dijet analysis we do not consider the effects of any possible non-global logarithms which could be generated from in-out of jet correlations of the collinear-soft modes and are not associated with any of the TMD matrix elements (TMD-soft and TMD-PDF).Also we expect their effect to the resummed cross-section to be particularly small for the kinematic region of interest: p T ∈ [5,40] GeV and for small jet radius R ∼ 0.4 [37].Thus, these effects could be incorporated into the "jetsmearing" effects, from the hadronization of the jets, which are expected to be larger or of the same size.To this end, other possible extensions of this work can improve on this aspect by implementing modern jet substructure techniques, such as grooming, to reduce the sensitivity of the jets to nonglobal logarithms and hadronization effects.Recently in [53] the angular de-correlation between a color-singlet boson (γ, Z, W µ ) and the winner-take-all (WTA) axis was studied in hadronic collisions, and it was shown to be free from non-global logarithms and, in addition, to have small sensitivity to the choice between charged-particles-only (tracks) or full jets.This last property of the angular decorrelation measurement using the WTA axis can be particularly useful when experimental limitations exist on the reconstruction of full jets.Therefore, extensions of [53] in dijet process in DIS are of great interest to both jet and gluon TMD studies.
This study focusses on the theoretical framework and the necessary elements for the resummed cross-section.The details of a numerical study can depend on various aspects, such as the schemes for the TMD evolution and treatment of power corrections.We postpone a more quantitative analysis for a future study.In addition, a natural extension of this work is to incorporate spin effects from a polarized target hadron.Such effects will give rise to spin asymmetries, and particularly interesting is the case of the Sivers asymmetry (recently studied in hadronic dijet production [62]).Thus, a simple generalization of our formalism can help formulate a factorization framework for processes with sensitivity to the gluon Sivers function.

B Anomalous dimensions
In this appendix we collect from the literature the anomalous dimensions of the elements of the crosssection needed to obtain the dijet soft function anomalous dimension at two and three loops for the γ * g-channel, as explained in (3.26).We follow the standard procedure of separating all anomalous dimensions into a term proportional to the cusp anomalous dimension and a non-cusp term.We first give the three-loop cusp which applies to all functions and then proceed to give the non-cusp term for each function separately.The anomalous dimensions we give here are the ones that correspond to the renormalization group equation in (2.53).

B.1 The cusp anomalous dimension and β function
where and for i = f, g we have C f = C F and C g = C A .The β-function is given by where The non-cusp anomalous dimensions for the hard function and the jet function are known up to two-loops [58] and are given by δγ Their sum is known up to three loops and is given by δγ

B.3 TMDPDF
The universal TMDPDF anomalous dimension is given in [39] by The non-cusp anomalous dimension for the TMDPDF is known up to three loops and is given by

B.4 Hard Function
The hard function H γ * g anomalous dimension is given in [51] by The non-cusp anomalous dimensions is given by δγ V . (B.12)

B.5 Dijet soft function
The dijet anomalous dimension for the γ * g-channel is given by eq.(3.27).The non-cusp three-loop anomalous dimension is given by δγ

C Dijet soft function integrals
In this section we give a pedagogical review of the integrals needed to obtain the dijet soft function result at one-loop order.The integrals are introduced in (2.28).The real diagrams are shown in fig. 2 and are the only ones contributing to the soft function.The virtual diagrams vanish, as we show in this section, and are not shown in this work.In the following, we use d = 4 − 2 .

Real diagram n -v J
The integral we need to compute is given by the expression Using the delta to integrate over k − we get Completing the square, the denominator can be written as We change variables k → k − k + v + J v J .In this way, the integral simplifies to the following expression, This allows us to perform each integral separately, which leads us to the final result of the integral, The integral we need to compute is given by Using the delta to integrate over k − we get We use Feynman parametrization in order to rewrite the denominator, (C.8) In our case, we identify (C.9) In this way, the denominator can be rewritten the following way: where We can perform a change of variables k = k − k + R µ 1 (x) and the integral simplifies as follows, (C.12) We perform a Mellin-Barnes transformation in order to be able to integrate over k and k + separately, (C. 13) In this way, we integrate over k and get The integral that has to be computed is given by First, we integrate over k − .We get a different result depending on the sign of k + , which leads to the addition of θ(k + ) and θ(−k + ), , (C.22) , (C.23) Both integrals are computed the same way, so we focus in one of them, I 1 for example.Here, we use Feynman parametrization in order to rewrite the denominator, (C.25) In our case, we identify In this way, the denominator can be rewritten the following way: where R 1 (x) = x, (C.28) We can perform a change of variables k = k − k + E 1 (x) − E 2 (x) and the integral simplifies to (C.29) We can integrate over k we get Finally, integrating over k + we have which evaluates to zero if we perform the integration over x.The term I 2 is zero for the same reason, we can compute it following the same steps.This means

1 <µ 2 <
2 j p O F a E t E v N Y d U P Q l D N J W 4 Y Z T r u J o i B C T j v h + G b m d x 6 p 0 i y W 9 2 a S 0 E D A U L K I E T B W 8 q v + E I S A h 3 M c n f X L F b f m z o H / E i 8 n F Z S j 2 S 9 / + o O Y p I J K Q z h o 3 f P c x A Q Z K M M I p 9 O S n 2 q a A B n D k P Y s l S C o D r L 5 z V N 8 Y p U B j m J l S x o 8 V 3 9 O Z C C 0 n o j Q d g o w I 7 3 s z c T/ v F 5 q o u s g Y z J J D Z V k s S h K O T Y x n g W A B 0 x R Y v j E E i C K 2 V s x G Y E C Y m x M J R u C t / z y X 9 K u 1 7 y L W v 3 u s t I 4 z e M o o i N 0 j K r I Q 1 e o g W 5 R E 7 U Q Q Q l 6 Q i / o 1 U m d Z + f N e V + 0 F p x 8 5 h D 9 g v P x D Y H X k J g = < / l a t e x i t > l a t e x i t s h a 1 _ b a s e 6 4 = " t G o j q v U + N Y c C L W a r x z 1 q 7 R g c S i o = " > A A A B 8 H i c d V B N S w M x E J 2 t X 7 V + V T 1 6 C R b F 0 7 J b R T 0 W v H i s Y D + k X U s 2 z b a h S X Z J s k J Z + i u 8 e F D E q z / H m / / G t F 1 B i z 4 Y e L w 3 w 8 y 8 M O F M G 8 / 7 d A p L y y u r a 8 X 1 0 s b m 1 v Z O e X e v q e N U E d o g M Y 9 V O 8 S a c i Z p w z D D a T t R F I u Q 0 1 Y 4 u p r 6 r Q e q N I v l r R k n N B B 4 I F n E C D Z W u k t 6 / n 3 W F e m k V 6 5 4 r j c D 8 t z z B e L n V g V y 1 H v l j 2 4 / J q m g 0 h C O t e 7 4 X m K C D C v D C K e T U j f V N M F k h A e 0 Y 6 n E g u o g m x 0 8 Q U d W 6 a M o V r a k Q T P 1 5 0 S G h d Z j E d p O g c 1 Q L 3 p T 8 S + v k 5 r o M s i Y T F J D J Z k v i l K O T I y m 3 6 M + U 5 Q Y P r Y E E 8 X s r Y g M s c L E 2 I x K N o T v T 9 H / p F l 1 / V O 3 e n N W q R 3 n c R T h A A 7 h B H y 4 g B p c Q x 0 a Q E D A I z z D i 6 O c J + f V e Z u 3 F p x 8 Z h 9 + w X n / A t v n k F w = < / l a t e x i t >p l a t e x i t s h a 1 _ b a s e 6 4 = " Z T W N n i B W a n f 5 a e q n 0 r G M T i q U h 0U = " > A A A B 8 H i c d V B N S w M x E J 2 t X 7 V + V T 1 6 C R b F 0 7 J b R T 0 W v H i s Y D + k X Us 2 z b a h S X Z J s k J Z + i u 8 e F D E q z / H m / / G t F 1 B i z 4 Y e L w 3 w 8 y 8 M O F M G 8 / 7 d A p L y y u r a 8 X 1 0 s b m 1 v Z O e X e v q e N U E d o g M Y 9 V O 8 S a c i Z p w z D D a T t R F I u Q 0 1 Y 4 u p r 6 r Q e q N I v l r R k n N B B 4 I F n E C D Z W u k t 6 1 f u s K 9 J J r 1 z x X G 8 G 5 L n n C 8 T P r Q r k q P f K H 9 1 + T F J B p S E c a 9 3 x v c Q E G V a G E U 4 n p W 6 q a Y L J C A 9 o x 1 K J B d V B N j t 4 g o 6 s 0 k d R r G x J g 2 b q z 4 k M C 6 3 H I r S d A p u h X v S m 4 l 9 e J z X R Z Z A x m a S G S j J f F K U c m R h N v 0 d 9 p i g x f G w J J o r Z W x E Z Y o W J s R m V b A j f n 6 L / S b P q + q d u 9 e a s U j v O 4 y j C A R z C C f h w A T W 4 h j o 0 g I C A R 3 i G F 0 c 5 T 8 6 r 8 z Z v L T j 5 z D 7 8 g v P + B d 1 x k F 0 = < / l a t e x i t > q µ = Q 2 (n µ nµ ) = (0, 0, 0, Q) < l a t e x i t s h a 1 _ b a s e 6 4 = " d R b W C J P V W g X N l + C u G n O u / G 5 j D 5 8 = " > A A A C I 3 i c d V B P S 8 M w H E 3 n v z n / T T 1 6 C Q 5 l A x 3 t n E 4 E Y e D F 4 w b O D d Y 5 0 i z d w t K 0 J q k w S r + L F 7 + K F w / K 8 O L B 7 2 K 2 V Z i i L w R e 3 n s / k j w n Y F Q q 0 / w w U g u L S 8 s r 6 d X M 2 v r G 5 l Z 2 e + d W + q H A p I F 9 5 o u W g y R h l J O G o o q R V i A I 8 h x G m s 7 w a u I 3 H 4 i Q 1 O c 3 a h S Q j o f 6 n L o U I 6 W l b v b i / i 6 y v T C G l 9 B 2 B c J R P Y 5 K c Z 4 n 6 j G 0 H S Q i H s / O B R 3 L m 0 e T V S 9 0 s z m r a E 4 B z e J p 5 a x c s T R J l G 8 r B x L U u t m x 3 f N x 6 B G u M E N S t i 0 z U J 0 I C U U x I 3 H G D i U J E B 6 i P m l r y p F H Z C e a / j G G B 1 r p Q d c X e n M F p + r 8 R I Q 8 K U e e o 5 M e U g P 5 2 5 u I f 3 n t U L n n n Y j y I F S E 4 9 l F b s

Figure 1 .
Figure 1.Example LO diagrams for the two processes.The momenta q µ and k µ (corresponding to the photon and incoming parton momenta respectively) are expressed in the Breit frame.
z s e 8 d c X J Z 4 7 g D 5 z P H 1 F F j R I = < / l a t e x i t > (+1) v 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " 5 L 0 b + k l c p 0 v J L H e + F r 0 9 Y 2 a 6 v 7 Y = " > A 2 2 0 2 7 d L M J u 5 t C D f 0 l X j w o 4 t W f 4 s 1 / 4 7 b N Q a s P B h 7 v z T A z L 0 g 4 U 9 p x v q z C 2 v r G 5 l Z x u 7 S z u 7 d f t g 8 O 2 y p O J a E t E v N Y d g O s K G e C t j T T n H Y T S X E U c N o J x r d z v z O h U r F Y P O h p Q r 0 I D w U L G c H a S L 5 d r l 7 0 m Q j 1 9 N z P J n 5 9 5 t s V p + Y s g P 4 S N y c V y N H 0 7 c / + I C r S e r T f r f d l a s P K Z I / g F 6 + M b 7 0 a S i Q = = < / l a t e x i t > ( 1) n < l a t e x i t s h a 1 _ b a s e 6 4 = " w h 9 x Q 4 Y b X f s M R 5 a 5 b 1 u k t B + S N + g = " > A A A B 9 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 B I t S D 5 a k C n o s e P F Y w X 5 A G 8 p m u 2 m X b j Z x d 1 I o o b / D i w d F v P p j v P l v 3 L Y 5 a O u D g c d 7 M 8 z M 8 2 P B N T r O t 5 V b W 9 / Y 3 M p v F 3 Z 2 9 / Y P i o d H T R 0

Figure 2 .
Figure 2. Tree level soft function is shown in diagram (a).Each double line represents a Wilson line whose pointing direction is also reported as (±∞) u , u = n, v, v; a and a are color indices.Diagrams (b), (c) and (d) also contribute to the NLO soft function.Virtual contributions vanish and therefore they are not shown here.Mirror diagrams are also not shown.
) and b = (0, 0, b)/|b|.Note that A b is a function of the angle φ bJ only and it does not depend on |b|.
s h a 1 _ b a s e 6 4 = " / A C k 3 H E s t w h 3 2 i W 7 H v 8 w 4 x 4 + r 1 Q = " > A A A B 8 3 i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b F U 0 m q o M e C F 4 8 V T C 0 0 s W w 2 m 3 b p Z h N 2 N 4 U S 8 j e 8 e F D E q 3 / G m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I O V M a d v + t i p r 6 x u b W 9 X t 2 s 7 u 3 v 5 B / f C o q 5 J M E u q S h C e y F 2 B F O R P U 1 U x z 2 k s l x X H A 6 W M w v p 3 5 j x M q F U v E g 5 6 m 1 I / x U L C I E a y N 5 L m D f F I 8 5 V 6 I h 8 W g 3 r C b 9 h x o l T g l a U C J z q D + 5 Y U J y W I q N O F Y q b 5 j p 9 r P s d S M c a h b Y z o j g y q 9 5 c / M / r p D i 4 C z K h k h S 5 Y s t F g 1 Q S j M n 8 c 9 I X m j O U U 0 s o 0 8 L e S t i I a s r Q 5 l O y I X i r L 6 + T Z q 3 q X V d r j z e V + m U e R x H O 4 B y u w I N b q M M D N M A H B g K e 4 R X e H O W 8 O O / O x 7 K 1 4 O Q z p / A H z u c P / G m O t g = = < / l a t e x i t > h v ( + ) < l a t e x i t s h a 1 _ b a s e 6 4 = " R 3 l w 8 y R y v X Q J w C O M L j q G S c 5 2 P o I = " > A A A B 9 X i c d V B N S 8 N A E N 3 U r 1 q / q h 6 9 L B a l I p S k i n o s e P F Y w X 5 A G 8 N m O 2 m X b j Z h d 1 M p o f / D i w d F v P p f v P l v 3 L Y R t O i D g c d 7 M 8 z M 8 2 P O l L b t T y u 3 t L y y u p Z f L 2 x s b m 3