Transverse momentum dependent distributions in dijet and heavy hadron pair production at EIC

We discuss the measurement of gluon transverse momentum distribution (TMD) in dijet and heavy hadron pair (HHP) production in semi-inclusive deep inelastic scattering. The factorization of these processes in impact parameter space shows the appearance of a specific new soft factor matrix element on top of angular a complex valued anomalous dimensions. We show in detail how these features can be treated consistently and we discuss a scale prescription for the evolution kernel of the dijet soft function. As a result we obtain phenomenological predictions for unpolarized and angular modulated cross-sections for the electron-ion collider (EIC) using current available information on unpolarized TMD.


Introduction
The access to non-perturbative gluon distributions from experiments is notoriously challenging. This is also the case of gluon transverse momentum distributions (TMDs). Gluons enter directly in Higgs production in hadronic colliders [1][2][3][4][5] that has a relatively high mass and low production rates, and quarkonium production both at EIC and LHC [3,[6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25] that is sensitive also to the heavy quark hadronization effects [20,21]. Recent studies (see for example [26,27]) suggest that the experimental observation of the dijet imbalance is possible at the future EIC. In a recent work [28] we have proposed the dijet and hadron pair production at electron-ion colliders (EIC) to probe gluon TMD. Previous studies on these JHEP03(2022)047 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

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 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 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 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 / + o O Y p I J K Q z h o 3 f P c x A Q Z K M M I p 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 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 = " x n 1 t + / h v 2 J P Z 1 N 8 L v P o y g r k y processes were performed in refs. [29,30] for dijet and in refs. [31][32][33] for heavy-meson pair production in an electron-hadron collider. The relevant processes for our case are where and are the initial and final state leptons, h is the colliding hadron, J i and H/H are the jets and heavy mesons respectively and X represent undetected particles. At leading order (LO), and ignoring the intrinsic momentum of partons inside the target hadron, the two hard-scattering processes are schematically shown in figure 1. We have shown that these processes are factorizable when considering the cross-section dσ dxdη 1 dη 2 dp T dr T , (1.2) 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. The factorization condition is |r T | p T in the Breit frame, as the virtual photon and target-hadron directions are back-to-back. These conditions are certainly fulfilled in the Breit frame for p T ∈ [5,40] GeV and in the central rapidity region. We also demand that there are no hierarchies among partonic Mandelstam variables,ŝ ∼ |t| ∼ |û|. On the experimental side recent studies [26] suggest that the measurement of dijet imbalance is possible at EIC. For the heavy-meson case monte-carlo generator studies suggest that charmed mesons can be reconstructed [34,35]. The charm production rates have been considered at the LO and NLO QCD for ep → c/c + X in ref. [36]. For our case in principle we require the transverse momenta of the heavy mesons, At leading order (LO) the jet processes can be initiated by either a gluon or a quark, while in the heavy-meson case only the gluon initial state is relevant. The factorized cross-section results in products/convolutions of several fundamental functions as TMDs, jet/heavy quark distributions, soft function and a new evolution kernel, detailed in [28]. The purpose of the present paper is to provide a phenomenological study of these processes,

JHEP03(2022)047
including theoretical errors. In order to achieve this, we have used the code Artemide [37,38], introducing new moduli necessary to describe the present cases. The code already includes quark TMDPDF and the TMD evolution kernel extracted from Drell-Yan and semi-inclusive DIS experiments [39].
The factorization of the cross-section is produced in position space, described by the variable b, conjugate of the momentum r T . In order to Fourier transform the factorized cross-section from b space to momentum space one has to perform an angular integration that results non trivial because the anomalous dimensions of several functions depend on this angle and are complex valued. We show that this integration can be performed in resummed perturbation theory and that this addresses the problem of complex values in all anomalous dimensions. As a result the φ b -angle integrated cross-section is then factorized into distributions that are b-rotation invariant. The φ b -angle integrated dijet evolution kernel is derived integrating a system of coupled differential equations similar to the TMD case. We show and discuss here a specific scale choice prescription that is analogue to the ζ-prescription already discussed in [40], which was already implemented in Artemide.
The factorization theorem and the detailed definition of the observables is provided in section 2. The cross-sections subtleties discovered when comparing to other groups are discussed in section 3. The resummation of logs involves angular integrations as explained in section 4, which leads to the evolution kernels detailed in section 5. The phenomenological results obtained with the codes that we have developed are summarized in section 6, after which the conclusions are drawn.

Notation and kinematics
In order to define angles and to deduce a factorized cross-section we need to establish some kinematics. The direction of the beam is fixed along theẑ axis. It is useful to define the four-vectors n µ = 1 √ 2 (1, 0, 0, 1),n µ = 1 √ 2 (1, 0, 0, −1), so that n 2 =n 2 = 0,n · n = 1. Then any other four vector can be decomposed into its light-cone components, with The direction of the two jets are v 1 and v 2 , normalized as andv J are defined by reversing the sign of the spacial components. We have then the standard Lorentz-invariants,

JHEP03(2022)047
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 P µ = 1 2x (Q, 0, 0, −Q) . The ratio of the longitudinal momenta of the incoming parton and the target hadron is ξ = k + P + .with k µ the momentum of the parton entering 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, neglecting corrections from the target hadron mass, η ± = η 1 ±η 2

2
. The partonic Mandelstam variables can be written using the same variables, with p µ 1 and p µ 2 the momenta of the outgoing partons. At partonic level they satisfŷ Finally, the transverse momentum imbalance of the two jets, r T , and the hard transverse momentum, p T , are defined through where the sub-index 1, 2 refers to the final jets. At Born level p 1T = −p 2T and thus r T = 0. It must be taken into account that the hadronization of the outgoing partons will form jet-like configurations along similar directions and wide angle radiation that can escape the jet clustering algorithm, affecting the imbalance.

Factorization theorem for dijet and heavy hadron pair production
The factorization of dijet and heavy hadron pair production at leading power (LP) for semiinclusive deep inelastic experiments has already been provided in [28]. The cross-sections reported here do not take into account any leptonic fiducial cuts, which however could be implemented once the experimental conditions are established (especially at EIC). In this section we recall the main formulas that are used in our phenomenological description. We start with the dijet cross-section which can be written as a sum of terms depending on the parton that initiates the hard process (quark or gluon) dσ(γ * g) = dσ U (γ * g) + dσ L (γ * g). (2.10)

JHEP03(2022)047
For an unpolarized hadronic process, the cross-section is made out of hard contributions from unpolarized initial quarks dσ U (γ * f ), unpolarized initial gluons dσ U (γ * g), and linearly polarized gluons dσ L (γ * f ). The quark contribution to the cross-section is In this formula f f 1 is the unpolarized quark TMDPDF for flavor f , H U the hard factor for the unpolarized quark case. The perturbative calculations of TMDPDF has been performed recently at NNLO [5,[41][42][43][44] and N3LO [45,46]. The jets are described by the product of a collinear-soft function C (f,g) and a jet shape function J (f,g) specific for each partonic flavor. The calculation at NLO of these functions can be found in [47,48] for generic k T -type and cone jet algorithms.
The factor S γf is the dijet soft function for the fundamental representation of SU(3) C and calculated in [28] up to NLO. A corresponding soft factor, S γg , for the adjoint representation of SU(3) C is also necessary for the incoming gluon contribution.
The hard factor H µν (µ) accounts for contributions of unpolarized and linearly polarized gluons, The TMD tensor F g,µν can be also decomposed in terms of unpolarized and linearly polarized parts, with g µν T = g µν − n µnν −n µ n ν . The hard factors are evaluated up to NNLO in the unpolarized case in [49,50] and at LO for the linearly polarized case [51]. In these equations f g 1 and h ⊥ 1 represent the unpolarized and linearly polarized gluon TMD. Both of them are known perturbatively up to NNLO [5,42,44]. Combining eq. (2.10), (2.12), (2.13), (2.14) We use s b = sin φ b and c b = cos φ b for the sine and cosine of the angle φ b between the vectors b and v 1T , respectively. Each of dσ has a hard factor that describes the initiating interaction. The coefficients σ The case of heavy hadron pair is very similar. The measured imbalance r T is where the superscript H indicates a generic heavy meson andH the corresponding antiparticle. The imbalance is measured in the Breit frame and assuming the TMD factorization scaling, i.e., |r T | p H,H T . We also assume that 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. The cross-section reads with η H and ηH the pseudo-rapidities of the heavy mesons, J Q→H the heavy quark jetfunctions [52,53]. The hard, soft, and beam functions are the same as in the dijet case. In the hard function we do not consider corrections due to the quark mass and we define The heavy quark jet functions, J Q→H , can be partially evaluated in perturbation theory as shown in [28]. We work in the limit p T m H Λ QCD and the heavy quark jet function can be re-factorized using bHQET. We also have r T p T so that it is possible to find large logs of two parametrically different scales in the fragmentation process,

JHEP03(2022)047
that need to be resummed to ensure the convergence of the expansion. Following [28] the jet function can be firstly factorized into a short distance matching coefficient and a bHQET matrix element, and the two-dimensional shape function is defined in momentum space as Notice that v is a Euclidean, two dimensional, transverse component of the light-like fourvector v µ pointing along the direction of the boosted heavy meson. In position space J Q→H is obtained by Fourier transformation The one-loop expression for these quantities are calculated in [28].

Cross-sections used in phenomenology
The cross-sections presented in previous section are usually partially integrated in phenomenological observables. We discuss here these integrations, which also allow us to relate the normalization of our cross-section with the ones obtained in the literature.

Extracting the Born-level cross-sections
The tree level cross-sections for the dijet and hadron pair production were considered at tree level in ref. [54]. We start considering the gluon case, from which one can easily deduce also the quark case. The gluon hard contribution to the cross-section is described by and the azimuthal angles (φ r , φ p ) of vectors p T , r T are measured with respect to the lepton plane. However, our preferred frame is the one where the φ angle is measured in the plane defined by p T and q T , the sum of the lepton momenta, and φ r is the azimuthal angle between r T and p T . In this frame and integrating over the angle φ we are left with:

JHEP03(2022)047
with the factor 2π coming from φ integration. The LO expressions are obtained by separating the unpolarized and linearly polarized gluon contributions and Fourier transforming. The unpolarized partonic part has a similar form also for quarks, so that we find The same for the linearly polarized gluons gives and both functions can be related through eq. (2.20) in [3]. We obtain the σ (g,f )(U,L) 0 prefactors from the structure functions given in eqs. (3.3, 3.5) in [54] and we list them in appendix A.

Angle integrated and azimuthally modulated cross-section
The scalar cross-section that we finally consider in the phenomenological studies is obtained by integrating over the φ r angle where dΠ = dxdη 1 dη 2 dp T . Because the factorized cross-section is always expressed in position space one can write (here J 0,2 are Bessel functions) for the dijet case and dσ U,L = dσ U,L (γ * g) for the heavy hadron pair case.
In our phenomenological analysis we consider also the azimuthal angle average

JHEP03(2022)047
The denominator is what we have discussed in the previous section, so now we will only focus on the numerator, Eqs. (3.7)-(3.9) show the relation among the final cross-section in momentum space and the factorized cross-sections in position space, where we have distinguished between the unpolarized cross-sections generated by gluons and quarks, dσ U = dσ U g + dσ U f and the contribution from linearly polarized gluons dσ L .

Evolution kernels with angular dependent anomalous dimensions
The anomalous dimensions of soft, collinear-soft and heavy-jets matrix elements are angular dependent and complex valued, and both these features are not common in literature. The angular dependence is parameterized here by the angle φ b . In order to discuss the issue we show in section 4.1 the dijet soft function as obtained in [28] and then we extend the conclusions to the other functions that have a similar angular dependent structure. The angular dependence of the anomalous dimension is strictly correlated with the imaginary parts of the cross-section in position space. The treatment of the angular dependence in the evolution has been discussed also in [55][56][57] which propose an approximate treatment. In section 4.2 we perform an original analysis discussing in detail how passing from position space to momentum space in the cross-section allows to obtain a real valued cross-section, including resummation. The proof of this statement is here provided at one-loop, and the same mechanism can be conjectured to work at all orders in perturbation theory.

Dijet soft function and angle dependent anomalous dimensions
The dijet soft function presents several properties that result particularly interesting. The renormalized soft function is obtained from its bare expression where the factor Z S γi was calculated in [28]. For our purposes it is sufficient to report expression of the soft function in the MS scheme for the (γ * g)-channel

JHEP03(2022)047
and for the (γ * f )-channel . The expressions above depend explicitly on the angle φ b and we have usedb = b/|b|. The rapidity scale ζ 2 is responsible for the rapidity evolution of this factor and it is related to the TMD rapidity scale ζ 1 by the consistency constraint The values of ζ 1 and ζ 2 that minimize hard logs are Despite the fact that the scale ζ 2 is dimensionless there are some formal similarities in the evolution of this soft function and the TMD. The dijet soft function double scale evolution is dictated by where i is gluon or quark. The rapidity evolution kernel D i is the same as in the TMD case, while for the rest we need a special treatment.

Treatment of angular dependent anomalous dimensions and resummation
The dijet soft function, the collinear-soft function in dijets and the heavy meson jet function in hadron pair production have, as a common feature, an anomalous dimension that is φ b -dependent and that can include some imaginary parts. The general structure of the anomalous dimension for these cases is where c i and c i are the color coefficients that multiply the corresponding part of the cusp anomalous dimension. The construction of the evolution kernel involves an integration over the angle φ b and after integration all imaginary parts cancel consistently. Note that this does not mean that one can ignore the imaginary components of the anomalous dimensions or fixed order functions as contributions of imaginary terms yield real contributions to the final result.
In order to show the cancellation of the imaginary part we separate the angular dependent part of the evolution kernel from the rest. For instance for the soft function we have The evolution factor is so splitted in an angular dependent part and the rest. The splitting is clearly not unique, however this should not affect the final result once the angular integration is performed. We have A similar splitting is done for all other angular dependent evolution kernels, for which we have used 14) ) where the ± sign refers to quark and anti-quark jet respectively. The angular independent part of the anomalous dimension is simply obtained from the complete expression of the anomalous dimension

JHEP03(2022)047
with i = S γg , S γf , C g , C f , J for each channel. The complete list of one loop anomalous dimensions can be found in appendix B. In order to describe the implications of the angle integration we consider here the dijet production case, being the HHP one very similar. The separation of the evolution kernels of all functions in an angular dependent part and the rest allows to write the resummed cross-section in position space as where we have omitted global scale independent hard factors and non-perturbative contributions that we assume independent of the angle φ b . In this equation we have combined the evolution kernel of the functions that appear in the cross-section in which are independent of the factorization scale µ because of eq. (4.10), and an angle independent part R({µ k } → µ). The perturbative parts of all functions that appear in the cross-sections at one loop are collected in the factor 1 + k∈{H,F,J,S,C} a s (µ k )f [1] k . In eq. (4.20) the imaginary part of the cross-section is proportional to the odd function Θ(φ b ). We expect so that this part cancel at all orders in perturbation theory when the φ b the Fourier transform is performed like in eq. (3.7)-(3.9). In the next subsection we show that this is explicitly the case at one loop.
In order to organize the discussion we firstly consider the case of dσ U and then deduce the necessary modifications to get the Fourier transform of dσ L . The angular integration of the cross-section at one loop can always be expressed in terms of the following basic integrals:

JHEP03(2022)047
The results on the r.h.s. of eq. (4.23) are defined only for A > −1/2. In order to satisfy this condition for every value of b the scales µ S and µ C must be chosen appropriately, which at the current level of precision can be done. Next we show the contribution of each term of f [1] k to the cross-section: • Constant terms: in this case we consider the leading order term in eq. (4.20) together with all terms in the functions f [1] k that are cos φ b independent. The integral we need to perform is, and the imaginary part cancels because of parity.
• Single logarithmic terms: in this case we consider terms proportional to ln(−i cos φ b ) that appear in the soft and collinear-soft terms f [1] S and f [1] C respectively. The relevant integral that we need to perform is Expanding the logarithm as follows we can split the main integral into one part proportional to cos(Bπ) and a second part proportional to sin(Bπ). In terms of the integrals in eq.  • Logarithmic A b terms: in the soft function we observe the presence of terms that involve the following combination, Particularly this term appears in logarithms and di-logarithms. Here we address the double logarithmic terms, i.e., Integrating over φ b the above we obtain three contributions: a constant term, a single-log and a double-log. Thus we can immediately write,

JHEP03(2022)047
• Double logarithmic terms: in this case we consider terms proportional to ln 2 (−i cos φ b ) that appear in the soft and collinear-soft functions. The relevant integral that we need to perform is we can then expand the logarithm as follows then we can expand as before in terms proportional to sin(Bπ) and terms proportional to cos(Bπ) In addition we have poly-logarithmic terms which are proportional to Li 2 (1 − A b ), and where we have dropped the terms proportional to sin(Bπ) since they vanish after integration from parity. We can then use properties of the poly-logarithm in order to express this integral in the following form, We can then expand the poly-logarithms and the ln(1 −c cos 2 φ b ) in the region 0 < c ≤ 1. This allows us to perform the integral over φ b and obtain the following, Although this expression is exact it includes a sum extending to infinity. In numerical applications one can truncate this sum at the desired value to achieve a certain numerical precision. This concludes the discussion of all possible cases that appear in the Fourier transform of the unpolarized cross-section.

JHEP03(2022)047
The part of the cross-section relative to linearly polarized gluons can be treated similarly. In this case we need to incorporate an additional cos 2φ b term in the integrals but this is the only change since the soft and collinear-soft functions that appear in the two contributions are the same as for the dijet case. Using the trigonometric identity cos 2φ b = 2 cos 2 φ b − 1 the integrals for this case can be deduced from the discussion of the unpolarized cross-section by the replacement In all these cases one obtains a cancellation of the imaginary part of the cross-section.
Treating perturbatively the angular integration as discussed in this section leads to write eq. (3.7) for the dijet case as where R f,g are products of evolution kernels to be described in the next section, andσ U,L f,g are the result of φ b angular integration and can be written aŝ

JHEP03(2022)047
The functions C and S in eq. (4.44)- (4.46) are the result of the φ b integration in collinearsoft and dijet soft functions. For the heavy meson case we have just contributions from gluon scattering, (4.47) and we have to change J f,f → H + and C f → J Q→H in eq. (4.41)-(4.43). In the case of angular modulation the cross-sections can also be written as in eq. (4.40)-(4.47), with the correct values of the functions I const. , C and S. The non-perturbative effects are in all cases encoded in the evolution kernels, TMD and jet functions. In the next section we describe how the evolution kernels are defined.

Evolution kernels and scale choices
The evolution kernels appearing in eq. (4.40)-(4.47) are where R J f,g is a jet function kernel, R C f,g is the one of collinear-soft functions, R q,g F the one of TMD and finally R qg S is the one of the dijet soft function. In the heavy quark case the evolution kernels are parameterized like in eq. (5.1) with the usual changes J f → H + and C f → J Q→H . The kernels for single-scale evolution have a standard form and a review up to NLL is given in [58], where

JHEP03(2022)047
with r = α s (p T )/α s (µ i ) and Initial scales µ i choice is given in section 6. The TMD kernel is considered here in the ζ-prescription described in [40] and implemented in the code Artemide [37, 38] that we use, (5.7) In the next paragraph we define a ζ-prescription also for the dijet evolution kernel R g S (µ 0 , ζ 2,0 ) → (p T , 1), which is the only missing part.

ζ-prescription for dijet evolution kernel
The angular independent kernel of the dijet soft function is obtained as a solution of a coupled system of differential equations, reported in eq. (4.13), that are formally very similar to the TMD ones [59,60]. The anomalous dimensions are given bȳ (5.10) and δγ S are the non-cusp SF anomalous dimension, which is known up to three-loops for the gluon-channel and up to one-loop for the quark-channel and are reported in appendix. The anomalous dimension and the rapidity anomalous dimension (RAD) in eq. (4.6), (4.7) satisfy also The evolution for the SF takes the general form

JHEP03(2022)047
ensures that the evolution kernel is path only independent when one knows the complete perturbative expansion of the anomalous dimensions. Since this is not the case the path independence is broken. In order to partially restore the path independence we proceed as in [40] defining a ζ-prescription also for the dijet soft function evolution kernel. The ζ-prescription provides a way to choose the initial scale ζ i of the evolution kernel as a function of µ and b so that the SF does not depend on the initial scale µ i . This is done by taking the integration path through a null-evolution line in the {µ, ζ}-plane and then taking a fixed-µ evolution.
To find the null-evolution line we interpret the pair of differential equations (5.11) as a two-dimensional gradient equation ∇F = EF , where E = (γ S (µ, ζ), −D S (µ, b)). The null-evolution line is then an equipotential line of the field E. In particular, there is a special null-evolution line that passes through the saddle-point {µ saddle , ζ saddle } of the evolution field. We find that the saddle point is exactly µ saddle = µ 0 and ζ γi saddle = ζ γi 0 . If we parameterize the null-evolution line as {µ, ζ µ (b)}, the value of ζ µ is given by which is solved perturbatively order by order in α s . The perturbative solution takes the form

JHEP03(2022)047
and we are using the following notation with C i = C F , C A for quark and gluon channel respectively. Notice that v 0 vanishes as it is proportional to the LO non-cusp AD, which is zero for the SF. The non-cusp AD is not known beyond LO for the quark-channel. The RAD is a function of b and therefore has important non-perturbative corrections in the large-b region. These corrections can be implemented as a model. The way to proceed is to solve (5.13) for a generic non-perturbative RAD. The equation is solvable but it is difficult to obtain the cancellation of perturbative logarithms in the small-b region. Following [40] we use the perturbative solution for the small-b region and move to the exact (generic RAD) solution for large-b: with B NP being the b value where non-perturbative (NP) effects become important (∼ 2.5 GeV −1 ). We have already discussed the perturbative solution to eq. (5.13). For the exact solution we find where g γi (a s , D S ) = 1 a s Γ 0 2β 2 0 ∞ n=0 a n s g γi n (D S ), and p = 2β 0 D S /Γ 0 . Finally, the evolution kernel that provides the evolution from the null-evolution line and that passes through the saddle-point to the final ζ point is given by

JHEP03(2022)047
and if we consider the evolution from an arbitrary initial scale we take . (5.30) with i = q, g. This discussion concludes the analysis of all terms that appear in the factorization theorem and the scale prescription. We are now ready for the implementation in the code Artemide [37, 38].

Dijet and heavy hadron pair (HHP) production at EIC
In order to test the phenomenology developed in the previous sections we consider the case of the EIC. In [28] we already studied the coverage of the EIC and we concluded that the most favourable case is given for a value of mass energy for dijet production around √ s = 140 GeV and central rapidity, η 1 = η 2 = 0. Typical values for jet radii and momenta at EIC are respectively R ∼ 0.7 and p T ∼ Q/2 ∼ 20 GeV. In order to simplify the discussion we show plots integrated over Bjorken variable x (the longitudinal fraction of momentum ξ that enters in the TMDPDFs is ξ ∼ 2x) in the allowed kinematic intervals. For the case of central rapidity we have x ∈ (0.0859, 0.5). The cross-sections that we plot are xmax x min dx dσ dxdη 1 dη 2 dp T dr T η 1 , η 2 , p T (6.1) and its value is presented as a function of the small transverse momentum r T . The crosssections and the error bands are obtained by using and preparing specific moduli for the code Artemide [37, 38]. In particular we use the TMD and the TMD evolution kernels already coded in Artemide, that come from the fit [39], while the new functions studied in this work are included in this code for the first time. The gluon TMD is not fitted yet, however in Artemide there is a parameterization for it. The code takes into account that the contribution of linearly polarized gluons is highly suppressed because in the small-b regime the matching of the linearly polarized gluon TMD onto the gluon PDF starts at order α 1 s and not at order α 0 s like other distributions. In ref. [5] the cross section obtained in this way agrees with Pythya 8 and current experimental results for the Higgs transverse momentum spectrum, which are however not very precise. The non-perturbative effects are expected to be important in the high-b region and they should not alter the small-b behavior of this distribution. Notice also that the non-perturbative effects play a role to control the behavior of the distribution around the Landau pole at large-b, which means a further suppression effect at large-b (as we also observe in the case of unpolarized distributions). Summing up, given the current perturbative and non-perturbative knowledge of TMDs, at this stage we prefer not to push for a hypothetical non-perturbative enhancement of the contribution of linearly polarized gluon TMD.
The factorization that we propose in general needs information of the non-perturbative effects in several functions. For the dijet case we have where the functions with suffix pert refer to their perturbative part in the MS scheme which is currently known at one loop. Similarly for the HHP case we need The non-perturbative effects are parameterized as The values of B i NP define thee non-perturbative model and we have tested several combinations as shown in figure 2. Higher values of B i NP are more sensitive to the perturbative series in the low transverse momentum spectrum, and in general provide higher values of the observables. In unpolarized TMD cases we have usually that typical values of B i NP are around 1-3 GeV −1 so we have found reasonable to fix their values as in table 1.
The factorization scales µ C for dijet and µ J for HHP are chosen to minimize perturbative logarithms and to not hit the Landau pole of the strong coupling constant, This scale choice deserves some comments. In the dijet case the scale choice does not include the dependence on the jet radius R. Similarly, the mass of the ratio m Q /p T does not enter the collinear-soft function and heavy meson jet. In all cases, this means that there is not a complete cancellation of the logarithms of these functions. The reason is that the φ b integration imposes some constraints on the choice of scales. In fact, the function A({µ i }) defined in eq. (4.21), which depends on the initial scale choice for the soft function, collinear-soft function and heavy meson jet, needs to be A > −1/2 in order to have a well defined angular integration. Because of this constraint some scale choices which could be considered like for instance can not be used. As a result in our approach we only partially resum the logs in the collinear-soft function and the heavy meson jet in order to maintain the structure of ζprescription and double scale evolution in the soft function that is described in section 5. This leads to the initial scales in eq. (6.6), (6.7).

JHEP03(2022)047
Finally for the dijet soft function we use the b * -prescription in the same way as for the TMDPDF: (6.9) Concerning the theoretical errors, the scale variations in collinear-soft and heavy meson jet function are the main source. This is due to the non-cancellation of logs in the functions by the choice of the initial scales. The choice of the values b max for collinear-soft function and heavy meson jet account for the convergence of our perturbative result. A more consistent way to treat the resummation of these scales is left for a future work, involving the refactorization of these functions.
For functions that do not depend on b the initial scale choice does not require a prescription or NP-model and it is dictated by the cancellation of the logarithms. For the jet function and the H + matching coefficient we have We use a NP-model for the rapidity anomalous dimension that enters the exact solution for the null-evolution ζ µ line as it is explained in section 5. In particular, we use the same model that has been used for TMDPDF in [39] D NP F,S = c 0 bb * , c 0 = 2.5 · 10 −2 . (6.11) This model dictates how the rapidity anomalous dimension behaves in the large-b region and is used for both dijet soft function and TMDPDF when performing double scale evolution. While a color re-scaling of the non-perturbative models for gluon TMDPDF and gluon channel soft function with respect to their quark analogues is possible, we observe that this change does not have a significant impact on the cross-section and, therefore, we choose to keep the same model for both quarks and gluons.

Results
In this section we show our results for the differential cross-section for both dijet and heavy hadron pair production processes. Differential cross-sections are shown with error bands coming from scale variation of the different final and initial scales of the functions appearing in our factorization formulas. Scale variation bands are obtained by changing the considered scale by a factor of 2 up and down relative to its central value.

Results for dijet production
In figure 3 we show the impact of the change of jet radius, jet transverse momentum (hard scale) and jet pseudorapidity over total dijet cross-section. We show that for the variation of the jet radius we see a change of around 20% on the cross-section from the central value when taking the jet radius to be ±0.2 from R = 0.7. For p T there is a variation of an order of magnitude in the total cross-section when taking ±5 GeV from 20 GeV. This corresponds to Q = 30, 40, 50 GeV respectively. Finally, for pseudorapidity variation we obtain an order of magnitude difference above and below when compared to the central JHEP03(2022)047 Figure 2. Impact of B NP variation over dijets and heavy meson total cross-section. Legend correspond to (B S NP , B C NP ) and (B S NP , B J NP ) for dijet and HHP production respectively Figure 3. Impact of the variation of the jet radius (R), hard scale jet transverse momentum (p T ) and jet pseudorapidity (η i ) for dijet production. For pseudorapidity variation legend is shown referring to (η 1 , η 2 ) pair, dashed and dotted lines correspond to negative and positive rapidity respectively. rapidity case. Positive rapidities (η i = 0.5) correspond to Q 53 GeV while negative rapidities (η i = −0.5) correspond to Q 32 GeV, so both p T and η plots are consistent. Notice that total dijet cross-section is not symmetrical for both jet rapidities as for quark channel we have both a quark and gluon jet in the final state. Every other plot is obtained taking R = 0.7, p T = 20 GeV and η i = 0.
In figure 4 the result for the cross-section including quark and gluon channels is shown. We consider the contribution of linearly polarized gluons in a separate panel to show that their contribution is completely negligible, being a factor 10 3 -10 4 smaller. This leads to the conclusion that the contribution from the linearly polarized gluons can be neglected when considering the unpolarized cross-section.
The angular modulation asymmetry is shown in figure 6, being around 5%.

Results for heavy hadron production
The analysis for HHP has followed similar steps of the dijet case when possible. The differential cross-section including all channels is plotted in figure 5.  Figure 5. Cross-sections for HHP production at EIC with error-bands coming from scale dependence in hard factor (Hard), heavy meson jet function (bHQET), heavy meson jet function matching coefficient (Hard+) and Wilson coefficients (OPE). The rows correspond to contributions from total cross-section (top) and linearly polarized gluons (bottom). √ s = 140 GeV, p T = 20 GeV, the contribution of linearly polarized gluons show also in this case that they are completely negligible being suppressed by a factor 10 2 -10 3 . The angular modulation asymmetry is shown in figure 7, being around 5%. Figure 6. Angular modulation contribution for dijet production at EIC with error-bands coming from scale dependence in collinear-soft factor (CSF), hard factor (Hard), jet distributions (Jet) and Wilson coefficients (OPE). √ s = 140 GeV, R = 0.7, p T = 20 GeV, η 1 = η 2 = 0. Figure 7. Angular modulation contribution for HHP production at EIC with error-bands coming from scale dependence in hard factor (Hard), heavy meson jet function (bHQET), heavy meson jet function matching coefficient (Hard+) and Wilson coefficients (OPE). √ s = 140 GeV, p T = 20 GeV, η 1 = η 2 = 0.

Conclusions
Dijet and HHP in SIDIS experiments present an opportunity to study the gluon TMD. In this work we have considered the case of the Electron Ion Collider (EIC), as an example. The processes have been proven to factorize consistently and this result has been checked at least at one loop [28]. Nevertheless the evolution of the functions that appear in the factorization theorem is non-trivial and we propose an original solution, which is generic and independent of the resummation framework. We also note that it is consistent with the ζ-prescription of TMD [40] which we implement in this work. The used prescription allows to separate the evolution kernels from other scale independent factors in the cross-section, so that in our final computations we can use some results already coded in the literature. This is the case for the TMD and their respective evolution kernels extracted from DY and SIDIS data and presented in the code Artemide [37,38].
The phenomenological analysis that we have performed has revealed several issues that need further study in the future. Estimating the errors due to scale dependence, we have found that several functions need a higher loop calculation to achieve sufficient precision for the low energy jets will be available at EIC. This is the case for instance of the collinearsoft function that appears in the dijet process and more urgently on the heavy hadron jet function. The fact that the perturbative convergence of these function is limited to small JHEP03(2022)047 values of the b may lead to consider also a re-factorization of these functions, such that the small-b effects are separately resummed. This possibility can eventually be considered in future works.
In all cross-sections we have found a contribution of unpolarized and linearly polarized gluon TMD. For both of these distributions we have used their re-factorization in coefficient functions and collinear PDF studied at higher loops in the literature [5,42,44,46]. In accordance to this well tested procedure, the linearly polarized gluon contribution results to be particularly suppressed in all considered cases, because its matching to collinear PDF starts at order α 1 s , instead of α 0 s like the unpolarized distributions. The effect of this suppression is particularly evident in the estimate of the angular modulation of cos 2φ r asymmetry that is here estimated to be around 5%. The study of next-to-leading power effects is beyond the purpose of the present work, so that further study is necessary to confirm a value of this asymmetry.
A source of uncertainty in our prediction comes from the usage of models for many functions that have not been yet compared against data. In this case we have studied several possibilities with simple Gaussian models, assigning values to the non-perturbative parameters according to an educated-guess. The models do not alter the overall-conclusions about scale choices or precision, but can have some effect on the shape of the curves that we have computed. Only an strict comparison with data or eventual lattice calculations can finally resolve this issue.
As a result of this study we can see that the extraction of gluon TMDs from dijet and HHP processes at the EIC is conditioned yet by the possible theoretical and experimental precision. In particular, the linearly polarized gluon TMD appears generally too suppressed and hardly accessible if one uses the usual matching of TMDs onto their collinear counterpart distributions. Nevertheless, the discussed theoretical issues can potentially be solved or improved in future studies.
where D is defined as (A.12)

B Anomalous dimensions
For the two channels in the dijet process the relevant anomalous dimensions up to oneloop are, γ [1] Hγg = 4 C F ln ŝ 2 µ 4 − 2γ q + C A ln tû sµ 2 , γ [1] H γf = 4 C F ln û 2 µ 4 − 2γ q + C A ln ŝt u µ 2 , γ [1] Sγg = 4 − C A ln ζ 2 + 2C F ln(Bµ 2 e 2γ E ) − lnŝ + ln p 2 T + ln(4c 2 b ) , γ [1] S γf = 4 (C F + C A ) ln(Bµ 2 e 2γ E ) − lnŝ + ln p 2 T + ln(4c 2 b ) The imaginary component in the soft and collinear-soft anomalous dimension is denoted by These anomalous dimensions, except the soft function which we calculated here, can be found in [47][48][49][50][51]61]. We also 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.
For the heavy hadron pair case the one-loop hard function, H + is, and the corresponding anomalous dimension is The heavy jet functions anomalous dimension is Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.