Higher-twist corrections to gluon TMD factorization

We calculate power corrections to TMD factorization for particle production by gluon-gluon fusion in hadron-hadron collisions.


Introduction
Particle production in hadron-hadron scattering with transverse momentum of produced particle much smaller than its invariant mass is described in the framework of TMD factorization [1][2][3][4][5]. The typical example is the Higgs production at LHC through gluon-gluon fusion. Factorization formula for particle production in hadron-hadron scattering looks like [1,6] dσ where η is the rapidity, D f /A (x, z ⊥ , η) is the TMD density of a parton f in hadron A, and σ(f f → H) is the cross section of production of particle H of invariant mass m 2 H = Q 2 in the scattering of two partons. (For simplicity, we consider the scattering of unpolarized hadrons.) In this paper we calculate the first power corrections ∼ q 2 ⊥ Q 2 in a sense that we represent them as a TMD-like matrix elements of higher-twist operators. It should be noted that our method works for arbitrary relation between s and Q 2 and between q 2 ⊥ and hadron mass m 2 (provided that pQCD is applicable), but in this paper we only present the result for the physically interesting region s Q 2 q 2 ⊥ m 2 . To obtain formula (1.1) with first corrections we use factorization in rapidity [7]. We denote quarks and gluons with rapidity close to the rapidity of the projectile and target protons as A-fields and B-fields, respectively. We call the remaining fields in the central region of rapidity by the name C-fields and integrate over them in the corresponding functional integral. At this step, we get the effective action depending on A and B fields. The subsequent integration over A fields gives matrix elements of some TMD-like operators switched between projectile proton states and integration over B fields will give matrix elements between target states. 1 The paper is organized as follows. In Sect. 2 we derive the TMD factorization from the double functional integral for the cross section of particle production. In Sect 2, which is central to our approach, we explain the method of calculation of higher-twist power corrections based on a solution of classical Yang-Mills equations. In Sect. 4 we find the leading higher-twist correction to particle production in the region s Q 2 q 2 ⊥ . Finally, in Sect. 5 we compare our calculations in the small-x limit to the classical field resulting from the scattering of two shock waves. The Appendices contain proofs of some necessary technical statements.

TMD factorization from functional integral
We consider production of an (imaginary) scalar particle Φ in proton-proton scattering. This particle is connected to gluons by the vertex This is a m H mt 1 approximation [12,13] for Higgs production via gluon fusion at LHC with g H = 1 48π 2 v 1 + 11 4π α s + ...
where α s = g 2 4π as usual. 2 The differential cross section of Φ production has the form where we defined the "hadronic tensor" W (p A , p B , q) as It should be noted that due to the kinematics Q 2 Q 2 ⊥ , m 2 we will not need the explicit form of the high-energy effective action which is much sought after in the small-x physics but not known up to now except a couple of first perturbative terms [7][8][9][10][11]. 2 For finite mt the constant gH should be multiplied by  As usual, X denotes the sum over full set of "out" states. It can be represented by double functional integral × Ψ * p B ( Ã (t i ),ψ(t i ))e −iS QCD (Ã,ψ) e iS QCD (A,ψ)F 2 (x)F 2 (0)Ψ p A ( A(t i ), ψ(t i ))Ψ p B ( A(t i ), ψ(t i )) Here the fields A, ψ correspond to the amplitude X|F 2 (0)|p A , p B , fieldsÃ,ψ correspond to complex conjugate amplitude p A , p B |F 2 (x)|X and Ψ p ( A(t i ), ψ(t i )) denote the proton wave function at the initial time t i . The boundary conditionsÃ(t f ) = A(t f ) andψ(t f ) = ψ(t f ) reflect the sum over all states X, cf. Refs. [15], [16], [17]. We use Sudakov variables p = αp 1 + βp 2 + p ⊥ and the notations x • ≡ x µ p µ 1 and x * ≡ x µ p µ 2 for the dimensionless light-cone coordinates (x * = s 2 x + and x • = s 2 x − ). Our metric is g µν = (1, −1, −1, −1) so that p · q = (α p β q + α q β p ) s 2 − (p, q) ⊥ where (p, q) ⊥ ≡ −p i q i . Throughout the paper, the sum over the Latin indices i, j... runs over the two transverse components while the sum over Greek indices runs over the four components as usual.
To derive the factorization formula, we separate the (quark and gluon) fields in the functional integral (2.4) into three sectors: "projectile" fields A µ , ψ a with |β| < σ a , " target" fields with |α| < σ b and "central rapidity" fields C µ , ψ with |α| > σ b and |β| > σ a : 3 Our goal is to integrate over central fields and get the amplitude in the factorized form, as a (sum of) products of functional integrals over A fields representing projectile matrix elements (TMDs) and functional integrals over B fields representing target matrix elements. In the spirit of background-field method, we "freeze" projectile and target fields (and denote them theĀ,ξ a , ξ a andB,ξ b , ξ b respectively) and get a sum of diagrams in these external fields. Since |β| < σ a in the projectile fields and |α| < σ b in the target fields, at the tree-level one can set with power accuracy β = 0 for the projectile fields and α = 0 for the target fields -the corrections will be O m 2 σas and O m 2 σ b s . Beyond the tree level, one should expect that the integration over C fields will produce the logarithms of the cutoffs σ a and σ b which will cancel with the corresponding logs in gluon TMDs of the projectile and the target.
As usual, diagrams disconnected from the vertices F 2 (x) and F 2 (0) ("vacuum bubbles" in external fields) exponentiate so the result has the schematic form where O µν (q, x; A, ψ A ; B, ψ B ) is a sum of diagrams connected toF 2 (x)F 2 (0). Since rapidities of central fields and A, B fields are very different, one should expect the result of integration over C-fields to be represented in terms of Wilson-line operators constructed form A and B fields. The effective action has the form where Wilson lines U are made from projectile fields and Wilson lines V from target fields and similarly forŨ andṼ in the left sector. The explicit form of "Lipatov vertices" L i (U, V ) is presented in [20]. Unfortunately, the effective action beyond the first two terms in (2.7) is unknown, but we will demonstrate below that for our purposes we do not need the explicit form of the effective action.
After integration over C fields the amplitude (2.4) can be rewritten as Note that due to boundary conditions at t f in the above integral, the functional integral over C fields in Eq. (2.6) should be done in the background of the A and B fields satisfying Our approximation at the tree level is that β = 0 for A,Ã fields and α = 0 for B,B fields which corresponds to . Now comes the important point: because of boundary conditions (2.9), for the purpose of calculating the integral (2.6) over central fields one can set Indeed, because A, ψ andÃ,ψ do not depend on x * , if they coincide at x * = ∞ they should coincide everywhere. Similarly, if B, ψ b andB,ψ b do not depend on x • , if they coincide at x • = ∞ they should be equal. It should be emphasized that the boundary conditions (2.9) mean the summation over all intermediate states in corresponding projectile and target matrix elements in the functional integrals over projectile and target fields. Without the sum over all intermediate states the conditions (2.10) are no longer true. For example, if we would like to measure another particle or jet in the fragmentation region of the projectile, the second condition in Eq. (2.10) breaks down.
Next important observation is that due to Eqs. (2.10) the effective action (2.7) vanishes for background fields satisfying conditions (2.9). For the first two terms displayed in (2.7) it is evident, but it is easy to see that the effective action in the background fields satisfying (2.10) should vanish due to unitarity. Indeed, let us consider the functional integral (2.4) without sourcesF 2 (x)F 2 (0). It describes the matrix element (2.11) without Φ production, that is (modulo appropriate normalization of |p A and |p B states). If we perform the same decomposition into A, B, and C fields as in Eq. (2.4) we will see integral (2.8) without O µν (q, x, y; A, ψ a ,Ã,ψ a ; B, ψ b ,B,ψ b ) which can be represented as which means that the effective action should vanish for the Wilson-line operators constructed from the fields satisfying Eqs. (2.10). Summarizing, we see that at the tree level in our approximation . It is known that in the tree approximation the double functional integral (2.13) is given by a set of retarded Green functions in the background fields [21][22][23] (see also Appendix A for the proof). Since the double functional integral (2.13) is given by a set of retarded Green functions (in the background field A + B), the calculation of tree-level contributions to, say, F 2 (x) in the r.h.s. of Eq. (2.13) is equivalent to solving YM equation for A µ (x) (and ψ(x)) with boundary conditions that the solution has the same asymptotics at t → −∞ as the superposition of incoming projectile and target background fields. The hadronic tensor (2.8) can now be represented as whereÔ(q, x;Â,ψ a ;B,ψ b ) should be expanded in a series inÂ,ψ a ;B,ψ b operators and evaluated between the corresponding (projectile or target) states: if (where c µν m,n are coefficients and Φ can be any of A µ , ψ orψ) then 4 As we will demonstrate below, the relevant operators are quark and gluon fields with Wilsonline type gauge links collinear to either p 2 for A fields or p 1 for B fields.

Power counting for background fields
As we discussed in previous Section, to get the hadronic tensor in the form (2.14) we need to calculate the functional integral (2.13) in the background of the fields (2.10). To understand the relative strength of Lorentz components of these fields, let us compare the typical term in the leading contribution to W and some typical higher-twist terms. As we mentioned, we consider W (p A , p B , q) in the region where s, Q 2 Q 2 ⊥ , m 2 while the relation between Q 2 ⊥ and m 2 and between Q 2 and s may be arbitrary. So, for the purpose of counting of powers of s, we will not distinguish between s and Q 2 (although at the final step we will be able to tell the difference since our final expressions for higher-twist corrections will have either s or Q 2 in denominators). Similarly, for the purpose of power counting we will not distinguish between m and Q ⊥ and will introduce m ⊥ which may be of order of m or Q ⊥ depending on matrix element.
The estimate of the leading-twist matrix element between projectile states is Our logic here is the following: to get the expression forÔ in Eq. (here we assume normalization p A |p A = 1 for simplicity). The typical higher-twist correction is proportional to (see e.g. Eq. (4.4)) we see that an extraF µi in the matrix element between projectile states brings p 1µ m ⊥ which means thatÛ * i ∼ sm ⊥ .
Next, some of the higher-twist matrix elements have an extra U kl like Since we consider only unpolarized projectile and target hadrons and, comparing this to Eq. (3.3), we see that an extraF kl can bring an extra m 2 ⊥ . Combining this with an estimate U * i ∼ sm ⊥ we see that the typical fieldĀ * is of order s whilē A i ∼ m ⊥ . Similarly, for the target fields we getB • ∼ s,B i ∼ m ⊥ .
Some of the power corrections involve matrix elements like An extra field strength operatorF µν between the projectile states can bring Similarly, for the target we get B * ∼ m 2 ⊥ . Summarizing, the relative strength of the background gluon fields in projectile and The denominator pA · p2 is due to the fact that p2 enters only through the direction of Wilson line and therefore the matrix element should not change under rescaling p2 → λp2 To finish power counting, we need also the relative strength of quark background fields ψ a and ψ b . From classical equations for projectile and target Thus, to find TMD factorization at the tree level (with higher-twist corrections) we need to calculate the functional integral (2.4) in the background fields of the strength given by Eqs. (3.10) and (3.12).

Approximate solution of classical equations
As we discussed in Sect 2, the calculation of the functional integral (2.13) over C-fields in the tree approximation reduces to finding fields C µ and ψ c as solutions of Yang-Mills equations for the action As we discussed above (see also Appendix A) the solution of Eq. (3.13) which we need corresponds to the sum of set of diagrams in background fieldĀ +B with retarded Green functions (see Fig. 3). The retarded Green functions (in the background-Feynman gauge) are defined as and similarly for quarks. The solutions of Eqs. (3.13) in terms of retarded Green functions give fields C µ and ψ c that vanish at t → −∞. Thus, we are solving the usual classical YM equations following from C µ , ψ c t→−∞ → 0. These boundary conditions reflect the fact that at t → −∞ we have only incoming hadrons with "A" and "B" fields.
The solution of YM equations (3.16) in general case is yet unsolved problem, especially important for scattering of two heavy nuclei in semiclassical approximation. Fortunately, for our case of particle production with q ⊥ Q 1 we can construct the approximate solution of (3.16) as a series in this small parameter. However, before doing this, it is convenient to perform a gauge transformation so that the incoming projectile and target fields will no longer have large components ∼ s asĀ * andB • in Eq. (3.10). Let us perform the gauge transformation of Eq. (3.16) and initial conditions (3.17) with the gauge matrix Ω(x) such that The existence of such matrix is proved in Appendix B by explicit construction. After such gauge transformation, the YM equation of course stays the same but the initial conditions (3.17) turn to and Σ a , Σ b are defined as The initial conditions (3.19) look like the projectile fields in the light-like gauge p µ 2 A µ = 0 and target fields in the light-like gauge p µ 1 A µ = 0 so our construction of matrix Ω in a way proves that we can take the sum of projectile fields in one gauge and target fields in another gauge as a zero-order approximation for iterative solution of the YM equations. Note also that our power counting discussed in previous Section means that so we do not have large background fields ∼ s after this gauge transformation. Finally, the classical equations for projectile and target fields in this gauge read 6 : and similarly for V fields. We will solve Eqs. (3.16) iteratively, order by order in perturbation theory, starting from the zero-order approximation in the form of the sum of projectile and target fields and improving it by calculation of Feynman diagrams with retarded propagators in the background fields (3.24). The first step is the calculation of the linear term for the trial configuration (3.24). We rewrite field strength components as Here we consider only u, d, and s quarks which can be regarded as massless.
⊥ while all other components are not large. The linear term has the form for gluon fields (in the background-Feynman gauge) and for quarks where are operators in external zero-order fields (3.24). Hereafter we use Schwinger's notations for propagators in external fields normalized according to (x|F (p)|y) ≡ d −4 p e −ip(x−y) F (p) . Moreover, when it will not lead to a confusion, we will use short-hand notation . Next iterations will give us a set of tree-level Feynman diagrams in the background field U µ + V µ and Σ a + Σ b .
Let us consider the fields in the first order in perturbation theory: Here α, β, and p ⊥ are understood as differential operators α = i ∂ ∂x• , β = i ∂ ∂x * and p i = i ∂ ∂x i . Now comes the central point of our approach. Let us expand quark and gluon propagators in powers of background fields, then we get a set of diagrams shown in Fig. 3. The typical bare gluon propagator in Fig. 3 is Since we do not consider loops of C-fields in this paper, the transverse momenta in tree diagrams are determined by further integration over projectile ("A") and target ("B") fields in Eq. (2.8) which converge on either q ⊥ or m. On the other hand, the integrals over α converge on either α q or α ∼ 1 and similarly the characteristic β's are either β q or ∼ 1.
Since α q β q s = Q 2 Q 2 ⊥ , one can expand gluon and quark propagators in powers of The explicit form of operators 1 α+i , 1 β+i , and After the expansion (3.33), the dynamics in the transverse space effectively becomes trivial: all background fields stand either at x or at 0. (This validates the reasoning in the footnote on page 3).
One may wonder why we do not cut the integrals in Eq. (3.34) to |α| > σ b and |β| > σ a according to the definition of C fields in Sect. 2. 7 The reason is that in the diagrams like Fig. 3 with retarded propagators (3.34) one can shift the contour of integration over α and/or β to the complex plane away to avoid the region of small α or β. 8 Note that the background fields are also smaller than typical p 2 ∼ s. Indeed, from Eq.
The only exception is the fields V •i or U * i which are of order of sm ⊥ but we will see that effectively the expansion in powers of these fields is cut at the second term with our accuracy.

Twist expansion of classical gluon fields
Now we expand the classical gluon fields in powers of s . It is clear that for the leading higher-twist correction we need to take into account only the first two terms (3.28) of the perturbative expansion of classical field. The expansion (3.28) of gluon field A • takes 7 Such cutoffs for integrals over C fields are introduced explicitly in the framework of soft-collinear effecive theory, see the review [24]. 8 This may be wrong if there is pinching of poles in the integrals over α or β but we will see that in our integrals for the tree-level power corrections the pinching of poles never occurs. In the higher orders in perturbation theory the pinching does occur so one needs to formulate a subtraction program to avoid double counting.

the form
Similarly, from Eq. (3.28) one obtains The corresponding expansion of field strengths reads Power corrections to hadronic tensor are proportional to and the leading higher-twist correction is proportional to 4 Leading higher-twist correction at s Q 2 Q 2 ⊥ m 2 As we mentioned in the Introduction, our method is relevant for calculation of higher-twist corrections at any s, Q 2 Q 2 ⊥ , m 2 . However, the expressions become manageable in the physically interesting case s Q 2 Q 2 ⊥ m 2 which we consider in this Section. 9 We will demonstrate that the leading correction in this region comes from the following part of Eq. (3.41) The higher-twist correction coming from the second term in the r.h.s. will be ∼ βqs all of which are small (see the footnote 9). In this approximation we get where the first term is the leading order and the second is the higher-twist correction. Substituting our approximation (4.1) to Eq. (2.3) and promoting background fields to operators as discussed in Sect. 2 we get (note that α q β q s = Q 2 Q 2 ): where we used formula [25,26] Since an extra U * k (or V •k ) brings s x i We parametrize gluon TMD for unpolarized protons as (cf. Ref. [27]) where σ b is the cutoff in α integration in the target matrix elements, see the discussion in Ref. [18]. The normalization here is such that D g (β q , k 2 ⊥ ; σ b ) is an unintegrated gluon distribution: 10 To see this, we compared matrix elements of leading-twist operator pA|U mi * (x•, x ⊥ )U mj * (0)|pA and higher-twist operator pA|U a where D g (β q , µ 2 ) is the usual gluon parton density (this formula is correct in the leading log approximation, see the discussion in Ref. [18]). Next, the three-gluon matrix element in Eq. (4.4) for unpolarized hadrons can be parametrized as At large k 2 ⊥ gluon TMDs in the r.h.s. of Eq. (4.6) behave as D g (β q , k 2 ⊥ ) ∼ 1 It is well known that in our kinematic region s Q 2 Q 2 ⊥ gluon TMDs (4.6) possess Sudakov logs of the type Let us now demonstrate that the terms in (F 2 (x)) (0) (see Eq. (3.42)) which we neglected give small contributions. For example, consider the following contribution to F 2 (x)F 2 (0): The corresponding contribution to hadronic tensor W has the form Note that unlike Eq. (4.4), the factor in the denominator is α q s Q 2 so the contribution (4.11) is power suppressed in comparison to Eq. (4.4) in our kinematic region. 11 As a less trivial example, consider the following term in The corresponding contribution to hadronic tensor W reads where we used In both examples (4.11) and (4.13) the factor 1 αq comes from an extra integration over 14) The way to figure out such integrations is very simple: take α q → 0 and check if there is an infinite integration of the type Thus, to get the terms ∼ 1 Q 2 we need to find contributions which satisfy both of the above conditions which singles out the contribution (4.3).

Small-x limit and scattering of shock waves
Let us consider the hadronic tensor in the small-x limit s → ∞, Q 2 and q 2 ⊥ -fixed. At first, let us not impose the condition Q 2 q 2 ⊥ which means that the relation between x 2 and x 2 ⊥ is arbitrary (later we will see that Q 2 q 2 ⊥ corresponds to x 2 x 2 ⊥ ). The small-x limit may be obtained by rescaling s → λ 2 s ⇔ p 1 → λp 1 , p 2 → λp 2 . As discussed in Refs. [7,20,28], the only components of field strength surviving in this rescaling are U * i (x • , x ⊥ ) and V •i (x * , x ⊥ ). Moreover, if we study classical fields at longitudinal distances which does not scale with λ, we can replace the projectile and target fields by infinitely thin "shock waves" However, since we need to compare the classical fields in the small-x limit to our expressions (3.40) at small longitudinal distances, we will keep x * and x • dependence for a while. As described above, to find the classical fields we can start with the trial configuration with the linear term At longitudinal distances x • , x * ∼ 1 these expressions agree with Eq. (52) from Ref. [7] after the replacement (5.2). Now let us compare the fields (5.8) at small longitudinal distances to our approximate solution (3.40). Let us start with F ij (x) in the last line in Eq. (5.8). If (x − z) 2 is smaller than the characteristic transverse distances in the integral over z ⊥ one can replace On the other hand, the first line in Eq.
which agrees with Eq. (5.11) after integration over z ⊥ . Similarly, one can check the consistency of two expressions for F * i .

Conclusions and outlook
We have formulated the approach to TMD factorization based on the factorization in rapidity and found the leading higher-twist contribution to the production of a scalar particle (e.g. Higgs) by gluon-gluon fusion in the hadron-hadron scattering. Up to now our results are obtained in the tree-level approximation when the question of exact matching of cutoffs in rapidity does not arise. However, this question will become crucial starting from the first loop. In our previous papers we calculated the evolution of gluon TMD with respect to our rapidity cutoff so we need to match it to the coefficient functions in front of TMD operators. The work is in progress. Also, we obtained power corrections for particle production only in the case of gluongluon fusion. It would be interesting (and we plan) to find power corrections to Drell-Yan process. There is a statement that for semi-inclusive deep inelastic scattering (SIDIS) the leading-order TMDs have different directions of Wilson lines: one to +∞ and another to −∞. We think that the same directions of Wilson lines will be in the case of power corrections and we plan to study this question in forthcoming publications.

Appendix A
In this Section we prove that the field C µ created by a source J µ in the presence of external fieldsĀ µ andB µ is given by a set of Feynman diagrams with retarded Green functions (note that Eq. (7.1) implies that J µ ,Ā µ , andB µ are the same in the right and left part of the amplitude). Hereafter we use the notation µν ≡P 2 g µν + 2iḠ µν . First, we consider gluon propagators for the double functional integral over C fields in the background fileldsĀ =Ā,B =B and prove that Note that we define O in this Section as To prove Eq. (7.2), we write down µν = p 2 g µν + O µν , O µν ≡ {p ξ ,Ā ξ +B ξ } + (Ā +B) 2 g µν + 2iḠ µν (7.4) and expand in powers of O µν . In the trivial order Eqs. (7.2) immediately follow from the bare propagators for the double functional integral (7.3) 12 For simplicity, in this section we disregard quarks so in our case Jµ is Eq. (3.26) without quark terms.
In the first order in O µν we get Similarly, it is easy to see that In the second order in O µν so using the results (7.8) and (7.10) we get One can prove now Eq. (7.2) by induction using formulas Now we are in a position to prove Eq. (7.1). In the leading order in g it is trivial: using Eqs. (7.2) one immediately sees that In the first order in g (with one three-gluon vertex) we obtain C a µ (x) which is the desired result. Similarly, in the g 2 order one obtains after some algebra To find matrix Ω(x) satisfying Eqs. (3.18) we will solve the following auxiliary problem: we fix x ⊥ as a parameter and find the solution of Yang-Mills equations Since 2-dimensional gluodynamics is a trivial theory, the solution of the equation (8.1) will be a pure-gauge field A µ = Ωi∂ µ Ω † with Ω(x * , x • ) being the sought-for matrix satisfying Eqs. (3.18). Let us first demonstrate that the solution A µ (x * , x • ) of the YM equations (8.1) with boundary conditions (8.2) in two longitudinal dimensions is a pure gauge. To this end, we will construct A µ (x * , x • ) order by order in perturbation theory (see Fig. 3, but now in two dimensions) and prove that F a µν (A) = 0. We are looking for the solution of Eq. (8.1) in the form we get the equation The boundary conditions (8.2) in terms of C fields read It is convenient to rewrite the equation (8.5) in components as 2(P •P * ) abCb We will solve this equation by iterations inF * • and prove that F * • = 0 in all orders.
In the first order we get the equation The solution satisfying boundary conditions (8.6) has the form Using the explicit form of the propagators in externalĀ * andB • fields we getC (1) in the form where we used Eq. (8.12) to reduce the r.h.s. Again, from the explicit form of the propagators (8.10) we get • (x * , z • )][z • , x • ] A * (8.15) from which it is clear thatC (2) µ satisfy boundary conditions (8.6) (recall that we already proved thatC (1) µ satisfy Eq. (8.6)). Next, we usē (1)c • (8.16) to prove that F * • vanishes in the second order: In the third order we get  One can continue and prove by induction that F * • vanishes in an arbitrary order in G n * • and therefore the field A µ is a pure gauge Now we shall demonstrate that the matrix Ω satisfies our requirement (3.18). Since C * (x * → −∞, x • ) = 0 due to Eq. (8.2), we get Similarly, One can also construct the expansion of matrix Ω in powers ofĀ * andB • . For example, up to the fifth power of theĀ µ fields Ω(x * , x • ) (8.25) Now, for each x ⊥ we solve auxiliary 2-dimensional classical problem (8.1) and find Ω(x * , x • , x ⊥ ) satisfying the requirement (3.18).