QCD Corrections to Pair Production of Type III Seesaw Leptons at Hadron Colliders

If kinematically accessible, hadron collider experiments provide an ideal laboratory for the direct production of heavy lepton partners in Seesaw models. In the context of the Type III Seesaw Mechanism, the $\mathcal{O}(\alpha_s)$ rate and shape corrections are presented for the pair production of hypothetical, heavy $SU(2)_L$ triplet leptons in $pp$ collisions at $\sqrt{s}=$ 13, 14, and 100 TeV. The next-to-leading order (NLO) $K$-factors span, approximately, $K^{NLO}=1.1 - 1.4$ for both charged current and neutral current processes over a triplet mass range $m_T = 100~\text{GeV}-2~\text{TeV}$. Total production cross sections exhibit a $^{+5\%}_{-6\%}$ scale dependence at 14 TeV and $\pm1\%$ at 100 TeV. The NLO differential $K$-factors for heavy lepton kinematics are largely flat, suggesting that na\"ive scaling by the total $K^{NLO}$ is reasonably justified. The resummed transverse momentum distribution of the dilepton system is presented at leading logarithmic (LL) accuracy. The effects of resummation are large in TeV-scale dilepton systems. Discovery potential to heavy lepton pairs at 14 and 100 TeV is briefly explored: At the High-Luminosity LHC, we estimate a $4.8-6.3\sigma$ discovery potential maximally for $m_T = 1.5-1.6~\text{TeV}$ after 3000 fb$^{-1}$. With 300 (3000) fb$^{-1}$, there is $2\sigma$ sensitivity up to $m_T = 1.3-1.4~\text{TeV}~(1.7-1.8~\text{TeV})$ in the individual channels. At 100 TeV and with 10 fb$^{-1}$, a $5\sigma$ discovery can be achieved for $m_T=1.4-1.6~\text{TeV}$. Due to the factorization properties of Drell-Yan-type systems, the fixed order and resummed calculations reduce to convolutions over tree-level quantities.


Introduction
The origin of sub-eV neutrino masses is a central issue in particle physics. As right-handed neutrinos do not exist in the Standard Model (SM), which thus predicts massless neutrinos, new particles are necessary to explain neutrino masses [1], e.g., gauge singlet fermions in the Type I [2][3][4][5][6][7] Seesaw Mechanism, or scalar and fermionic SU(2) L triplets in the Types II [8][9][10][11] and III [12] scenarios. Searches for these degrees of freedom constitute an important component of hadron collider programs; see Refs. [13][14][15][16] and references therein. Furthermore, the maturity of the formalism underlying QCD corrections in hadron collisions, which are required for predicting accurate production rates and distribution shapes, readily permit their application to beyond the SM (BSM) processes.
In the Type I Seesaw, production cross section is known at next-to-leading order (NLO) in QCD [26] and estimated at next-to-next-to-leading order (NNLO) via a K-factor * [25]. For the Type II case, rates are known at NLO [27,28], NLO with next-to-leading logarithm (NLL) recoil and threshold resummations [29], and automated at NNLO [30,31].
Pair production of heavy Type III Seesaw leptons has, until now, been evaluated only to leading order (LO) accuracy. For m T = 100 GeV − 2 TeV, we report the NLO K-factors: with scale uncertainty of +5% −6% at 14 TeV and ±1% at 100 TeV, and are comparable to other DY-type processes in Seesaw models. The NLO differential K-factors † for heavy lepton kinematics are largely flat for TeV-scale m T , suggesting that naïve scaling by the total K N LO is reasonably justified.
In this study, production rates of TeV-scale Type III Seesaw lepton pairs at O(α s ) accuracy are presented for pp collisions at √ s = 13, 14, and 100 TeV. Differential distributions at NLO and NLO with leading logarithm (LL) resummation of TeV-scale lepton kinematics are presented for the first time at 14 TeV. The fixed order (FO) calculation is carried out via phase space slicing (PSS) [33][34][35][36]. The calculation of the dilepton system's transverse momentum, q T , follows the Collins-Soper-Sterman (CSS) formalism [37][38][39]. This text continues in the following order: In section 2, we summarize the Type III Seesaw model and comment on experimental constraints. The PSS and CSS formalisms are briefly introduced in section 3. Due to the factorization properties of DY-type systems, the fixed order and resummed results reduce to convolutions over tree-level quantities; technical details are relegated to appendices A and B. Results are reported in section 4. We summarize and conclude in section 5.
2 Type III Seesaw Mechanism

Model Lagrangian
The Type III Seesaw [12] generates tree-level neutrino masses via couplings to SU(2) L triplet leptons with zero hypercharge. In terms of Pauli matrices σ a , the left-handed (LH) fields are denoted by where Σ ± L have U(1) EM charges Q = ±1, and the right-handed (RH) conjugate fields are (2.2) Chiral conjugates are related by ψ c R ≡ (ψ c ) R = (ψ L ) c , where ψ L/R ≡ P L/R ψ = 1 2 (1 ∓ γ 5 )ψ. For a single generation (but generalizable to more), the model's Lagrangian is where L SM is the SM Lagrangian, the triplet's covariant derivative and mass are given by 4) and the SM LH lepton (L) and Higgs (Φ) doublet fields couple to Σ L via the Yukawa coupling Dirac masses are then spontaneously generated after electroweak symmetry breaking (EWSB): (2.6) leading to the neutral fermion mass matrix The Seesaw mechanism proceeds by supposing m T ≫ y T Φ , leading to light/heavy mass eigenvalues Thus, tiny neutrino masses follow from mixing with heavy states, whereby light (heavy) mass eigenstates align with the doublet (triplet) gauge states. For y T comparable to the electron's SM Yukawa, sub-eV m light can be explained by sub-TeV m T , a scale within the LHC's kinematic reach. We combine the Σ fields and their conjugates into physical Dirac and Majorana fields: In the gauge basis and in terms ofT , the triplet interaction Lagrangian L T is written as Our aim is to report the O(α s ) corrections to heavy lepton pair production, which are independent of the mixing between the gauge statesT and mass states T, ℓ (ℓ = e, µ, τ ), and ν m (m = 1, 2, 3). For the remainder of the study, we generically denote the mixing as Y and ε, and writẽ The resulting interaction Lagrangian in the mass eigenbasis relevant to our study is
• Collider Production and Decay: CMS experiment searches for T 0 T ± production and decay into ℓ + ℓ − ℓ ′ ± E T have restricted the cross section and branching ratio to [46] [47], depending on mixing parameters, m T < 325 − 540 GeV at the 95% C.L. (2.15) Throughout this study, we take T 0 and T ± to be mass degenerate. Electroweak (EW) corrections at one loop induce a mass splitting of ∆m T ≈ 160 MeV for m T > 100 GeV [42,48,49], and is thus negligible. For differential distributions, we use representative mass m T = 500 GeV. (2.16) As in the LO case, the total partonic and hadronic cross sections at NLO and NLO+LL factorize into a product of the mixing parameter |Y | and a mixing-independent "bare" cross section σ 0 : Therefore, we express our results in terms of σ 0 and do not choose any particular |Y |. Furthermore, factorization implies that that total and differential NLO K-factors are independent of |Y |.

Heavy Lepton Pair Production at O(α s ) in Hadron Collisions
Here we outline the PSS [33][34][35][36] and CSS [37][38][39] formalisms, which we use to calculate the processes We note that these corrections are not unique but are well-known and general for the production of any SU (2) L triplet color-singlet, e.g., [50,61]. However, unlike previous studies, we investigate the O(α s ) effects on the kinematic distributions of TeV-scale leptons.

Phase Space Slicing
To evaluate T T production at NLO, we follow the usual procedure: evaluate virtual and radiative corrections to the LO process in d = 4 − 2ε dimensions; collect soft divergences, which cancel exactly; collect collinear divergences, which cancel partially; and subtract residual collinear poles from parton distribution functions (PDFs). For an n-body LO process, we divide, or slice, the phase space of its (n + 1)-body correction into soft and collinear kinematic regions. For radiation energy E j , partonic c.m. energy √ŝ , and small dimensionless cutoff parameters δ S , δ C , a volume of the (n + 1)body phase space is soft if For partonic-level invariant masses and momentum transferŝ where indices i, k run over initial-and final-state momenta, a region of phase space is collinear ifŝ A volume is hard (non-collinear ) if not soft (collinear). Exact choices of δ S , δ C do not matter: dependences on δ S , δ C cancel for sufficiently inclusive processes [36]. However, so soft and collinear factorization remain justified, one needs The hard-non-collinear T T j process is then finite everywhere and given by whereσ B is the Born-level T T j partonic cross section. For a, b ∈ {q, q ′ , g} with q ∈ {u, d, c, s}, the PDF f a/p (ξ i , µ 2 f ) is the likelihood of parton a carrying away longitudinal momentum fraction ξ i from proton p evolved to a factorization scale µ f . The c.m. beam energy √ s and partonic c.m. energy are related byŝ = ξ 1 ξ 2 s, and we denote the threshold at which T T production occurs by τ 0 : In the soft/collinear limits, amplitudes for soft, soft-collinear, and hard-collinear radiation factorize into divergent expressions proportional to the (color-connected) Born amplitude. The poles are grouped with virtual corrections and the PDFs, resulting in a finite expression given by [36] σ (2) wereσ B is the Born-level T T partonic cross section,σ V is its O(α s ) virtual correction, σ S andσ SC are the soft and soft-collinear radiation terms, and σ HC is the hard-collinear radiation correction. Inclusive triplet lepton production at NLO is now reduced to a sum of two-and three-body processes:   (2) for δ S O(0.5) is due to large hard-collinear PDF subtractions. The NLO result is insensitive to δ S for δ S 1−3×10 −3 , reflecting the large but fine cancellation of Eqs. (3.7) and (3.9). For example: for δ S = 1×10 −5 , the three-and two-body calculations are approximately +10.2 and −9.0 times the LO cross section. For δ S 3 × 10 −3 , the Kfactor plummets, indicating the importance of terms proportional to powers of δ S in the two-body expression, and hence a breakdown of soft/collinear factorization.

Collins-Soper-Sterman Transverse Momentum Resummation
As a color-singlet process, colored initial state radiation (ISR) is the dominant contribution to the T T system's q T spectrum. For dilepton invariant mass M V * , the distribution has the power series and indicates a breakdown of the perturbative description in the q 2 That is, α s (M 2 V * ) must be perturbative and the scales associated with the process must be comparable. As the Seesaw mass scale is pushed higher [46,47], so too does the scale at which artificially large logarithms appear. For m T = 500 GeV (1 TeV), FO predictions breakdown at q T ∼ 55 (95) GeV, and resummation of recoil logarithms become necessary to describe q T below this threshold.
Fortunately, as q T /M V * → 0, gluon radiation factorizes. This permits one to reorganize, sum, and exponentiate large logarithms in Eq. (3.12), resulting in an all-orders expression in terms of the Born process. The resummed distribution with respect to q 2 T , M 2 V * , and dilepton rapidity y, is [39] dσ 14) The integral is over the impact parameter b and is the Fourier transform of q T ; the zeroth order Bessel function J 0 emerges as a simplification. W expands to a perturbative and nonperturbative set of universal Sudakov form factors, and a process-dependent luminosity weightW : Expressions for S N P and S P are in appendix B, in Eqs. (B.10) and (B.14). The resummed result describes well the q T ≪ M V * behavior due to the Sudakov suppression. But because of neglected terms proportional to powers of q T , it underestimates the spectrum at q T M V * , precisely where the FO calculation becomes reliable. To describe accurately q T everywhere, one introduces the auxiliary function dσ Asymp that matches the asymptotic FO (resummed) behavior at small (large) q T . Combining the three expressions, the total, matched spectrum is given by [51,52] dσ The area bound by the dσ Matched curve is then normalized to the total σ N LO rate of Eq. (3.11) [67]. Individual terms of Eq.

Results
Tree-level results are calculated using helicity amplitudes. The Cuba library [53]  The Cabbibo-Kobayashi-Masakawa (CKM) matrix is taken to be unity, introducing a percent-level error that is no larger than the estimated O(α 2 s ) contributions. The CJ12mid NLO parton distribution functions (PDFs) [59] with α s (M Z ) = 0.118 are used. Since 2m T > M W , M Z , the factorization (µ f ) and renormalization (µ r ) scales are fixed to the sum of the triplet lepton masses . Setting n f = 6 increases the threebody channel by 2 − 3%, but the total NLO cross section by less than +1%. For total cross section calculations, we choose soft and collinear cutoffs PSS involves fine cancellation of large numbers; for differential distributions, δ S is relaxed to Born-level T T jX events can be generated efficiently by implementing soft and collinear cuts Eqs.   The CC rate is approximately twice as large as the NC channel due to W boson charge multiplicity. In the low-(high-)mass range, transitioning from 14 to 100 TeV increases the total cross section roughly by a factor of 10 (1000) for both processes.  At lower collider energies, K-factors are larger for heavier m T due to the rarity of antiquarks possessing sufficiently large momentum at LO. At NLO, this is compensated by large Bjorken-x gluons undergoing high-p T g → q splitting. CC and NC K-factors are appreciable and, due to their color structures, comparable to those of the Seesaw Types I [25] and II [27][28][29][30][31]. Table 1 summarizes these results for representative m T at √ s = 13, 14, and 100 TeV pp collider configurations.

Scale Dependence
Higher order QCD corrections are necessary to further reduce theoretical uncertainty. To quantify and estimate the size of these contributions, the default scale µ 0 = 2m T is varied over the range The scale variation is then defined as the ratio of the NLO rate evaluated at scale µ to the same rate at µ = µ 0 . Figures 4(a) and 4(b) show, respectively, the scale variation band of   the CC and NC processes at NLO as a function of m T ; the lower panels show the NLO K-factor for the three scale choices. The default scale is denoted by a solid line at 1. The high (low) scale scheme is denoted by right-side up (upside down) triangles and is found to decrease (increase) the total cross section. This suggests that the renormalization scale evolves α s (µ 2 ) to smaller values faster than the factorization scale evolves PDFs to larger values, and also leads (accidentally) to a vanishing scale dependence for m T ∼ 100 GeV.
The behavior is consistent with other TeV-scale Seesaws mechanisms [25]. For the m T studied, the CC and NC calculations exhibit maximally a +5% −6% scale dependence at 14 TeV; this reduces to ±1% at 100 TeV for the same mass range considered. The 14 TeV NLO K-factor varies maximally +7% −8% for the CC process and +8% −8% for the NC process. The slightly larger scale dependence in the K-factors than the total cross sections is due to the larger scale dependence of the LO result. The size of the µ-dependence suggests that O(α 2 s ) effects are small, consistent with Refs. [25,30,31]. Results are summarized in table 2 for representative m T . Figure 5(a) shows the 14 TeV NLO differential distribution, divided by the mixing parameter |Y | 2 , with respect to the p T of T ± in T 0 T ± production. The panel shows the differential NLO K-factors. At low (high) p T the change at NLO is small (large) and follows from the T 0 T ± j channel at O(α s ). The transverse recoil of dilepton system from hard ISR propagates to individual leptons, thereby providing an additional transverse boost. As the jet energy softens or is radiated more collinearly to its progenitor, its p T vanishes, and kinematics at NLO approach those at LO. In figure 5(b), the rapidity (y) distribution of T ± in T 0 T ± production is presented. We observe that QCD corrections have a small impact on the rapidity distribution shape as indicated by the mostly flat K-factor, with only a slight upwards bump at small y, where p T ≫ p z . Similar to the p T spectrum, new kinematic channels at O(α s ) all involve high-p T ISR and do not induce longitudinal boosts, leaving the y distribution shape largely unchanged.

14 TeV Kinematic Distributions at NLO and NLO+LL
Similar p T and y behavior are observed for T 0 in T 0 T ± and T ± in the NC processes. Aside from the Majorana nature of light (ν) and heavy (T 0 ) neutral leptons, and the relative CC and NC production rates, i.e., σ(pp → T 0 T ± )/σ(pp → T + T − ) ∼ 2, observing vector-like coupling of heavy leptons to electroweak gauge bosons is a critical test of the Type III Seesaw mechanism [42]. As in the SM, this done by measuring the polar distribution made by, for example, T ± in T 0 T ± production in the dilepton rest frame with respect to dilepton system's direction of propagation in the lab frame. Symbolically, the observable is given by where p * T is the 3-momentum of lepton T in the T T frame and q = p T + p T in the lab frame. Figures 6(a) and 6(b) show, respectively, the CC and NC cos θ * distribution. At LO and NLO, the vector coupling structure is clear. The uniform K-factor follows from O(α s ) corrections involving only initial-state partons and amount simply to a boost of the dilepton system in the lab frame. The O(α s ) effects are unraveled in constructing Eq. (4.14) and subsequently affect only the normalization.
As discussed in section 3.2, the transverse momentum distribution q T of the T ± T 0 dilepton system is ill-defined at q T ≪ M V * in FO perturbation theory and requires resum-  Figure 6. Same as figure 5 but with respect to cos θ * , as defined in Eq. (4.14), of (a) T ± in T 0 T ± production and (b) T + in T + T − production.
mation of logarithms associated with soft gluon radiation. For representative masses (a) m T = 500 GeV and (b) 1 TeV, figure 7 shows the various contributions to q T spectrum: the asymptotic (dash) and FO (solid) terms at O(α s ), the resummed rate at LL (dot), and combination of the pieces. The lower panel shows the ratio of the combined result to the FO result. For m T = 500 GeV (1 TeV), the FO calculation overestimates the combined differential rate for q T 25 GeV. At q T ∼ 55 (95) GeV, the estimate from Eq. (3.13), the FO result remains about 35% (70%) below the combined result. The largeness of the resummation corrections is consistent with other recoil resummations at high-scales [67]. The combined distributions peak at q T ≈ 5 − 6 GeV, below which the Sudakov suppression from multiple soft gluon emissions overtakes the divergent nature of soft emissions.

Discovery Potential at 14 TeV High-Luminosity LHC and 100 TeV
In the high luminosity LHC (HL-LHC) scenario [60], detector experiments aim to each collect 1−3 ab −1 of data. We briefly address the maximum sensitivity to heavy triplet lepton pair production in this scenario. TeV-scale triple leptons decay dominantly to longitudinally polarized weak bosons and the Higgs [42,44], a consequence of the Goldstone Equivalence Theorem, implying For visible decay modes Z → jj and h → bb, gg, heavy lepton pairs can decay into fully reconstructible final-states with four jets and two high-p T leptons that scale like p ℓ T ∼ m T /2: The corresponding branching fractions are BR T 0 T ± → ℓℓ ′ + 4j/2j + 2b ≈ 11.5%, (4.19) BR (T + T − → ℓℓ ′ + 4j/2j + 2b/4b) ≈ 11.6%. (4.20) Taking |Y | 2 = 1, m T = 1.5 TeV, and acceptance efficiency of A = 0.75 [42], then after 3 ab −1 one expects tens of heavy lepton pairs across both channels To a good approximation, the kinematics of TeV-scale T T decays render the SM background negligible [42,44]. Using a Gaussian estimator, the statistical significances are at the 3 − 5σ level: Summing in quadrature, the combined significance surpasses the 6σ level (over the null hypothesis): For m T = 1.6 TeV, the combined significance is approximately 4.8σ, demonstrating a maximum HL-LHC discovery potential to Type III Seesaw leptons in the m T = 1.5 − 1.6 TeV range.
Fixing the branching fractions and acceptance rates, we plot in figure 8 the required luminosity as a function of m T for a 5σ discovery (dash-star) and 2σ sensitivity (dashdiamond) of the (a) T 0 T ± and (b) T + T − channels at 14 and 100 TeV. At 14 TeV and after 300 (3000) fb −1 , one finds sensitivity to triplet pairs up to m T = 1.3 − 1.4 TeV (1.7 − 1.8 TeV) in the individual CC and NC channels. At 100 TeV, we observe that with 10 fb −1 a 5σ discovery can be achieved in the CC channel for m T ≈ 1.6 TeV and in the NC channel for m T ≈ 1.4 TeV. At large m T , however, taking A = 0.75 is a not justified as a different search methodology is required to account for the boosted kinematics of the T T decay products.

Summary and Conclusions
The existence of tiny neutrino masses has broad impact in cosmology, particle, and nuclear physics. Hence, understanding the origin of their sub-eV masses is a pressing issue. We report the leading QCD corrections to the production rate and kinematic distributions of hypothetical, heavy Type III Seesaw lepton pairs, including soft gluon resummation of the heavy dilepton system. We find: 1. The pp → T 0 T ± and T + T − K-factors at NLO in QCD (see 2. For the range of m T considered, the total T 0 T ± and T + T − production cross sections exhibit a +5% −6% scale dependence at 14 TeV and ±1% dependence at 100 TeV; see section 4.2. 3. Differential K-factors with respect to rapidity, transverse momentum, and polar distributions of heavy triplet leptons are largely flat, indicating that naïve rescaling of Born-level result is a largely justified estimate of kinematics at O(α s ); see section 4.3.
4. The resummed transverse momentum of the T ± T 0 dilepton system illustrates that the impact of TeV-scale systems recoiling off soft radiation is large; see figure 7.

A Triplet Lepton Pair Production at NLO in QCD via PSS
The PSS technique and explicit examples, including the DY process, are reviewed authoritatively in Ref. [36]. Terms relevant to triplet lepton pair production in hadron collisions are collected here for convenience. We emphasize that the NLO calculation greatly simplifies to simply convolving PDFs with tree-level amplitudes.

A.1 QCD-Corrected Two-Body Final State
In the soft/collinear limits, radiation in T T j production becomes unresolvable and the three-body kinematics approach those of the two-body process Subsequently, the singular propagators in the amplitude factor and the soft/collinear contributions to the T T j cross section can be expressed as a divergent, but process-independent, piece and the LO T T cross section. Loop corrections cancel soft and soft-collinear poles rendering finite the quantity PDFs are redefined to subtract residual hard-collinear divergence. We now give each term explicitly.

A.1.1 Virtual Corrections
For massless fermions, O(α s ) corrections to EW vertices, i.e., figure 1 (b), factorize [61] into a product of the Born matrix element, M B , and a universal form factor: The vertex correction is UV-finite, so its counter term is zero. Counter terms cancel external quark self-energy corrections in the on-shell and MS schemes. Hence, the summed, squared amplitude is The Born and virtual terms at O(α s ) for heavy lepton production are then (A.10)

A.1.2 Soft Corrections
For colored partons a, b, and color-connected Born cross sectionσ 0 ab , the soft radiation expression is where dS is the soft one-particle phase space in n = 4 − 2ε dimensions, and we make use of the fact that the only two colored partons in the DY process are the initial-state quark and antiquark. For the present case, the color-connected and plain Born cross sections are related by We take p q (p q ′ ) to propagate in the +ẑ (−ẑ) direction, implying Using the tabled integrals of Ref. [36], the soft contribution to DY processes is

A.1.3 Soft Collinear Corrections
For qq ′ initial-states, the soft-collinear correction to the Born amplitude is The soft-collinear splitting functions, A(1 → 2 3), are derived by integrating over the Altarelli-Parisi (AP) splitting functions. The q → q and q ′ → q ′ functions are equal by CP-invariance and therefore are combined in the last step.

A.1.4 Assembly of Two-Body Corrections
A finite result emerges after summing the virtual, soft, and soft-collinear corrections, With poles removed, we send ε → 0. The non-hard-collinear, two-body contribution of Eq. (3.10) is

A.2 Hard-Collinear Subtraction
Remaining hard-collinear ISR poles are associated with PDFs, which themselves are all orders resummations of hard-collinear splittings. The divergences are removed by subtracting out from PDFs redundant O(α s ) splittings. The hard-collinear subtraction is given by where the redefined, O(α s )-corrected PDFsf a/p (ξ, µ 2 ) arẽ The summation is over all partons that give rise to q, q ′ after one radiation, δ ba is the Kronecker δ-function, and the modified AP i → j splitting functionsP ji (z) have the form In n = 4 − 2ε dimensions, P and P ′ are related by P ji (z, ε) = P ji (z) + εP ′ ji (z), with

A.3 Efficient Event Generation of Three-Body, Hard, Non-Collinear Radiation
Imposing soft/collinear cuts in the partonic c.m. frame on the processes regulates all diverges. T T j event generation can be handled efficiently by implementing Eqs.
For the processes q q ′ → V * → T T , V ∈ {γ, Z, W }, (B.1) the resummed and matched q T distribution of the T T system decomposes into the FO, resummed, and asymptotic pieces; see Eq. (3.17). The oscillatory behavior of J 0 (q T b) and the fine numerical cancellation of the resummed and asymptotic terms when q T M V * give rise to an unnecessary computational burden for a contribution that is small compared to the FO piece. In practice [62,63], one introduces f Match (q T ) that is unity in the q T /M V * → 0 limit and vanishes the asymptotic-resummed difference at q T M V * : We set q Match [52], and write the matched, triply differential distribution [51,52] dσ Matched dM 2 V * dy dq 2 In the q T /M V * → 0 limit, the longitudinal momentum fractions of initial-state partons q, q ′ are (B.4) ‡ We report and clarify a slight but unfortunate typo in Ref. [52]: the resummed result is to LL accuracy, not NLL.
Though untrue for the resummed and asymptotic pieces at large q T /M V * , their combined contribution vanishes by construction. For the FO calculation in this limit, M V * → √ŝ and the above momentum fractions are for initial-state partons q, q ′ , and g. With this, one may write Defining dσ ∆ /dξ 1 dξ 2 dq T as the quantity in the parentheses, one at last has To account for the appropriate normalization at O(α s ), the area under the above distribution is scaled to equal σ N LO in Eq. (3.11). The resummed and asymptotic terms are now discussed.

B.1 Recoil q T Resummation
Recalling the resummed expression, but in terms of ξ i , . This is ameliorated [39,64] by S N P , which is valid for all b and approaches unity when b ≪ 1/Λ QCD , and by introducing This approach, however, requires that S N P be extracted from data and dependent on b max . We choose the non-perturbative Sudakov factor given by the BLNY parameterization [65,66] S N P (b, Q, ξ 1 , ξ 2 ) = b 2 g 1 + g 2 log Q 2Q BLNY + g 1 g 3 log(100ξ 1 ξ 2 ) . To numerically integrate over the domain of b, we exploit the b-dependence of S N P and define: The b integral now takes the manageable form The perturbative Sudakov-like form factor is given by where the low-scale integration limit is µ = c 0 /b * , with c 0 = 2e −γ E and γ E ≈ 0.577 is the Euler-Mascheroni constant. The functions A and B can be expanded perturbatively in powers of α s (µ 2 ). For a resummed calculation at LL, we take A and B to O(α n=1 s ). The expressions are [39] At one-loop, α s is given by and allows Eq. (B.14) to be evaluated analytically [67]. The result is , t(µ 2 ) = log (B.20) In the perturbative limit b ≪ 1/Λ QCD , the TMD F T factorize into universal Wilson coefficient functions for i → j splitting, C ji , and the usual transverse momentum-independent PDFs: (B.21) The summation is over all partons i, k that can split into q, q ′ . The convolution notation denotes Formally, C ji (z) (1) is an O(α 2 s ) contribution except at q T = 0, where it O(α s ). Hence, our q T spectrum is accurate to O(α s ) in shape everywhere but the origin [52], an acceptable error considering the uncertainty in our knowledge of m T . Evaluating the convolutions, the product of TMDs is qq ⊗ f q/p = f q/p (ξ 1 , µ 2 f )f q/p (ξ 2 , µ 2 f ).