Transverse Parton Distribution and Fragmentation Functions at NNLO: the   Quark Case

Ming-Xing Luo; Xing Wang; Xiaofeng Xu; Li Lin Yang; Tong-Zhi Yang and; Hua Xing Zhu

arXiv:1908.03831·hep-ph·January 8, 2020

Transverse Parton Distribution and Fragmentation Functions at NNLO: the Quark Case

Ming-Xing Luo, Xing Wang, Xiaofeng Xu, Li Lin Yang, Tong-Zhi Yang and, Hua Xing Zhu

PDF

TL;DR

This paper introduces an exponential regulator for calculating quark transverse momentum dependent functions at NNLO, simplifying the process and enabling higher-order precision in QCD calculations.

Contribution

The paper presents a novel exponential regulator that streamlines NNLO calculations of TMD functions and provides the two-loop quark jet function for N$^3$LL resummation.

Findings

01

Simplified NNLO calculations of TMD functions using the exponential regulator.

02

Obtained the two-loop quark jet function for Energy-Energy Correlator.

03

Enabled higher-order resummation in QCD analyses.

Abstract

We revisit the calculation of perturbative quark transverse momentum dependent parton distribution functions and fragmentation functions using the exponential regulator for rapidity divergences. We show that the exponential regulator provides a consistent framework for the calculation of various ingredients in transverse momentum dependent factorization. Compared to existing regulators in the literature, the exponential regulator has a couple of advantages which we explain in detail. As a result, the calculation is greatly simplified and we are able to obtain the next-to-next-to-leading order results up to $O (ϵ^{2})$ in dimensional regularization. These terms are necessary for a higher order calculation which is made possible with the simplification brought by the new regulator. As a by-product, we have obtained the two-loop quark jet function for the Energy-Energy…

Equations384

Q^{2} \equiv - q^{2}, x \equiv \frac{Q ^{2}}{2 P _{1} \cdot q}, z \equiv \frac{P _{1} \cdot P _{2}}{P _{1} \cdot q} .

Q^{2} \equiv - q^{2}, x \equiv \frac{Q ^{2}}{2 P _{1} \cdot q}, z \equiv \frac{P _{1} \cdot P _{2}}{P _{1} \cdot q} .

k^{μ} = \overset{n}{ˉ} \cdot k \frac{n ^{μ}}{2} + n \cdot k \frac{n ˉ ^{μ}}{2} + k_{⊥}^{μ} \equiv k_{+} \frac{n ^{μ}}{2} + k_{-} \frac{n ˉ ^{μ}}{2} + k_{⊥}^{μ} .

k^{μ} = \overset{n}{ˉ} \cdot k \frac{n ^{μ}}{2} + n \cdot k \frac{n ˉ ^{μ}}{2} + k_{⊥}^{μ} \equiv k_{+} \frac{n ^{μ}}{2} + k_{-} \frac{n ˉ ^{μ}}{2} + k_{⊥}^{μ} .

p \cdot k = \frac{p _{+} k _{-} + p _{-} k _{+}}{2} + p_{⊥} \cdot k_{⊥} .

p \cdot k = \frac{p _{+} k _{-} + p _{-} k _{+}}{2} + p_{⊥} \cdot k_{⊥} .

P_{1} = (P_{1 +}, 0, 0_{⊥}), P_{2} = (0, P_{2 -}, 0_{⊥}), q = (q_{+}, q_{-}, q_{⊥}),

P_{1} = (P_{1 +}, 0, 0_{⊥}), P_{2} = (0, P_{2 -}, 0_{⊥}), q = (q_{+}, q_{-}, q_{⊥}),

q_{-} = \frac{P _{2 -}}{z} = \frac{Q ^{2}}{x P _{1 +}}, q_{+} q_{-} + q_{⊥}^{2} = - Q^{2},

q_{-} = \frac{P _{2 -}}{z} = \frac{Q ^{2}}{x P _{1 +}}, q_{+} q_{-} + q_{⊥}^{2} = - Q^{2},

W^{μν} \equiv X \sum (2 π)^{4} δ^{(4)} (P_{1} + q - P_{2} - P_{X}) \times ⟨ N_{1} (P_{1}) ⟩ J^{μ} (0) N_{2} (P_{2}), X ⟨ N_{2} (P_{2}), X ⟩ J^{ν} (0) N_{1} (P_{1}) .

W^{μν} \equiv X \sum (2 π)^{4} δ^{(4)} (P_{1} + q - P_{2} - P_{X}) \times ⟨ N_{1} (P_{1}) ⟩ J^{μ} (0) N_{2} (P_{2}), X ⟨ N_{2} (P_{2}), X ⟩ J^{ν} (0) N_{1} (P_{1}) .

W^{μν} = i, j \sum H_{ij}^{μν} (Q, q_{⊥}, μ) ϕ_{i / N_{1}} (x, μ) d_{N_{2} / j} (z, μ) + O (Λ_{QCD}^{2} / Q^{2}),

W^{μν} = i, j \sum H_{ij}^{μν} (Q, q_{⊥}, μ) ϕ_{i / N_{1}} (x, μ) d_{N_{2} / j} (z, μ) + O (Λ_{QCD}^{2} / Q^{2}),

ϕ_{q / N_{1}}^{bare} (x)

ϕ_{q / N_{1}}^{bare} (x)

d_{N_{2} / q}^{bare} (z)

W^{μν}

W^{μν}

\times B_{i / N_{1}} (x, b_{⊥}, μ) D_{N_{2} / j} (z, b_{⊥}, μ) S_{ij} (b_{⊥}, μ) + O (q_{T}^{2} / Q^{2})

= i, j \sum H_{ij}^{' μν} (Q, μ) \int d^{2} k_{1 ⊥} d^{2} k_{2 ⊥} d^{2} k_{s ⊥} δ^{(2)} (k_{1 ⊥} + q_{⊥} - k_{2 ⊥} - k_{s ⊥})

\times \tilde{B}_{i / N_{1}} (x, k_{1 ⊥}, μ) \tilde{D}_{N_{2} / j} (z, k_{2 ⊥}, μ) \tilde{S}_{ij} (k_{s ⊥}, μ) + O (q_{T}^{2} / Q^{2}),

B_{q / N_{1}}^{bare} (x, b_{⊥})

B_{q / N_{1}}^{bare} (x, b_{⊥})

\equiv \int \frac{d t}{2 π} e^{- i x t \overset{n}{ˉ} \cdot P_{1}} ⟨ N_{1} (P_{1}) ⟩ \overset{χ}{ˉ}_{n} (t \overset{n}{ˉ} + b_{⊥}) \frac{\neq n ˉ}{2} χ_{n} (0) N_{1} (P_{1})

= \int \frac{d b _{-}}{4 π} e^{- i x b_{-} P_{1 +} /2} ⟨ N_{1} (P_{1}) ⟩ \overset{χ}{ˉ}_{n} (0, b_{-}, b_{⊥}) \frac{\neq n ˉ}{2} χ_{n} (0) N_{1} (P_{1}) .

D_{N_{2} / q}^{bare} (z, b_{⊥})

D_{N_{2} / q}^{bare} (z, b_{⊥})

\equiv X \sum \frac{1}{z} \int \frac{d t}{2 π} e^{i t n \cdot P_{2} / z} Tr ⟨ 0 ⟩ \frac{\neq n}{2} χ_{\overset{n}{ˉ}} (t n + b_{⊥}) N_{2} (P_{2}), X ⟨ N_{2} (P_{2}), X ⟩ \overset{χ}{ˉ}_{\overset{n}{ˉ}} (0) 0

= X \sum \frac{1}{z} \int \frac{d b _{+}}{4 π} e^{i b_{+} P_{2 -} / (2 z)} Tr ⟨ 0 ⟩ \frac{\neq n}{2} χ_{\overset{n}{ˉ}} (b_{+}, 0, b_{⊥}) N_{2} (P_{2}), X ⟨ N_{2} (P_{2}), X ⟩ \overset{χ}{ˉ}_{\overset{n}{ˉ}} (0) 0 .

F_{N_{2} / q}^{bare} (z, b_{⊥} / z) \equiv \int d^{d - 2} P_{2 ⊥} e^{i b_{⊥} \cdot P_{2 ⊥} / z} \tilde{F}_{N_{2} / q}^{bare} (z, P_{2 ⊥}),

F_{N_{2} / q}^{bare} (z, b_{⊥} / z) \equiv \int d^{d - 2} P_{2 ⊥} e^{i b_{⊥} \cdot P_{2 ⊥} / z} \tilde{F}_{N_{2} / q}^{bare} (z, P_{2 ⊥}),

\tilde{F}_{N_{2} / q}^{bare} (z, P_{2 ⊥}) = X \sum \frac{1}{z} \int \frac{d b _{+}}{4 π} \frac{d ^{d - 2} b _{⊥}}{( 2 π ) ^{d - 2}} e^{i b_{+} P_{2 -} / (2 z)} \times Tr ⟨ 0 ⟩ \frac{\neq n}{2} χ_{\overset{n}{ˉ}} (b_{+}, 0, b_{⊥}) N_{2} (P_{2}), X ⟨ N_{2} (P_{2}), X ⟩ \overset{χ}{ˉ}_{\overset{n}{ˉ}} (0) 0 .

\tilde{F}_{N_{2} / q}^{bare} (z, P_{2 ⊥}) = X \sum \frac{1}{z} \int \frac{d b _{+}}{4 π} \frac{d ^{d - 2} b _{⊥}}{( 2 π ) ^{d - 2}} e^{i b_{+} P_{2 -} / (2 z)} \times Tr ⟨ 0 ⟩ \frac{\neq n}{2} χ_{\overset{n}{ˉ}} (b_{+}, 0, b_{⊥}) N_{2} (P_{2}), X ⟨ N_{2} (P_{2}), X ⟩ \overset{χ}{ˉ}_{\overset{n}{ˉ}} (0) 0 .

\tilde{F}_{N_{2} / q}^{bare} (z, P_{2 ⊥}) = \tilde{D}_{N_{2} / q}^{bare} (z, - P_{2 ⊥} / z) .

\tilde{F}_{N_{2} / q}^{bare} (z, P_{2 ⊥}) = \tilde{D}_{N_{2} / q}^{bare} (z, - P_{2 ⊥} / z) .

F_{N_{2} / q}^{bare} (z, b_{⊥} / z) = z^{2 - 2 ϵ} D_{N_{2} / q}^{bare} (z, b_{⊥}) .

F_{N_{2} / q}^{bare} (z, b_{⊥} / z) = z^{2 - 2 ϵ} D_{N_{2} / q}^{bare} (z, b_{⊥}) .

S_{q \overset{q}{ˉ}}^{bare} (b_{⊥}) \equiv \frac{1}{N _{c}} Tr ⟨ 0 ⟩ S_{\overset{n}{ˉ}}^{†} (b_{⊥}) S_{n} (b_{⊥}) S_{n}^{†} (0) S_{\overset{n}{ˉ}} (0) 0,

S_{q \overset{q}{ˉ}}^{bare} (b_{⊥}) \equiv \frac{1}{N _{c}} Tr ⟨ 0 ⟩ S_{\overset{n}{ˉ}}^{†} (b_{⊥}) S_{n} (b_{⊥}) S_{n}^{†} (0) S_{\overset{n}{ˉ}} (0) 0,

S_{n}(x)\equiv\operatorname{\mathcal{P}}\exp\bigg{(}ig_{s}\int_{-\infty}^{0}ds\,n\cdot A_{s}(x+sn)\bigg{)}\,,

S_{n}(x)\equiv\operatorname{\mathcal{P}}\exp\bigg{(}ig_{s}\int_{-\infty}^{0}ds\,n\cdot A_{s}(x+sn)\bigg{)}\,,

collinear: p_{n}

collinear: p_{n}

anti-collinear: p_{\overset{n}{ˉ}}

soft: p_{s}

B_{q / N_{1}} (x, b_{⊥}, μ) D_{N_{2} / q} (z, b_{⊥}, μ) S_{q \overset{q}{ˉ}} (b_{⊥}, μ) = (\frac{b _{T}^{2} Q ^{2}}{b _{0}^{2}})^{- F_{q \overset{q}{ˉ}} (b_{⊥}, μ)} B_{q / N_{1}} (x, b_{⊥}, μ) D_{N_{2} / q} (z, b_{⊥}, μ),

B_{q / N_{1}} (x, b_{⊥}, μ) D_{N_{2} / q} (z, b_{⊥}, μ) S_{q \overset{q}{ˉ}} (b_{⊥}, μ) = (\frac{b _{T}^{2} Q ^{2}}{b _{0}^{2}})^{- F_{q \overset{q}{ˉ}} (b_{⊥}, μ)} B_{q / N_{1}} (x, b_{⊥}, μ) D_{N_{2} / q} (z, b_{⊥}, μ),

b_{T}^{2} = - b_{⊥}^{2}, b_{0} = 2 e^{- γ_{E}} .

b_{T}^{2} = - b_{⊥}^{2}, b_{0} = 2 e^{- γ_{E}} .

\int \frac{d ^{d} k}{( 2 π ) ^{d}} (2 π) δ_{+} (k^{2}) \to \int \frac{d ^{d} k}{( 2 π ) ^{d}} \frac{ν ^{2 η}}{∣2 k _{z} ∣ ^{2 η}} (2 π) δ_{+} (k^{2}),

\int \frac{d ^{d} k}{( 2 π ) ^{d}} (2 π) δ_{+} (k^{2}) \to \int \frac{d ^{d} k}{( 2 π ) ^{d}} \frac{ν ^{2 η}}{∣2 k _{z} ∣ ^{2 η}} (2 π) δ_{+} (k^{2}),

\int \frac{d ^{d} k}{( 2 π ) ^{d}} (2 π) δ_{+} (k^{2}) \to τ \to 0 lim \int \frac{d ^{d} k}{( 2 π ) ^{d}} (2 π) δ_{+} (k^{2}) exp (- b_{0} τ k^{0}) .

\int \frac{d ^{d} k}{( 2 π ) ^{d}} (2 π) δ_{+} (k^{2}) \to τ \to 0 lim \int \frac{d ^{d} k}{( 2 π ) ^{d}} (2 π) δ_{+} (k^{2}) exp (- b_{0} τ k^{0}) .

i = 1 \prod n \int \frac{d ^{d} k _{i}}{( 2 π ) ^{d}} (2 π) δ_{+} (k_{i}^{2}) \to τ \to 0 lim i = 1 \prod n \int \frac{d ^{d} k _{i}}{( 2 π ) ^{d}} (2 π) δ_{+} (k_{i}^{2}) exp (- b_{0} τ i \sum n k_{i}^{0}) .

i = 1 \prod n \int \frac{d ^{d} k _{i}}{( 2 π ) ^{d}} (2 π) δ_{+} (k_{i}^{2}) \to τ \to 0 lim i = 1 \prod n \int \frac{d ^{d} k _{i}}{( 2 π ) ^{d}} (2 π) δ_{+} (k_{i}^{2}) exp (- b_{0} τ i \sum n k_{i}^{0}) .

\mathcal{S}^{\text{bare}}_{q\bar{q}}(b_{\perp},\nu)\equiv\frac{1}{N_{c}}\lim_{\tau\to 0}\operatorname{Tr}\braket{0}{[S^{\dagger}_{\bar{n}}S_{n}](-ib_{0}\tau,-ib_{0}\tau,b_{\perp})\,[S^{\dagger}_{n}S_{\bar{n}}](0)}{0}\Big{|}_{\tau\equiv 1/\nu}\,,

\mathcal{S}^{\text{bare}}_{q\bar{q}}(b_{\perp},\nu)\equiv\frac{1}{N_{c}}\lim_{\tau\to 0}\operatorname{Tr}\braket{0}{[S^{\dagger}_{\bar{n}}S_{n}](-ib_{0}\tau,-ib_{0}\tau,b_{\perp})\,[S^{\dagger}_{n}S_{\bar{n}}](0)}{0}\Big{|}_{\tau\equiv 1/\nu}\,,

\mathcal{B}^{\text{bare}}_{q/N_{1}}(x,b_{\perp},\nu)\equiv\frac{1}{\mathcal{S}_{\text{0b}}}\lim_{\tau\to 0}\int\frac{db_{-}}{4\pi}\,e^{-ixb_{-}P_{1+}/2}\\ \times\braket{N_{1}(P_{1})}{\bar{\chi}_{n}(-ib_{0}\tau,b_{-}-ib_{0}\tau,b_{\perp})\frac{\not{\bar{n}}}{2}\chi_{n}(0)}{N_{1}(P_{1})}\Big{|}_{\tau\equiv 1/\nu}\,,

\mathcal{B}^{\text{bare}}_{q/N_{1}}(x,b_{\perp},\nu)\equiv\frac{1}{\mathcal{S}_{\text{0b}}}\lim_{\tau\to 0}\int\frac{db_{-}}{4\pi}\,e^{-ixb_{-}P_{1+}/2}\\ \times\braket{N_{1}(P_{1})}{\bar{\chi}_{n}(-ib_{0}\tau,b_{-}-ib_{0}\tau,b_{\perp})\frac{\not{\bar{n}}}{2}\chi_{n}(0)}{N_{1}(P_{1})}\Big{|}_{\tau\equiv 1/\nu}\,,

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

aainstitutetext: Zhejiang Institute of Modern Physics, Department of Physics, Zhejiang University, Hangzhou, 310027, Chinabbinstitutetext: School of Physics and State Key Laboratory of Nuclear Physics and Technology, Peking University, Beijing 100871, Chinaccinstitutetext: Center for High Energy Physics, Peking University, Beijing 100871, China

Transverse Parton Distribution and Fragmentation Functions at NNLO: the Quark Case

Ming-Xing Luo

[email protected] b

Xing Wang

[email protected] b

Xiaofeng Xu

[email protected] b,c

Li Lin Yang

[email protected] a

Tong-Zhi Yang

[email protected] a

and Hua Xing Zhu

[email protected]

Abstract

We revisit the calculation of perturbative quark transverse momentum dependent parton distribution functions and fragmentation functions using the exponential regulator for rapidity divergences. We show that the exponential regulator provides a consistent framework for the calculation of various ingredients in transverse momentum dependent factorization. Compared to existing regulators in the literature, the exponential regulator has a couple of advantages which we explain in detail. As a result, the calculation is greatly simplified and we are able to obtain the next-to-next-to-leading order results up to $\mathcal{O}(\epsilon^{2})$ in dimensional regularization. These terms are necessary for a higher order calculation which is made possible with the simplification brought by the new regulator. As a by-product, we have obtained the two-loop quark jet function for the Energy-Energy Correlator in the back-to-back limit, which is the last missing ingredient for its N3LL resummation.

Keywords:

SCET, beam function, QCD corrections

1 Introduction

Parton distribution functions (PDFs) and fragmentation functions (FFs) describe the partonic contents of hadrons. They are of fundamental importance in quantum chromodynamics (QCD) Collins:2011zzd ; Angeles-Martinez:2015sea . They enter factorization formulas for scattering processes involving hadrons, and are essential for comparing theoretical predictions for the cross sections against experimental measurements. In most cases, the observables are only sensitive to the longitudinal momenta of the partons, and the transverse momenta can be integrated over, leading to the so-called “collinear” PDFs and FFs. However, in certain regions of phase space, the transverse momenta of the partons become relevant, and one needs the transverse momentum dependent PDFs (TMDPDFs) and FFs (TMDFFs) in the corresponding factorization formulas. This is case for the small transverse momentum ( $Q_{T}$ ) region in the Drell-Yan process Dokshitzer:1978yd ; Parisi:1979se ; Collins:1984kg ; Arnold:1990yk ; Ladinsky:1993zn ; Bozzi:2010xn ; Becher:2011xn ; Bizon:2019zgf ; Bertone:2019nxa , and also for similar regions in, e.g., semi-inclusive deep-inelastic scattering (SIDIS) Ji:2004wu ; Ji:2004xq ; Su:2014wpa ; Kang:2015msa ; Liu:2018trl , electron-positron annihilation to hadrons and jets Collins:1981uk ; Collins:1981va ; Neill:2016vbi ; Gutierrez-Reyes:2018qez ; Gutierrez-Reyes:2019vbx ; Gutierrez-Reyes:2019msa , Higgs boson production Berger:2002ut ; Bozzi:2005wk ; Gao:2005iu ; Echevarria:2015uaa ; Neill:2015roa ; Bizon:2017rah ; Chen:2018pzu ; Bizon:2018foh , top quark pair production Zhu:2012ts ; Li:2013mia ; Catani:2014qha ; Catani:2018mei , as well as Energy-Energy Correlator (EEC) in the back-to-back limit at both lepton and hadron colliders Moult:2018jzp ; Gao:2019ojf . To improve the theoretical predictions for these observables, it is desirable to have precise knowledges about these basic objects.

TMDPDFs and TMDFFs can be defined as hadronic matrix elements of bilinear quark or gluon field operators with a measured transverse momentum $\vec{k}_{\perp}$ (or a transverse separation $\vec{b}_{\perp}$ in position space). If the transverse momentum $\vec{k}_{\perp}\sim\Lambda_{\text{QCD}}$ , the TMDPDFs and TMDFFs are essentially non-perturbative, and can only be extracted from experimental data or calculated using lattice methods. On the other hand, if $\vec{k}_{\perp}\gg\Lambda_{\text{QCD}}$ , the TMDPDFs and TMDFFs can be related to the collinear PDFs and FFs via perturbatively calculable matching coefficients. These coefficients are known at the next-to-next-to-leading order (NNLO) for the TMDPDFs Catani:2011kr ; Gehrmann:2012ze ; Gehrmann:2014yya ; Echevarria:2016scs and TMDFFs Echevarria:2016scs . They have played an important role in a number of cutting-edge calculations, including precision predictions for the Drell-Yan process and Higgs boson production at small transverse momentum Chen:2018pzu ; Bizon:2018foh ; Bizon:2019zgf , and NNLO calculations for top quark pair production using the $Q_{T}$ subtraction method Catani:2019iny ; Catani:2019hip .

In this work, we revisit the calculation of the matching coefficients for TMDPDFs and TMDFFs at NNLO. We consider the quark TMDPDFs and TMDFFs in this paper, while the gluon case is left to a forthcoming article. There are several new elements in our calculation compared to those in the literature:

•

We employ the exponential regulator for rapidity divergences Li:2016axz . Rapidity divergences or “collinear anomalies” appear in the calculation of individual TMD functions in factorization formulas, which are cancelled in physical observables. These divergences are not regularized by dimensional regularization, and additional regulators need to be introduced Collins:1984kg ; Ji:2004wu ; Collins:2011zzd ; Becher:2010tm ; Becher:2011dz ; Chiu:2012ir ; Chiu:2009yx ; Echevarria:2015byo ; Ebert:2018gsn . The exponential regulator has been shown to be particularly suitable in the calculation of TMD soft functions, as demonstrated in the recent next-to-next-to-next-to-leading order (N3LO) calculation Li:2016ctv . We show in this work that the exponential regulator can also be used to calculate TMDPDFs and TMDFFs, which are more complicated objects than TMD soft functions. Our results show that the exponential regulator is a consistent rapidity regulator in both the soft and collinear sectors.

•

We develop systematic calculation method based on modern techniques for loop integrals, such as integration-by-parts (IBP) identities Chetyrkin:1981qh ; Laporta:2001dd and differential equations Bern:1993kr ; Gehrmann:1999as ; Henn:2013pwa . Our method paves the way to calculate TMDPDFs and TMDFFs at N3LO.

•

We obtain the bare NNLO TMDPDFs and TMDFFs up to $\mathcal{O}(\epsilon^{2})$ , where $\epsilon$ is the dimensional regulator. They directly contribute to TMDPDFs and TMDFFs at N3LO upon renormalization.

•

Our results for TMDPDFs agree with previous calculations Gehrmann:2012ze ; Gehrmann:2014yya ; Echevarria:2016scs , but we find a small discrepancy for the TMDFFs compared to those presented in Ref. Echevarria:2016scs . We have performed several consistency checks on our results to make sure that they are correct.

•

As a by-product, we obtain the NNLO quark jet function relevant for the resummation of EEC in the back-to-back limit. This is the last missing ingredient for this resummation at the next-to-next-to-next-to-leading logarithmic (N3LL) accuracy.

This paper is organized as follows. In Section 2 we introduce the definitions of quark TMDPDFs and TMDFFs in the context of the SIDIS process, and discuss the exponential regulator for rapidity divergences. In Section 3 and 4 we perform the calculation of the quark TMDPDFs and TMDFFs at NNLO using the exponential regulator. In Section 4 we also use the results for TMDFFs to compute the two-loop jet function for EEC in the back-to-back limit. This is by itself a new result of our paper, and also serves as a cross-check of our results. We conclude in Section 5.

2 Transverse momentum dependent factorization

2.1 Kinematics and factorization

In this section, we briefly review the formalism of transverse momentum dependent factorization and introduce the definitions of TMDPDFs and TMDFFs. For our purpose, it is easiest to consider (unpolarized) SIDIS which involves hadrons in both the initial state and the final state. In SIDIS, a hadron $N_{1}$ with momentum $P_{1}^{\mu}$ is probed by a virtual photon $\gamma^{*}$ with momentum $q^{\mu}$ and produces a jet containing a specific hadron $N_{2}$ with momentum $P_{2}^{\mu}$ . We define the kinematic invariants

[TABLE]

We introduce two light-like 4-vectors $n$ and $\bar{n}$ satisfying $n^{2}=\bar{n}^{2}=0$ and $n\cdot\bar{n}=2$ , such that we can decompose any 4-vector $k^{\mu}$ as

[TABLE]

When quoting the components of a 4-vector, we use $k=(k_{+},k_{-},k_{\perp})$ . The scalar product of two 4-vectors is given by

[TABLE]

In the hadron frame and ignoring the hadron masses, we have

[TABLE]

where

[TABLE]

and we define $q_{\perp}^{2}\equiv-q_{T}^{2}$ .

The hadronic tensor is defined as

[TABLE]

In the region $q_{T}\sim Q\gg\Lambda_{\text{QCD}}$ , the hadronic tensor can be factorized into products of hard kernels with collinear PDFs and FFs:

[TABLE]

where we have suppressed the dependence of the hard kernel on other kinematic variables. In the language of soft-collinear effective theory (SCET) Bauer:2000ew ; Bauer:2000yr ; Bauer:2001yt ; Bauer:2002nz ; Beneke:2002ph , the collinear PDFs and FFs can be defined as matrix elements of gauge-invariant collinear fields. For example, the bare quark collinear PDF and FF are defined by Collins:1981uw ; Bauer:2002nz ; Collins:2011zzd

[TABLE]

where $\chi_{n}$ and $\chi_{\bar{n}}$ are the gauge-invariant collinear quark fields along the $n$ and $\bar{n}$ directions, respectively. We have assumed dimensional regularization with $d=4-2\epsilon$ . The collinear PDF $\phi_{q/N_{1}}(x,\mu)$ describes (in a sense) the probability distribution of finding the quark $q$ with momentum fraction $x$ inside the fast-moving hadron $N_{1}$ . The collinear FF $d_{N_{2}/q}(z,\mu)$ , on the other hand, describes the probability distribution of finding the hadron $N_{2}$ with momentum fraction $z$ inside the jet initiated by the quark $q$ .

If $q_{T}\ll Q$ , however, the above picture of collinear factorization breaks down due to the appearance of large logarithms of $q_{T}/Q$ in the hard kernel $H_{ij}$ . One should instead rely on TMD factorization of the form

[TABLE]

where $\mathcal{B}_{i/N_{1}}$ , $\mathcal{D}_{N_{2}/j}$ and $\mathcal{S}_{ij}$ are TMDPDFs, TMDFFs and TMD soft functions in the impact parameter space, with $b_{\perp}$ the impact parameter; while $\tilde{\mathcal{B}}_{i/N_{1}}$ , $\tilde{\mathcal{D}}_{N_{2}/j}$ and $\tilde{\mathcal{S}}_{ij}$ are their counterparts in the transverse momentum space. For our purpose, we only consider $i,j$ being quarks and anti-quarks.

The quark TMDPDF $\tilde{B}_{q/N_{1}}(x,k_{1\perp},\mu)$ describes the probability distribution of find a quark with momentum fraction $x$ and transverse momentum $k_{1\perp}$ inside the hadron $N_{1}$ , as depicted in the left plot of Figure 1. Naively, the bare quark TMDPDF can be defined by

[TABLE]

Similarly, the bare quark TMDFF may be defined as

[TABLE]

Note that in the above definition, $k_{2\perp}$ represents the transverse momentum of the quark in the hadron frame (where $N_{2}$ has zero transverse momentum). In practice, it is also useful to define the TMDFFs in the parton frame where the quark has zero transverse momentum. In the parton frame, $N_{2}$ now has a non-zero transverse momentum $P_{2\perp}$ which is related to $k_{2\perp}$ by $P_{2\perp}=-zk_{2\perp}$ . The parton frame quark TMDFF is then

[TABLE]

with

[TABLE]

Here, $P_{2\perp}$ is defined with respect to the axis chosen such that the total transverse momentum of $N_{2}$ and $X$ is zero. It is easy to show that

[TABLE]

The function $\tilde{\mathcal{F}}^{\text{bare}}_{N_{2}/q}(z,P_{2\perp})$ represents the probability distribution of finding a hadron $N_{2}$ with momentum fraction $z$ and transverse momentum $P_{2\perp}$ inside the jet initiated by the quark $q$ , as depicted in the right plot of Figure 1. From the above definitions, it is easy to see that

[TABLE]

Finally, the quark TMD soft function is given by the vacuum expectation value of a soft Wilson loop

[TABLE]

where the soft Wilson line is defined by

[TABLE]

with $A_{s}$ the soft gluon field in SCET.

2.2 Rapidity divergences and the exponential regulator

While the TMD factorization formula (9) makes some sense, the TMDPDF (10), TMDFF (11) and TMD soft function (16) are actually ill-defined due to the appearance of rapidity divergences which are not regularized in dimensional regularization. These divergences cancel when one combines the 3 functions in the factorization formula (9) to calculate physical observables. However, they also carry important information, just like the relationship between ultraviolet (UV) divergences and the renormalization group.

The rapidity divergences arise due to the fact that the collinear modes and soft modes have the same typical off-shellness around $q_{T}^{2}$ . More precisely, in the $q_{T}\ll Q$ limit we have the relevant momentum regions

[TABLE]

where $\lambda=q_{T}/Q\ll 1$ . The effective field theory describing these modes are sometimes called $\text{SCET}_{\text{II}}$ . The collinear modes and the soft mode are related by a boost in the $n$ or $\bar{n}$ direction. As a result, they cannot really be separated by a boost invariant regulator such as dimensional regularization. A brute-force separation as done in Eq. (9) then leads to inconsistencies manifesting themselves as rapidity divergences.

To deal with the rapidity divergences, one needs to introduce a regulator in addition to dimensional regularization. This however leads to another subtle issue. Any such regulator necessarily reintroduces a logarithmic dependence on the hard scale $Q$ into integrals in the collinear and anti-collinear regions through $\bar{n}\cdot P_{1}$ and $n\cdot P_{2}$ , which was supposed to be factorized out into the hard function $H^{\prime}$ in Eq. (9). This fact is sometimes called “collinear anomaly” or “factorization anomaly” in the literature MartinBeneke ; Becher:2010tm . Nevertheless, using the structure of the rapidity divergences, it can be shown that these $Q$ -dependence can be extracted and exponentiated to all orders. After such a “re-factorization”

[TABLE]

where

[TABLE]

The functions $B_{q/N_{1}}$ and $D_{N_{2}/q}$ can be regarded as the “genuine” quark TMDPDF and TMDFF which are free from rapidity divergences and are also independent of $Q$ . The exponent function $F_{q\bar{q}}$ is closely related to the so-called Collins-Soper kernel Collins:2011zzd . It has been known perturbatively to three loops Li:2016ctv ; Vladimirov:2016dll . Very recently, there are proposals to compute it non-perturbatively on the lattice Ebert:2018gzl ; Ebert:2019okf .

In the literature, there are a variety of ways to regularize the rapidity divergences Ji:2004wu ; Chiu:2009yx ; Collins:2011zzd ; Chiu:2012ir ; Becher:2011dz ; Echevarria:2015byo ; Ebert:2018gsn . In this paper, we consider the so-called exponential regulator Li:2016axz which was used to calculate the TMD soft function to the N3LO. We will show that it is a consistent regularization scheme also for the TMDPDFs and TMDFFs. Before discussing the exponential regulator, we briefly review the $\eta$ -regulator of Ref. Chiu:2012ir which shares many similarities. At the next-to-leading order (NLO), the $\eta$ -regulator amounts to the subsitution

[TABLE]

for the phase-space integrals over the real gluon momentum $k^{\mu}$ , where $\delta_{+}(k^{2})=\theta(k^{0})\delta(k^{2})$ . The rapidity divergences appear as $1/\eta$ poles which can be subtracted in the same way as renormalizing the UV divergences. After the subtraction, the TMDPDFs, TMDFFs and TMD soft functions still depend on the “rapidity scale” $\nu$ . For the TMDPDFs and TMDFFs, the natural rapidity scale is $\nu\sim Q$ , while for the TMD soft functions $\nu\sim q_{T}$ . The evolution equations of these functions with respect to $\nu$ can be used to exponentiate the rapidity logarithms $\ln(b_{T}^{2}Q^{2})$ leading to the refactorization in Eq. (19).

While the above $\eta$ -regulator is conceptually simple, it is not easy to implement in higher order calculations beyond NLO. For example, the regulator has to be carefully applied to maintain non-Abelian exponentiation in the soft sector Chiu:2012ir ; Luebbert:2016itl . In particular, when there are two real gluon emissions with momenta $k_{1}$ and $k_{2}$ , it is different to apply the regulator on $k_{1z}+k_{2z}$ as a whole, or on $k_{1z}$ and $k_{2z}$ separately. Recently, a new regulator for rapidity divergences called “exponential regulator” has been proposed in Ref. Li:2016axz , which leads to the same rapidity evolution equations as the $\eta$ -regulator, and is easier for higher order calculations. In momentum space, the new rapidity regulator is simply multiplying each soft/collinear phase space measure by an exponential factor

[TABLE]

Note that the $\tau\to 0$ limit has to be taken after integration. Beyond NLO, when there are multiple soft/collinear partons, the regularization simple becomes

[TABLE]

Due to the exponential form, the multiple emission case naturally factorizes into products of single emissions. Therefore, non-Abelian exponentiation is manifestly preserved by this regulator. An important feature of the exponential regulator is that it leads to enormous simplification in perturbative calculations, as demonstrated by the calculation of TMD soft functions at N3LO in Ref. Li:2016ctv .

The exponential regulator also admits simple operator definitions for the TMD functions. For example, the quark TMD soft function is defined as

[TABLE]

where the rapidity regularization procedure is understood as keeping non-vanishing terms in the limit of $\tau\to 0$ (including the $\log\tau$ terms which are the manifestation of rapidity divergences), and then identify the rapidity scale as $\nu=1/\tau$ . No subtraction is needed and the rapidity divergence are now renormalized. The remaining results depend on logarithms of the rapidity scale $\nu$ .

Similarly, the exponentially regularized quark TMDPDF and TMDFF are defined as

[TABLE]

and

[TABLE]

Note that for both the TMDPDF and TMDFF, we need to perform a zero-bin subtraction to avoid double-counting between the collinear sectors and the soft sector. The zero-bin soft function is the same as the TMD soft function

[TABLE]

Having operator definitions Eqs. (24), (25) and (26) for the TMD functions could be advantageous for studying non-perturbative aspects of TMD physics. In this work we focus on the perturbative part of the TMDPDF and TMDFF.

2.3 Renormalization and perturbative matching

For large impact parameter $b_{T}\sim 1/\Lambda_{\text{QCD}}$ , the TMDPDFs and TMDFFs are dominated by long distance contributions and are genuine non-perturbative objects. In this work, we are interested in the semi-perturbative region $b_{T}\ll 1/\Lambda_{\text{QCD}}$ . In this region the TMDPDF admits an operator product expansion

[TABLE]

where $\phi_{i/N}$ is the (renormalized) collinear PDF of parton $i$ , and $\mathcal{I}_{qi}$ is a perturbatively calculable matching coefficient function describing the splitting of the parton $i$ into the quark $q$ . Similarly, the TMDFF can also be factorized as

[TABLE]

with perturbatively calculable coefficient functions $\mathcal{C}_{iq}$ describing the fragmentation of the quark $q$ into the parton $i$ .

The functions $\mathcal{I}_{qi}$ and $\mathcal{C}_{iq}$ will be the main objects we are going to study in this work. As indicated by the superscript “bare” in Eqs. (28) and (29), there are UV divergences which require renormalization. For the TMDPDF, we define the renormalization factor according to

[TABLE]

where we have used $\otimes$ to denote the convolution in Eq. (28). The matching coefficients $\mathcal{I}_{qi}$ do not depend on the external state $N$ , and can therefore be calculated with $N$ replaced by a partonic state $j=q$ or $g$ . We can then extract $\mathcal{I}_{qi}$ by calculating $\mathcal{B}^{\text{bare}}_{q/j}$ , performing the renormalization and subtracting the partonic collinear PDFs $\phi_{i/j}$ . Up to the NNLO, the partonic collinear PDFs are given by

[TABLE]

where $P_{ij}^{(0)}$ is the LO splitting kernel and $\beta_{0}$ is the LO beta function.

After renormalization, the TMDPDF obeys a renormalization group equation (RGE)

[TABLE]

where $\Gamma^{\text{cusp}}$ is the usual cusp anomalous dimension and $\gamma^{B}$ is the non-cusp anomalous dimension for the TMDPDF, whose perturbative expansions are collected in the Appendix. From the above equation and the famous DGLAP equation

[TABLE]

one can deduce the RGEs for the coefficient functions as

[TABLE]

Besides the normal RGE, the TMDPDF and the coefficient functions also satisfy the rapidity evolution equation Chiu:2012ir

[TABLE]

The rapidity anomalous dimension $\gamma^{R}$ is known to three loops in QCD Li:2016ctv ; Vladimirov:2016dll . For our purpose, we need the first two orders which are given by111Note that the convention here differ by a factor of 2 from Ref. Li:2016ctv .

[TABLE]

The renormalization equations (34) and (35) can be used to determine all the renormalization and rapidity scale dependent terms for the coefficient functions in perturbation theory. Throughout this paper, we organize perturbative expansions of various functions in powers of $\alpha_{s}/(4\pi)$ . For example

[TABLE]

Here and below we introduce two logarithms

[TABLE]

Up to $\mathcal{O}(\alpha_{s}^{2})$ , we then have

[TABLE]

Similarly for the TMDFF, the UV renormalization is given by

[TABLE]

Note that the procedure of renormalization and matching is easier to be done with the parton frame $\mathcal{F}_{N/q}$ instead of the hadron frame $\mathcal{D}_{N/q}$ (used in Echevarria:2015usa ). To extract the coefficient functions $\mathcal{C}_{iq}$ , we calculate the bare TMDFFs with external parton states, and the partonic collinear FFs up to NNLO are given by

[TABLE]

where $P^{T}_{ij}(z)$ are the time-like splitting kernels which will be presented in the Appendix. The renormalized $\mathcal{C}_{iq}$ functions satisfy the evolution equations

[TABLE]

and

[TABLE]

At this point, it is worth noting that the product of the TMDPDF, TMDFF and the TMD soft function is independent on the rapidity scale $\nu$ as expected, namely

[TABLE]

where we have used

[TABLE]

It can also be shown that the $\mu$ -dependence of this product is cancelled by that of the hard function (which does not know about the rapidity divergences), such that the physical observables are independent of the renormalization scale. To see that we recall the RGEs of the hard and soft functions

[TABLE]

We note that the cancellation happens since $\gamma^{B}+\gamma^{H}-\gamma^{S}=0$ and

[TABLE]

By using the evolution equations, we can derive the scale-dependent part of $\mathcal{C}_{iq}$ . Up to the NNLO we have

[TABLE]

Note that here we have used the same symbol $L_{Q}$ as in Eq. (38) to denote a different meaning:

[TABLE]

which can be regarded as the crossing $P_{1+}\to P_{2-}$ and $x\to 1/z$ .

2.4 Rapidity renormalization group and re-factorization

We now use the rapidity evolution equations of the TMDPDF, TMDFF and TMD soft function to derive the re-factorization formula (19). From the perturbative matching coefficients, it is evident that for the TMDPDF and TMDFF, the natural rapidity scale is $\nu\sim xP_{1+}\sim P_{2-}/z\sim Q$ , while for the TMD soft function the natural choice is $\nu\sim b_{0}/b_{T}\sim q_{T}$ . In order to reconcile these different choices, we may use the rapidity RGE (45) for the TMD soft function to evolve it from $\nu=b_{0}/b_{T}$ to $\nu=Q$ . The result is

[TABLE]

where

[TABLE]

where we have used $\gamma_{0}^{R}=0$ . The “genuine” quark TMDPDF and TMDFF which are free from rapidity divergences and are independent of the hard scale can then be defined as

[TABLE]

These are essentially the functions appearing in the re-factorization formula (19).

3 Quark TMDPDF with the exponential regulator

In this section, we calculate the perturbative matching coefficients of the quark TMDPDF at NLO and NNLO using the exponential regulator. While these results are known to order $\epsilon^{0}$ in the literature Gehrmann:2012ze ; Gehrmann:2014yya ; Echevarria:2016scs , we are able to obtain higher order terms in $\epsilon$ . The calculation with the exponential regulator is also much simpler and more systematic, which makes it possible to be extended to N3LO.

3.1 Quark TMDPDF at NLO

In this subsection, we briefly discuss the NLO results with the exponential regulator. While the calculation is straightforward, it illustrates the basic procedure and some interesting features of the regularization scheme.

We begin with the bare TMDPDFs before zero-bin subtraction. According to the definition in Eqs. (25), the TMDPDFs at NLO are given by the cut diagrams in Fig. 2. The result can be written as

[TABLE]

where we have changed to the notation that $\vec{k}_{T}$ denotes the transverse components of $k_{\perp}^{\mu}$ , but with Euclidean signature such that

[TABLE]

The $d$ -dimensional splitting amplitudes are given by

[TABLE]

Using the delta function for $k_{+}$ and the on-shell condition, we can write the exponential regulator as

[TABLE]

At this stage we can already drop the second term proportional to $\tau(1-x)$ in the exponent, as it gives no contribution in the limit $\tau\to 0$ . They might be relevant for subleading power corrections Ebert:2018gsn . The first term involving $\tau/(1-x)$ in the exponent provides the main service of regularizing the rapidity divergences. To see how that happens, we note that the rapidity divergence appears here as a singularity as $x\to 1$ in the $q\to q$ splitting amplitude in Eq. (55). The exponential regulator provides a suppression in the $x\to 1$ limit, and turns this singularity into a regularized distribution according to

[TABLE]

Applying the above equation to $\mathcal{B}_{q/q}^{{(1)},\text{bare},\text{unsub}}$ , we find

[TABLE]

The $\vec{k}_{T}$ integral can be easily performed with the help of the generating integral

[TABLE]

In particular, we have

[TABLE]

where $\psi(x)=\Gamma^{\prime}(x)/\Gamma(x)$ . Applying the above formulas to Eq. (58) then gives $\mathcal{B}_{q/q}^{{(1)},\text{bare},\text{unsub}}$ exact in $\epsilon$ , and similarly for $\mathcal{B}_{q/g}^{{(1)},\text{bare},\text{unsub}}$ . It is easy to expand the results to any order in $\epsilon$ . In particular, the $\mathcal{O}(\epsilon^{2})$ terms will be used for the NNLO calculations later, while the $\mathcal{O}(\epsilon^{4})$ terms are relevant for the calculations at N3LO.

We now need to subtract the collinear PDFs in Eq. (2.3) to obtain the coefficient functions $\mathcal{I}^{{(1)},\text{bare},\text{unsub}}_{qi}$ , and then perform the zero-bin subtraction, where the NLO zero-bin contribution is given by

[TABLE]

where $L_{\nu}=\ln(\nu^{2}/\mu^{2})$ . After the subtraction, we find up to $\mathcal{O}(\epsilon^{0})$

[TABLE]

where $L_{Q}=2\ln(xp_{+}/\nu)$ , and

[TABLE]

We now renormalize the UV divergences in the $\overline{\text{MS}}$ scheme with the renormalization factor

[TABLE]

and find

[TABLE]

Comparing the above form with Eq. (39), we can extract the (renormalization and rapidity) scale-independent part of the NLO coefficients

[TABLE]

Remarkably, in the exponential regularization scheme, the scale independent coefficients are regular in the soft limit $x\to 1$ . As will be explicitly shown below, at NNLO there are $1/(1-x)_{+}$ distributions in the $\mu$ -independent part, but these terms are governed by the rapidity anomalous dimension and depend on the rapidity scale $\nu$ . In general this is true even at high orders in perturbation theory Echevarria:2016scs ; Lustermans:2016nvk .

3.2 Quark TMDPDF at NNLO

We now turn to the NNLO calculations. As before, we begin with the bare TMD coefficient functions before zero-bin subtraction. At NNLO, diagrammatically there are two kinds of contributions. One is the interference of the LO amplitude with the diagrams containing one loop and one real emission, i.e., the so-called real-virtual (RV) contribution. The other is the square of the diagrams with two real emissions, i.e., the so-called double real contribution (RR). We will discuss these two contributions one-by-one in the following.

3.2.1 The real-virtual contribution

We adopt the light-cone gauge $n\cdot A=0$ where the relevant cut diagrams for the real-virtual contribution are depicted in Figure 3. Note that with both the exponential regulator and the analytic regulator used in Gehrmann:2012ze ; Gehrmann:2014yya , the loop integral does not need to be regularized. Therefore the treatments of the loop amplitude are rather similar. After performing the Dirac algebras and partial fractioning, there remain two classes of scalar integrals as shown in Fig. 4, which are given by

[TABLE]

where $q=p-k$ and we make the $+i\epsilon$ prescription for all propagators implicit. The results of these integrals have already been given in Gehrmann:2014yya and we do not repeat them here. After the loop integration, the results are functions of $x=\bar{n}\cdot q/\bar{n}\cdot p$ and $k_{T}^{2}$ . The remaining integral over $k_{T}^{2}$ can be carried out in the same way as the NLO calculation using Eq. (59).

3.2.2 The double real contribution

At NNLO, the double real contribution is the most troublesome one to calculate. We will show that with the exponential regulator, we can apply many modern techniques for loop integrals. It is therefore possible to extend the calculation method to higher orders.

We use QGRAF Nogueira:1991ex to generate the relevant Feynman diagrams in the light-cone gauge, which are shown in Fig. 5. We then use FORM Vermaseren:2000nd to manipulate the squared amplitudes, and write them as integrals over the two cut momenta which we denote as $k_{1}$ and $k_{2}$ . We now need to apply the exponential regulator, and the integral measure then becomes

[TABLE]

It is useful to introduce an identity Zhu:2014fma

[TABLE]

and rewrite the integral measure as

[TABLE]

Now the integration over $k_{1}$ and $k_{2}$ does not produce rapidity divergences and can be performed with usual techniques. Note that this fact holds also beyond NNLO where more than two cut momenta are present, due to the exponential form of the regulator.

We can now use the delta function to integrate over $k_{2}$ , and rename $k_{1}$ as $l$ . The double real contribution can then be written in the form

[TABLE]

where ${\cal M}_{qi}$ is the squared amplitude. We will first integrate over $l$ using the methods of reverse unitarity Anastasiou:2002yz , integration-by-parts (IBP) Chetyrkin:1981qh and differential equations Bern:1993kr ; Gehrmann:1999as ; Henn:2013pwa . The relevant topologies are given by (the square of) the diagrams shown in Figure 6. There are 4 topologies for the $l$ -integrals, which are defined by

[TABLE]

where we use the subscript “cut” to label the cut propagators Anastasiou:2002yz . Integrals in each topology are further reduced to a set of Master Integrals (MIs) by IBP identities Chetyrkin:1981qh . In this work, we use FIRE5 Smirnov:2014hma and LiteRed Lee:2012cn to perform the reduction. In total we have 6 MIs which can be chosen as

[TABLE]

where the normalization factor

[TABLE]

We will use the method of differential equations to evaluate these MIs. For that purpose we have introduced the rescale factors $w_{i}$ to convert the MIs into a canonical basis Henn:2013pwa . They are given by

[TABLE]

where the dimensionless variables $x$ and $y$ are defined as

[TABLE]

The factors $\omega_{i}$ can be easily obtained using an in-house code or the program package CANONICA Meyer:2017joq which implements the algorithm of Meyer:2016slj .

Among all the MIs, $F_{1}$ , $F_{2}$ , $F_{4}$ and $F_{6}$ are easy to be evaluated in closed form

[TABLE]

where ${}_{2}F_{1}$ is the hypergeometric function. $F_{3}$ and $F_{5}$ depend on both $x$ and $y$ and are more difficult to calculate. We can construct the differential equations of them with respect to $y$

[TABLE]

which are in the so-called canonical form Henn:2013pwa . Given their boundary conditions at $y=0$

[TABLE]

it is easy to solve the differential equations order-by-order in $\epsilon$ in terms of Goncharov multiple polylogarithms (GPLs). We have obtained the solutions up to weight 6, which will be sufficient for a future N3LO calculation.

The next step is then to perform the remaining integration over $k$ in Eq. (71). The $k_{+}$ integral can be done using the delta function. And the $k_{-}$ integral can be changed to use the $y$ variable through

[TABLE]

Note that we now have singularities at $y\to 0$ or $x\to 1$ , which are both manifestations of rapidity divergences. These overlapping singularities often make high order perturbative calculations difficult due to the fact that the regularized integrand is often a complicated function of $x$ and $y$ . In our scheme, the regularization is provided by

[TABLE]

where the $y\to 0$ and $x\to 1$ limits are both exponentially suppressed. To perform the integration over $y$ , we expand the above exponential regulator in terms of delta functions and plus-distributions according to Eq. (57) and

[TABLE]

The $y$ integration can now be done using the package HyperInt Panzer:2014caa and the $\vec{k}_{T}$ integration can again be evaluated with the help of Eq. (59). After the integration, the results can be expressed in terms of Harmonic PolyLogarithms (HPLs) Remiddi:1999ew of the variable $x$ . We use the program package HPL Maitre:2005uu to deal with these functions.

3.2.3 Final results at NNLO

Combining the real-virtual and double-real contributions, we obtain the bare un-subtracted NNLO TMDPDF. We then perform the zero-bin subtraction to remove double-counting between the collinear and soft sectors, and apply the usual $\alpha_{s}$ renormalization and operator renormalization $Z^{B}_{q}$ to remove the UV divergences. We have reproduced all the renormalization and rapidity scale dependent parts in Eq. (39), and the scale independent NNLO coefficients $I^{(2)}_{qi}(x)$ are given by

[TABLE]

where $q^{\prime}$ is a light quark flavor different from $q$ , and we have used the shorthand notation

[TABLE]

with $H$ being HPLs. The $p_{qi}(x)$ functions are related to the DGLAP splitting kernels and are collected in the Appendix. We note that with the exponential regulator, the scale-independent part of the TMDPDF does not involve $\delta(1-x)$ terms, and the coefficients of $1/(1-x)_{+}$ is determined by the rapidity anomalous dimension $2\gamma_{1}^{R}$ given in Eq. (36). We have compared our results to those in the literature Gehrmann:2012ze ; Gehrmann:2014yya ; Echevarria:2016scs and found full agreement. We have also obtained the bare NNLO TMDPDFs through to $\mathcal{O}(\epsilon^{2})$ , which is required for a future N3LO calculation. Their expressions are quite lengthy and we choose to put them in an electronic file attached with the arXiv submission of this paper.

4 Quark TMDFF with the exponential regulator

We now turn to the quark TMDFF. Technically, it is very similar to the TMDPDF. The squared amplitudes are related via a crossing symmetry. The only subtlety is that one may perform the calculations in the hadron frame or in the parton frame. The two results should be related according to Eqs. (14) and (15). We have explicitly performed the two calculations and confirmed those relations.

In the hadron frame, one has, at a given order in perturbation theory

[TABLE]

where $\mathcal{M}_{iq}$ is the squared amplitude for the $q\to i$ splitting, $p$ is the momentum of the observed hadron in the $\bar{n}$ direction, $l_{i}$ are loop momenta and $k_{i}$ are momenta of real emissions, $\vec{k}_{T,\text{hf}}$ denotes the total transverse momentum of real emissions in the hadron frame. In the parton frame, one has instead

[TABLE]

where $\vec{k}_{T,\text{pf}}$ denotes the total transverse momentum of real emissions in the parton frame. The only differences with respect to the hadron frame formula are the additional factor of $1/z$ in the Fourier transform, and the different definition of $\vec{k}_{T}$ . They are related by $\vec{k}_{T,\text{pf}}=z\vec{k}_{T,\text{hf}}$ .

4.1 Quark TMDFF at NLO and NNLO

We now present some details about the calculation of the $q\to g$ fragmentation function at NLO in the hadron frame. In the light-cone gauge, there is only one cut diagram contributing, as shown in Fig. 7. The squared amplitude can be straightforwardly obtained from the diagram, or can be related to the $g\to q$ splitting amplitude in Eq. (55) via a crossing symmetry. We have

[TABLE]

with

[TABLE]

The above equation is the manifestation of the Gribov-Lipatov relation.

For the $q\to g$ fragmentation at NLO, there are no rapidity divergences, and therefore we do not need to introduce the exponential regulator. The bare TMDFF in the hadron frame then reads

[TABLE]

The above integral is similar to the ones appearing in the calculation of TMDPDFs. The result is

[TABLE]

We now need to proceed with the matching procedure (29), where one should pay attention to the prefactor $z^{2-2\epsilon}$ , which will produce logarithms of $z$ when expanding in $\epsilon$ . We have

[TABLE]

Performing matching and renormalization as in Eq. (29), we then obtain

[TABLE]

Note that the scale-dependent part agrees with Eq. (48).

We perform the calculation for the $q\to q$ fragmentation in a similar manner, where we need to use the exponential regulator for the rapidity divergences. The scale-independent coefficients at NLO are then given by

[TABLE]

where we use the shorthand notation

[TABLE]

The NNLO calculations proceed in an analogous way, and we do not repeat the details here. The results are

[TABLE]

Again, we find that the scale independent parts do not contain $\delta(1-z)$ terms, and the $1/(1-z)_{+}$ terms are determined by the rapidity anomalous dimension. We can convert our results to the convention of Ref. Echevarria:2015usa and compare with the results in that work. We find that the results agree for the splitting processes $q\to q^{\prime}$ , $q\to\bar{q}$ and $q\to g$ . However, for $q\to q$ , there is a small difference concerning a term $C_{A}C_{F}\,\pi^{4}\,\delta(1-z)$ . In our framework, this term comes from the TMD soft function, which is universal for the TMDPDF and TMDFF. To address this discrepancy, we have performed several independent checks. The strongest check of our calculation comes from the calculation of the two-loop jet function of EEC in the back-to-back limit, which we shall explain in the next subsection.

4.2 Jet function for the EEC in the back-to-back limit

The EEC measures the energy correlation of two detectors in $e^{+}e^{-}$ annihilation at an angle $\chi$ . The TMDFFs obtained in the last subsection can be used to calculate the jet function for the EEC in the back-to-back limit $\chi\to\pi$ . It has been known for a long time that resummation of large logarithms for the EEC in this limit is closely related to $q_{T}$ resummation in the Drell-Yan process Collins:1981uk ; Dokshitzer:1999sh ; deFlorian:2004mp . Recently, an all-order factorization formula in terms of operator matrix elements for the EEC in the back-to-back limit has been presented Moult:2018jzp . The factorization formula at leading power reads

[TABLE]

where $z=(1-\cos\chi)/2$ . In the back-to-back limit one has $z\to 1$ . The hard function and soft function in Eq. (96) is well known and can be found to two loops in Ref. Moult:2018jzp . The only missing ingredient for the resummation at N3LL accuracy is the two-loop jet function. In QCD, due to charge conjugation invariance, we have $J^{q}(b_{\perp},\mu,\nu)=J^{\bar{q}}(b_{\perp},\mu,\nu)$ . The quark jet function can be obtained from the second Mellin moments of the matching coefficients of quark TMDFFs Moult:2018jzp

[TABLE]

We expand the jet function in terms of $\alpha_{s}$ as

[TABLE]

Using the two-loop TMDFFs computed in this paper, the expansion coefficients are then given by

[TABLE]

The $\mu$ and $\nu$ dependence of the jet function are in full agreement with the RGE and rapidity evolution equation Moult:2018jzp . The new result from this paper is the two-loop constant term

[TABLE]

which represents the last missing ingredient for N3LL resummation of EEC in the back-to-back limit.

Using the two-loop jet function together with the two-loop hard and soft function, we obtain the full leading power prediction for the EEC in the back-to-back limit through two loops from the factorization formula in Eq. (96), including the $\delta(1-z)$ terms,

[TABLE]

where we have set $\mu=Q$ and $T_{F}=1/2$ for simplicity, and

[TABLE]

The two-loop $\delta(1-z)$ term is

[TABLE]

The two-loop plus distribution terms $D_{n}(z)$ are in full agreement with the analytical NLO calculation in Ref. Dixon:2018qgp , while the two-loop $\delta(1-z)$ term is new. Eq. (4.2) has already been used in a previous publication to extract the $\delta(z)$ term of EEC using the energy conservation sum rule Dixon:2019uzg . Two independent checks are made for Eq. (4.2). Firstly, EEC in the back-to-back limit obey the leading transcendental principle Korchemsky:2019nzm ; Henn:2019gkr , which states that the maximal transcendental part of the QCD results are identical to the same quantity in ${\cal N}=4$ supersymmetric Yang-Mills (SYM) theory, up to trivial overall color factor Kotikov:2002ab . In Eq. (4.2), the leading transcendental term is the $\zeta_{4}$ terms. To compare with the same quantity in ${\cal N}=4$ SYM theory, we replace $C_{F}\to C_{A}$ in Eq. (4.2) and found the leading transcendental piece to be $40C_{A}^{2}\zeta_{4}$ , which is in full agreement with an independent calculation in Korchemsky:2019nzm . Secondly, Besides the energy conservation sum rule, EEC in massless perturbation theory also obey a sum rule due to momentum conservation, which reads Kologlu:2019mfz ; Korchemsky:2019nzm

[TABLE]

where $\sigma_{\rm tot}$ is the total hadronic cross section for $e^{+}e^{-}$ including higher order QCD corrections. Using the analytical NLO formula of EEC for $0<z<1$ from Dixon:2018qgp , and the end point contribution in Eq. (101), we explicit verify the sum rule in Eq. (104).

The end point contributions in Eq. (101) are directly computed using the two-loop jet function in Eq. (99), which by itself reduces to moment of the TMDFFs. Therefore, the checks made for the end point contributions apply also to the TMDFFs computed in this paper, in particular to its $\delta(1-z)$ terms.

5 Conclusion

In this work, we have revisited the calculation of perturbative quark TMDPDFs and TMDFFs at NNLO using a new regulator for rapidity divergences. We use the SIDIS process to set-up our calculation, while our results are universal and can be used for other processes as well. We show that the exponential regulator provides a consistent framework to carry out the calculation of the TMD soft functions, TMDPDFs and TMDFFs.

Compared to existing regulators in the literature, the exponential regulator has a couple of advantages. Firstly, the regulator can be implemented at the level of operator definitions for the TMD functions, where it manifests itself as a small shift of the space-time coordinates. Secondly, the exponential regulator is applied to the total momentum of the extra emissions in the final state. Except for this last integration, the regulator does not change the structure of (cut)-propagators in the amplitudes. As a result, we can apply many modern techniques for loop integrals such as IBP identities and differential equations. This allows us to obtain the bare NNLO TMDPDFs and TMDFFs up to $\mathcal{O}(\epsilon^{2})$ in dimensional regularization, and can also be extended to a future N3LO calculation. Finally, the regulator can already be expanded in terms of delta-functions and plus-distributions at the integrand level, which makes the final round of integration easy to carry out.

Our results for the quark TMDPDFs up to $\mathcal{O}(\epsilon^{0})$ agree with the results in the literature, while our results for the quark TMDFFs have a small discrepancy with another calculation. To further check our results, we use the TMDFFs to calculate the NNLO jet function appearing in the factorization formula of EEC in the back-to-back-limit. This also serves as a new result of our paper, and is the last missing ingredient for an N3LL resummation. We have checked that our NNLO jet function produces the correct leading singular terms for the EEC in the back-to-back limit. This is a strong validation of our results for the TMDFFs.

Given the benefits provided by the exponential regulator, the calculation for the gluon TMDPDFs and TMDFFs (which was more difficult than the quark case using other regulators) can also be greatly simplified. We also believe that our method can be extended to the N3LO level. We leave these considerations to future publications.

Acknowledgements.

This work was supported in part by the National Natural Science Foundation of China (11135006, 11275168, 11422544, 11375151, 11535002, 11575004, 11635001), and the Zhejiang University Fundamental Research Funds for the Central Universities (2017QNA3007).

Appendix A Anomalous dimensions, splitting functions and the TMD soft function

In this Appendix, we list some necessary ingredients which enter our calculation.

A.1 Anomalous dimensions

For all the anomalous dimensions entering the RGEs of various TMD functions, we define the perturbative expansion according to

[TABLE]

where the coefficients up to $\mathcal{O}(\alpha_{s}^{2})$ are

[TABLE]

The cusp anomalous dimension $\Gamma^{\text{cusp}}$ can be found in Korchemsky:1987wg ; Korchemskaya:1992je . The hard and soft anomalous dimensions $\gamma^{H}$ and $\gamma^{S}$ can be extracted from the two-loop quark form factor Gehrmann:2005pd ; Moch:2005id , and can also be found in, e.g., Refs. Becher:2009qa ; Li:2014afw . Finally, the beam anomalous dimension $\gamma^{B}$ is related to $\gamma^{S}$ and $\gamma^{H}$ through $\gamma^{B}=\gamma^{S}-\gamma^{H}$ . And the renormalization factor for the quark TMDPDFs and TMDFFs up to $\mathcal{O}(\alpha_{s}^{2})$ reads

[TABLE]

The QCD beta function is defined by

[TABLE]

with Gross:1973id ; Politzer:1973fx ; Caswell:1974gg ; Jones:1974mm ; Egorian:1978zx

[TABLE]

A formula particularly useful for us is

[TABLE]

where

[TABLE]

A.2 Space-like splitting functions

The first order space-like splitting functions can be written as Altarelli:1977zs

[TABLE]

where

[TABLE]

The second order space-like splitting functions are Furmanski:1980cm ; Curci:1980uw

[TABLE]

A.3 Time-like splitting function

The first order time-like splitting functions are exactly the same as the space-like ones, while the second order time-like splitting functions are given by Furmanski:1980cm ; Curci:1980uw

[TABLE]

where we use the same convention as in Ref. Ritzmann:2014mka .

A.4 TMD soft function

The exponentially regularized TMD soft function is given by Li:2016ctv

[TABLE]

where $L_{R}=L_{\perp}+L_{\nu}$ with $L_{\nu}=\ln(\nu^{2}/\mu^{2})$ and the scale-independent terms are

[TABLE]

Bibliography102

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1(1) J. Collins, Foundations of perturbative QCD , Camb. Monogr. Part. Phys. Nucl. Phys. Cosmol. 32 (2011) 1–624.
2(2) R. Angeles-Martinez et al., Transverse Momentum Dependent (TMD) parton distribution functions: status and prospects , Acta Phys. Polon. B 46 (2015), no. 12 2501–2534, [ ar Xiv:1507.05267 ].
3(3) Y. L. Dokshitzer, D. Diakonov, and S. I. Troian, On the Transverse Momentum Distribution of Massive Lepton Pairs , Phys. Lett. 79B (1978) 269–272.
4(4) G. Parisi and R. Petronzio, Small Transverse Momentum Distributions in Hard Processes , Nucl. Phys. B 154 (1979) 427–440.
5(5) J. C. Collins, D. E. Soper, and G. F. Sterman, Transverse Momentum Distribution in Drell-Yan Pair and W and Z Boson Production , Nucl. Phys. B 250 (1985) 199–224.
6(6) P. B. Arnold and R. P. Kauffman, W and Z production at next-to-leading order: From large q(t) to small , Nucl. Phys. B 349 (1991) 381–413.
7(7) G. A. Ladinsky and C. P. Yuan, The Nonperturbative regime in QCD resummation for gauge boson production at hadron colliders , Phys. Rev. D 50 (1994) R 4239, [ hep-ph/9311341 ].
8(8) G. Bozzi, S. Catani, G. Ferrera, D. de Florian, and M. Grazzini, Production of Drell-Yan lepton pairs in hadron collisions: Transverse-momentum resummation at next-to-next-to-leading logarithmic accuracy , Phys. Lett. B 696 (2011) 207–213, [ ar Xiv:1007.2351 ].