Generalized Threshold Factorization with Full Collinear Dynamics

TL;DR

This paper introduces a generalized threshold factorization theorem for color-singlet processes, capturing soft and collinear effects more comprehensively than traditional methods, with applications to high-order calculations and PDF fits.

Contribution

It develops a new factorization theorem valid in a broader kinematic limit, including flavor-nondiagonal channels, and demonstrates its application to NNLO and N3LO predictions for rapidity spectra.

Findings

Explicit NNLO calculations for Z and Higgs rapidity spectra.

Predictions for N3LO contributions to rapidity spectra.

Identification of the correct soft threshold factorization among existing literature.

Abstract

Soft threshold factorization has been used extensively to study hadronic collisions. It is derived in the limit where the momentum fractions $x_{a,b}$ of both incoming partons approach $x_{a,b}\to 1$. We present a generalized threshold factorization theorem for color-singlet processes, which holds in the weaker limit of only $x_a \to 1$ for generic $x_b$ (or vice versa), corresponding to the limit of large rapidity but generic invariant mass of the produced color singlet. It encodes the complete soft and/or collinear singular structure in the partonic momentum fractions to all orders in perturbation theory, including in particular flavor-nondiagonal partonic channels at leading power. It provides a more powerful approximation than the classic soft threshold limit, capturing a much larger set of contributions. We demonstrate this explicitly for the Z and Higgs rapidity spectrum to NNLO, and we use it to predict a nontrivial set of its N3LO contributions. Our factorization theorem provides the relevant resummation of large-$x$ logarithms in the rapidity spectrum required for resummation-improved PDF fits. One of our factorization ingredients is a new beam function closely related to the N-jettiness beam function. As a byproduct, we identify the correct soft threshold factorization for rapidity spectra among the differing results in the literature.

Tables

Table 1. Table S1: Relevant EFT modes in the limit λ QCD ∼ λ 2 ∼ 1 − x − ≪ 1 similar-to subscript 𝜆 QCD superscript 𝜆 2 similar-to 1 superscript 𝑥 much-less-than 1 \lambda_{\mathrm{QCD}}\sim\lambda^{2}\sim 1-x^{-}\ll 1 in the lab (hadronic center-of-mass) frame and the leptonic frame where Y ^ = 0 ^ 𝑌 0 \hat{Y}=0 . In the right column we used that in the leptonic frame, q ± → q ^ ± = q + ​ q − ∼ Q → superscript 𝑞 plus-or-minus superscript ^ 𝑞 plus-or-minus superscript 𝑞 superscript 𝑞 similar-to 𝑄 q^{\pm}\to\hat{q}^{\pm}=\sqrt{q^{+}q^{-}}\sim Q .

Mode	Lab frame	Leptonic ($\hat{Y}=0$) frame
	$(+,-,⟂)$	$(+,-,⟂)$
$p_{\overline{n}}$	$(q^{+},{\lambda }^2{q}^{-},\lambda \sqrt{q^{-}{q}^{+}})$	$Q(1,{\lambda }^2,\lambda )$
$P_{\overline{n}}$	$(q^{+},\frac{{\Lambda }_{\mathrm{QCD}}^2}{q^{+}},{\Lambda }_{\mathrm{QCD}})$	$Q(1,{\lambda }_{\mathrm{QCD}}^2,{\lambda }_{\mathrm{QCD}})$
$P_n$	$(\frac{{\Lambda }_{\mathrm{QCD}}^2}{q^{-}},q^{-},{\Lambda }_{\mathrm{QCD}})$	$Q({\lambda }_{\mathrm{QCD}}^2,1,{\lambda }_{\mathrm{QCD}})$
$P_s$	$(\frac{1}{{\lambda }^2}\frac{{\Lambda }_{\mathrm{QCD}}^2}{q^{-}},{\lambda }^2{q}^{-},{\Lambda }_{\mathrm{QCD}})$	$Q(\frac{{\lambda }_{\mathrm{QCD}}^2}{{\lambda }^2},{\lambda }^2,{\lambda }_{\mathrm{QCD}})$
$P_{us}$	$(\frac{{\Lambda }_{\mathrm{QCD}}^2}{q^{-}},\frac{{\Lambda }_{\mathrm{QCD}}^2}{q^{+}},\frac{{\Lambda }_{\mathrm{QCD}}^2}{\sqrt{q^{+}{q}^{-}}})$	$Q({\lambda }_{\mathrm{QCD}}^2,{\lambda }_{\mathrm{QCD}}^2,{\lambda }_{\mathrm{QCD}}^2)$
$P_G$	$(\frac{{\Lambda }_{\mathrm{QCD}}^2}{q^{-}},\frac{{\Lambda }_{\mathrm{QCD}}^2}{q^{+}},{\Lambda }_{\mathrm{QCD}})$	$Q({\lambda }_{\mathrm{QCD}}^2,{\lambda }_{\mathrm{QCD}}^2,{\lambda }_{\mathrm{QCD}})$

Table 2. Table S2: Translation identities of two-dimensional plus distributions between the ( z , y ) 𝑧 𝑦 (z,y) and ( z a , z b ) subscript 𝑧 𝑎 subscript 𝑧 𝑏 (z_{a},z_{b}) parametrizations. Here, f ​ ( x ) 𝑓 𝑥 f(x) is an arbitrary function of x 𝑥 x , potentially distribution-valued for x → 1 → 𝑥 1 x\to 1 , while r ​ ( x ) = 𝒪 ​ [ ( 1 − x ) 0 ] 𝑟 𝑥 𝒪 delimited-[] superscript 1 𝑥 0 r(x)=\mathcal{O}[(1-x)^{0}] has at most an integrable singularity for x → 1 → 𝑥 1 x\to 1 . When indicated, the relations receive power corrections in 1 − z a 1 subscript 𝑧 𝑎 1-z_{a} starting at 𝒪 ​ ( 1 ) ≡ 𝒪 ​ [ ( 1 − z a ) 0 ] 𝒪 1 𝒪 delimited-[] superscript 1 subscript 𝑧 𝑎 0 \mathcal{O}(1)\equiv\mathcal{O}[(1-z_{a})^{0}] . Overall factors of θ ​ ( 1 − z ) ​ θ ​ ( y ) ​ θ ​ ( 1 − y ) = θ ​ ( 1 − z a ) ​ θ ​ ( 1 − z b ) 𝜃 1 𝑧 𝜃 𝑦 𝜃 1 𝑦 𝜃 1 subscript 𝑧 𝑎 𝜃 1 subscript 𝑧 𝑏 \theta(1-z)\,\theta(y)\,\theta(1-y)=\theta(1-z_{a})\,\theta(1-z_{b}) are understood on both sides.

$dzdy\text{LHS}$	$=d{z}_{a}d{z}_b\text{RHS}$
$f\left(z\right)\delta \left(y\right)$	$f\left(z_b\right)\delta \left(1-z_a\right)$
$r\left(z\right){ℒ}_0\left(y\right)$	$r\left(z_b\right)\left[{ℒ}_0\left(1-z_a\right)+\delta \left(1-z_a\right)\ln \frac{2z_b}{\left(1+z_b\right)\left(1-z_b\right)}\right]+𝒪\left(1\right)$
$f\left(z\right)\delta \left(1-y\right)$	$f\left(z_a\right)\delta \left(1-z_b\right)$
$r\left(z\right){ℒ}_0\left(1-y\right)$	$𝒪\left(1\right)$
$r\left(z\right)𝒪\left[y^0{\left(1-y\right)}^0\right]$	$𝒪\left(1\right)$
${ℒ}_0\left(1-z\right)\left[{ℒ}_0\left(y\right)+{ℒ}_0\left(1-y\right)\right]$	$-{ℒ}_1\left(1-z_a\right)\delta \left(1-z_b\right)+{ℒ}_0\left(1-z_a\right){ℒ}_0\left(1-z_b\right)-\delta \left(1-z_a\right){ℒ}_1\left(1-z_b\right)$
	$+\frac{{\pi }^2}{6}\delta \left(1-z_a\right)\delta \left(1-z_b\right)+\delta \left(1-z_a\right)\frac{1}{1-z_b}\ln \frac{2z_b}{1+z_b}+𝒪\left(1\right)$

Equations

$$p(P_a^{\mu })+p(P_b^{\mu })\to L(q^{\mu })+X(P_X^{\mu }),$$

$$x_a=\frac{Q}{E_{\mathrm{cm}}}{e}^{+Y},x_b=\frac{Q}{E_{\mathrm{cm}}}{e}^{-Y},\tau =x_a{x}_b,$$

$$\frac{\mathrm{d}\sigma }{\mathrm{d}{x}_a\mathrm{d}{x}_b}$$

$$\displaystyle\quad\times f^{\mathrm{thr}}_{i}\Bigl{[}x_{a}\Bigl{(}1+\frac{k^{-}}{Qe^{+Y}}\Bigr{)}\Bigr{]}f^{\mathrm{thr}}_{j}\Bigl{[}x_{b}\Bigl{(}1+\frac{k^{+}}{Qe^{-Y}}\Bigr{)}\Bigr{]}.$$

$${\hat{\sigma }}_{ij}(z_a,z_b)$$

$$x_{\mp }\equiv \frac{q^{\mp }}{P_{a,b}^{\mp }}=\frac{\sqrt{Q^2+q_T^2}}{E_{\mathrm{cm}}}{e}^{\pm Y}.$$

$${\lambda }_{\mathrm{QCD}}^2\ll {\lambda }^2\sim 1-x_{-}\ll 1\text{for generic x+},$$

$$p_X^{\mu }\sim (q^{+},P_a^{-}-q^{-},p_{X\perp })\sim (q^{+},{\lambda }^2{q}^{-},\lambda \sqrt{q^{+}{q}^{-}}).$$

$$\frac{\mathrm{d}\sigma }{\mathrm{d}{x}_{-}\mathrm{d}{x}_{+}}$$

$$\times B_j(t,x_{+}).$$

$$\displaystyle B_{j}(t,x)=\int\!\frac{\mathrm{d}z}{z}\,\mathcal{I}_{jk}(t,z)\,f_{k}\Bigl{(}\frac{x}{z}\Bigr{)}\,.$$

$$p_{\bar{n}}$$

$$P_n$$

$$\frac{\mathrm{d}\sigma }{\mathrm{d}{x}_{-}\mathrm{d}{x}_{+}\mathrm{d}{\vec{q}}_T}$$

$$\times B_j(t,{\vec{q}}_T,x_{+}).$$

$$\frac{\mathrm{d}\sigma }{\mathrm{d}{x}_a\mathrm{d}{x}_b\mathrm{d}{\vec{q}}_T}$$

$$\times B_j(t,{\vec{q}}_T,x_b).$$

$$\frac{\mathrm{d}\sigma}{\mathrm{d}x_{a}\mathrm{d}x_{b}}=H_{ij}(Q^{2})\!\int\!\mathrm{d}\tilde{t}\,f_{i}^{\mathrm{thr}}\Bigl{[}x_{a}\Bigl{(}1+\frac{\tilde{t}}{Q^{2}}\Bigr{)}\Bigr{]}\tilde{B}_{j}(\tilde{t},x_{b})\,.$$

$$\displaystyle\tilde{B}_{j}(\tilde{t},x)=\int\mathrm{d}^{2}\vec{\kern 0.0ptk}_{T}\,B_{j}\Big{(}\tilde{t}-\frac{k_{T}^{2}}{2},\vec{\kern 0.0ptk}_{T},x\Big{)}\,.$$

$$\tilde{B}_{j}(\omega k^{-}\!,x_{b})=\int\!\frac{\mathrm{d}k^{+}}{\omega}\,S(k^{-},k^{+})f_{j}^{\mathrm{thr}}\Bigl{[}x_{b}\Bigl{(}1+\frac{k^{+}}{\omega}\Bigr{)}\Bigr{]}$$

$$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}x_{a}\mathrm{d}x_{b}}=H_{ij}\Bigl{[}f_{i}^{\mathrm{thr}}\otimes\tilde{B}_{j}+\tilde{B}_{i}\otimes f_{j}^{\mathrm{thr}}-S\otimes f_{i}^{\mathrm{thr}}f_{j}^{\mathrm{thr}}\Bigr{]}\,,$$

$${\hat{\sigma }}_{ij}(z_a,z_b)$$

$$\displaystyle\quad\quad\quad-\delta_{ki}\,\delta_{\ell j}\,\hat{S}(z_{a},z_{b})\bigr{]}\,,$$

$${\hat{\mathcal{I}}}_{ij}(z_a,z_b)$$

$$\hat{S}(z_a,z_b)$$

$${\hat{\mathcal{I}}}_{ij}(z_a,z_b)$$

$${\hat{\sigma }}_{ij}(z_a,z_b)=\sum _{m_a,m_b}{\hat{\sigma }}_{ij}^{(m_a,m_b)}(z_a,z_b),$$

$${\hat{\sigma }}_{ij}^{(3)}$$

$$\displaystyle\quad+\mathcal{L}_{4}(1-z_{a})H^{(0)}_{ik}\biggl{[}\delta_{kj}\delta(1-z_{b})\Bigl{(}-\frac{2\beta_{0}}{3}-\frac{\gamma_{B\,0}^{i}}{2}\Bigr{)}$$

$$\displaystyle\quad+P^{(0)}_{kj}(z_{b})\biggr{]}\frac{5}{8}(\Gamma_{0}^{i})^{2}+\mathcal{L}_{3}(1-z_{a})\biggl{\{}H^{(1)}_{ij}\delta(1-z_{b})\frac{\Gamma_{0}^{i}}{2}$$

$$\displaystyle\quad+H^{(0)}_{ik}\biggl{[}\delta_{kj}\delta(1-z_{b})\Bigl{(}\Gamma_{1}^{i}-\frac{\pi^{2}}{6}(\Gamma_{0}^{i})^{2}+\frac{\beta_{0}^{2}}{3}$$

$$\displaystyle\quad+\frac{(\gamma_{B\,0}^{i})^{2}}{4}+\frac{5}{6}\beta_{0}\gamma_{B\,0}^{i}\Bigr{)}-P^{(0)}_{kj}(z_{b})\Bigl{(}\frac{5\beta_{0}}{3}+\gamma_{B\,0}^{i}\Bigr{)}$$

$$\displaystyle\quad+\bigl{(}P^{(0)}_{k\ell}\otimes P^{(0)}_{\ell j}\bigr{)}(z_{b})+\tilde{I}^{(1)}_{kj}(z_{b})\,\frac{\Gamma_{0}^{i}}{2}\biggl{]}\biggr{\}}\Gamma_{0}^{i}+\dotsb\,.$$

$$\frac{\mathrm{d}\sigma }{\mathrm{d}{x}_{-}\mathrm{d}{x}_{+}\mathrm{d}{\vec{q}}_T}$$

$${\lambda }^2\sim 1-x_{-}\ll 1\text{for generic }{x}_{+},{\lambda }_{\mathrm{QCD}}\equiv {\Lambda }_{\mathrm{QCD}}/Q\sim {\lambda }^2\ll \lambda .$$

Full text

††preprint:

DESY 19-135

NIKHEF 2019-021

I Introduction

Color-singlet processes play a central role in the LHC physics program. The $pp\to Z,W$ Drell-Yan processes are precision benchmarks, providing determinations of electroweak parameters and important inputs for fits of parton distribution functions (PDFs). Higgs and diboson processes provide strong sensitivity to possible contributions beyond the Standard Model.

We consider the production of a generic color-singlet final state $L$ together with hadronic radiation $X$,

[TABLE]

at hadronic center-of-mass energy $E_{\mathrm{cm}}^{2}=(P_{a}+P_{b})^{2}$. The key observables characterizing $L$ are its total invariant mass $Q\equiv\sqrt{q^{2}}$, rapidity $Y$, and transverse momentum $q_{T}\equiv\lvert\vec{\kern 0.0ptq}_{T}\rvert$. We define the momentum fractions

[TABLE]

which are equivalent to $Q$ and $Y$. The cross section differential in $x_{a,b}$ is given by Bodwin (1985); Collins et al. (1985, 1988)

[TABLE]

where $\hat{\sigma}_{ij}(z_{a},z_{b})$ denotes the perturbatively calculable partonic cross section and $f_{i,j}(x)$ are the standard PDFs. We always implicitly sum over parton indices $i,j$, and keep the dependence on renormalization scales $\mu$ implicit.

In the soft threshold limit $\tau\to 1$, which implies that both $x_{a,b}\to 1$, Eq. (3) factorizes further Sterman (1987); Catani and Trentadue (1989); Ravindran et al. (2007); Westmark and Owens (2017); Banerjee et al. (2018a, b),

[TABLE]

In this limit, the hadronic final state is forced to be soft, and is described by the soft function $S$, which encodes soft-gluon emissions from the colliding hard partons. Furthermore, only the hard Born processes, e.g. $q\bar{q}\to Z$ or $gg\to H$, contribute. They are encoded in the hard function $H_{ij}$, including the Born-like virtual corrections. Any nondiagonal partonic channels like $qg\to Lq$ vanish for $\tau\to 1$. The threshold PDF $f^{\mathrm{thr}}_{i}(x)$ encodes the extraction of parton $i$ from the proton for $x\to 1$. At partonic level, Eq. (4) implies that for $z\equiv z_{a}z_{b}\to 1$, which requires both $z_{a,b}\to 1$,

[TABLE]

up to power corrections in $1-z$.

While taking $\tau\to 1$ forces $z\to 1$, even for typical LHC values of $\tau\ll 1$, the $z\sim 1$ region often numerically dominates the cross section. Soft threshold factorization has thus been widely used for decades. It enables the all-order resummation of the leading terms in $1-z$, see e.g. Refs. Sterman (1987); Catani and Trentadue (1989); Ravindran et al. (2007); Westmark and Owens (2017); Banerjee et al. (2018a, b); Appell et al. (1988); Magnea (1991); Korchemsky and Marchesini (1993); Contopanagos et al. (1997); Catani et al. (1996); Belitsky (1998); Moch and Vogt (2005); Laenen and Magnea (2006); Idilbi et al. (2006); Mukherjee and Vogelsang (2006); Bolzoni (2006); Becher et al. (2008); Bonvini et al. (2011, 2015); Fuks et al. (2013); Bonvini and Marzani (2014); Schmidt and Spira (2016); A H et al. (2019). The resummation at next-to-leading power (NLP) in $1-z$ has also received recent interest Bonocore et al. (2016); Del Duca et al. (2017); Beneke et al. (2019). Another important use is to approximate the fixed-order cross section by expanding in $1-z$, e.g. at N3LO Anastasiou et al. (2014); Ahmed et al. (2014a, b); de Florian et al. (2014); Li et al. (2015); Anastasiou et al. (2016); Dulat et al. (2018, 2019).

In this letter, we derive the factorization that generalizes Eqs. (4) and (5) to the weaker limit where only one of $x_{a,b}$ (or $z_{a,b}$) approaches $1$ while keeping the exact dependence on the other variable. This corresponds to the kinematic limit $\lvert Y\rvert\to Y_{\mathrm{max}}=\ln(E_{\mathrm{cm}}/Q)$ for generic (including small) $Q$ values.

II Generalized threshold factorization

We use light-cone coordinates $p^{\mu}\equiv(n\cdot p,\bar{n}\cdot p,p_{\perp})\equiv(p^{+},p^{-},p_{\perp})$ with respect to lightlike vectors $n^{\mu}\equiv(1,\hat{z})$ and $\bar{n}^{\mu}\equiv(1,-\hat{z})$ along the beam axis $\hat{z}$. We first consider the observables $q^{\mp}$ instead of $Q$ and $Y$, with corresponding momentum fractions

[TABLE]

We consider the generalized threshold limit

[TABLE]

where $\lambda_{\mathrm{QCD}}\equiv\Lambda_{\mathrm{QCD}}/Q$ and $\lambda$ are power-counting parameters. In this limit, illustrated in Fig. 1, $L$ has large $Y$ while the emissions in $X$ become collimated in the opposite direction with typical momenta

[TABLE]

In this situation, the following factorization theorem holds to leading power in $1-x_{-}$,

[TABLE]

Here, $B_{j}(t,x)$ is the inclusive beam function that also appears in the factorization for $N$-jettiness Stewart et al. (2010a, b). It depends on the transverse virtuality $t$ and momentum fraction $x$ of the colliding parton $j$. Since $t\sim p_{X}^{2}\sim\lambda^{2}Q^{2}\gg\Lambda_{\mathrm{QCD}}^{2}$, it can be calculated perturbatively in terms of standard PDFs as Stewart et al. (2010a, c)

[TABLE]

The matching coefficients $\mathcal{I}_{jk}$ are known to $\mathcal{O}(\alpha_{s}^{2})$ Stewart et al. (2010c); Berger et al. (2011); Gaunt et al. (2014a, b), with progress at $\mathcal{O}(\alpha_{s}^{3})$ Melnikov et al. (2019a, b).

To derive Eq. (II), we use Soft-Collinear Effective Theory (SCET) Bauer et al. (2000, 2001); Bauer and Stewart (2001); Bauer et al. (2002a, b), which is the effective theory of QCD in the limit $\lambda\ll 1$. The derivation proceeds in a standard fashion Stewart et al. (2010a), with more details given in *[SeeSupplementalMaterialattheendofthispreprint][]supplement. The key elements are the necessary degrees of freedom (modes) in the effective theory,

[TABLE]

The $p_{\bar{n}}$ collinear modes describe the QCD final state at the scale $\lambda Q$, which due to Eq. (8) is ${\bar{n}}$-collinear. The $P_{n,{\bar{n}}}$ collinear modes describe the PDFs at the scale $\lambda_{\mathrm{QCD}}Q$. The $P_{s}$ modes describe the soft interactions between the $p_{\bar{n}}$ and $P_{n}$ modes. Possible ultrasoft/Glauber modes $P_{us,G}\sim Q(\lambda_{\mathrm{QCD}}^{2},\lambda_{\mathrm{QCD}}^{2},\lambda_{\mathrm{QCD}}^{(2)})$ cancel as in Eq. (3) because the measurement is the same and its scale $\lambda Q\gg\lambda_{\mathrm{QCD}}Q$.

The hard function $H_{ij}$ in Eq. (II) is the same as in Eq. (4). It arises from matching the electroweak current for $L$ onto the corresponding SCET current. There are no interactions between the modes in the leading-power SCET Lagrangian, so the complete matrix element factorizes into separate ones in each sector. The matrix element of the combined $p_{\bar{n}}$ and $P_{\bar{n}}$ modes yields $B_{j}(t,x)$, and their separation leads to Eq. (10) Stewart et al. (2010a, c). The combined matrix element of the $P_{n}$ and $P_{s}$ modes yields $f_{i}^{\mathrm{thr}}$ Fleming and Labun (2015); Hoang et al. (2016). The convolution structure in Eq. (II) follows from momentum conservation sup .

Next, we consider also measuring $\vec{\kern 0.0ptq}_{T}$. From Eq. (8), it follows that generically $q_{T}\sim p_{X\perp}\sim\lambda Q$, so the $q_{T}$ dependence is entirely described by the $p_{\bar{n}}$ modes, which yields the factorization theorem

[TABLE]

Here, $B_{j}(t,\vec{\kern 0.0ptk}_{T},x)$ is the double-differential beam function Jain et al. (2012); Gaunt and Stahlhofen (2014) that also occurs in the joint resummation of $q_{T}$ and [math]-jettiness Procura et al. (2015); Lustermans et al. (2019).

We can now change variables in Eq. (II) to $x_{a,b}$, using Eq. (6) and expanding in $\lambda$, which yields

[TABLE]

Crucially, when expanding $x_{-}=x_{a}[1+q_{T}^{2}/(2Q^{2})+\mathcal{O}(\lambda^{4})]$, we have to keep the $q_{T}^{2}/(2Q^{2})\sim\lambda^{2}$ term in the PDF argument because it is of the same order as $t/Q^{2}\sim\lambda^{2}$. Integrating Eq. (II) over $\vec{\kern 0.0ptq}_{T}$, we obtain

[TABLE]

The factorization theorems in Eqs. (II) and (14) hold at leading power in the generalized threshold limit $\lambda^{2}\sim 1-x_{a}\ll 1$ for generic $x_{b}$. They are our key new results.

In Eq. (14) we changed variables to $\tilde{t}=t+q_{T}^{2}/2$, and defined the new modified beam function

[TABLE]

It has the same $\mu$ evolution as $B_{j}(t,x)$ but different constant terms. It obeys a matching relation analogous to Eq. (10). Using the known results for $B_{j}(t,\vec{\kern 0.0ptk}_{T},x)$ Jain et al. (2012); Gaunt and Stahlhofen (2014), we have calculated its matching coefficients $\tilde{\mathcal{I}}_{jk}(\tilde{t},z)$ to $\mathcal{O}(\alpha_{s}^{2})$ for $j=q$ and $\mathcal{O}(\alpha_{s})$ for $j=g$ sup .

The factorization structure in Eqs. (II) and (14) turns out to be analogous to deep-inelastic scattering (DIS) at large Bjorken $x$ Sterman (1987); Catani and Trentadue (1989); Manohar (2003); Becher et al. (2007); Fleming and Labun (2015); Chay and Kim (2013); Hoang et al. (2016), which factorizes as $H_{ij}f_{i}^{\rm thr}\otimes J_{j}$, with the jet function $J_{j}$ describing collimated final-state radiation. For Eqs. (II) and (14) to be consistent, the $\mu$ dependence of the functions must cancel between them, and it does so in the same way as for DIS. The beam function is known to have the same $\mu$ evolution as the jet function Stewart et al. (2010c). The existence of the above factorization theorems provides an independent confirmation of this. The key differences to DIS are the additional dependence on $x_{b}$ and the nontrivial $q_{T}$ dependence. The latter bears some resemblance to different $1$-jettiness definitions in exclusive DIS Kang et al. (2013).

Since Eq. (14) is valid for $x_{a}\to 1$ and arbitrary $x_{b}$, it must contain the soft threshold factorization in Eq. (4) for $x_{b}\to 1$ as a special case. This implies sup

[TABLE]

to leading power in $1-x_{b}$, and identically for $B_{j}$.

We can now combine Eq. (14) with the analogous result in the opposite limit $x_{b}\to 1$ for generic $x_{a}$ by adding the two and subtracting their overlap, which is precisely given by the soft limit. This yields our main result, the generalized threshold factorization theorem

[TABLE]

with the convolutions as in Eqs. (4) and (14). Analogous results hold for $x_{\pm}$ and differential in $\vec{\kern 0.0ptq}_{T}$. Replacing $f_{i}^{\mathrm{thr}}[x(1+1-z)]$ by $f_{i}(x/z)/z$, which is justified at leading power in $1-z$, and comparing to Eq. (3), we obtain the corresponding partonic factorization theorem

[TABLE]

where we changed variables to $z_{a,b}$ and defined

[TABLE]

As illustrated in Fig. 2, $\hat{\mathcal{I}}_{ij}(z_{a},z_{b})$ captures all terms in $\hat{\sigma}(z_{a},z_{b})$ that are singular as $z_{a}\to 1$ with their exact dependence on $z_{b}$, and vice versa, with the overlap given by $\hat{S}(z_{a},z_{b})$. This includes flavor-nondiagonal contributions, e.g., $ij=qg$ with $k\ell=q{\bar{q}}$ or $gg$ in Eq. (II).

III Validation

A nontrivial validation is to check Eq. (17) against the full fixed-order result. First, Eq. (16) implies

[TABLE]

which we verified analytically to $\mathcal{O}(\alpha_{s}^{2})$. It then suffices to check that Eq. (14) reproduces the fixed-order result for $x_{a}\to 1$ (or $z_{a}\to 1$) for generic $x_{b}$ (or $z_{b}$).

At NLO, we can check analytically for Drell-Yan and $gg\to H$ against the results of Refs. Anastasiou et al. (2003a, b), which are given as distributions in $z=z_{a}z_{b}$ and a variable $y(z_{a},z_{b})$. Their translation to $(z_{a},z_{b})$ is nontrivial sup . We find complete agreement for all partonic channels. We also combined all singular terms from Eq. (II) with the regular terms from Refs. Anastasiou et al. (2003a, b) to construct the full $\hat{\sigma}_{ij}(z_{a},z_{b})$ at NLO. The result is given in sup and agrees with Refs. Kubar et al. (1980); Mathews et al. (2005); Ravindran et al. (2007). 111 As discussed further in sup , several soft threshold factorizations differential in rapidity Bolzoni (2006); Mukherjee and Vogelsang (2006); Becher et al. (2008); Bonvini et al. (2011, 2015) differ from Eq. (4) and do not reproduce the correct soft limit already at NLO. .

At NNLO, we numerically validate our own implementation of Eq. (14) in SCETlib Ebert et al. (2018) against Vrap Anastasiou et al. (2004). We use flat PDFs, $f_{i}^{\mathrm{thr}}(x)=f_{i}(x)=\theta(1-x)$, which amounts to taking cumulant integrals of the partonic cross section and provides the strongest possible numerical check. In Fig. 3, we compare the $\mathcal{O}(\alpha_{s}^{2})$ contribution for Drell-Yan as a function of $1-x_{a}$ at fixed $x_{b}=10^{-2}$. We find perfect agreement. The breakdown into partonic channels is given in sup . We also find similar agreement for other $x_{b}$ and for $pp\to Z/\gamma^{\ast}$ on the resonance.

IV Illustrative applications

The immediate question that arises is how well the generalized threshold limit approximates the full fixed-order result for physical PDFs, particularly in comparison to the soft limit. We use the MMHT2014nnlo68cl Harland-Lang et al. (2015) PDFs, Vrap Anastasiou et al. (2004) to obtain the full NNLO result222 The public Vrap 0.9 assumes $f_{q}(x)=f_{\bar{q}}(x)$ for $q=s,c,b$. We modified it to allow for different sea quark and antiquark PDFs.

and SCETlib Ebert et al. (2018) to implement Eq. (II). In Fig. 4 we compare the $\mathcal{O}(\alpha_{s})$ and $\mathcal{O}(\alpha_{s}^{2})$ contributions to the Drell-Yan rapidity spectrum at $Q=m_{Z}$, separated into quark channels ($q\bar{q}+qq^{\prime}$) and channels involving a gluon ($qg+gq+gg$). Analogous results at $\mu=Q/2$ and for $gg\to H$ are provided in sup . The generalized threshold limit approximates the full result well for all channels and all $Y$. As expected, it works particularly well toward large $Y$. It works significantly better than the soft limit, which only provides a poor approximation for the $q\bar{q}$ channel and none for the others. Currently, all ingredients are available to perform the resummation in the generalized threshold limit to N3LL, which we leave to future work. We expect it to be much more powerful for improving the precision of perturbative predictions than the soft threshold resummation.

Next, we may ask how well the power expansion around the generalized threshold limit works and how it relates to the soft expansion. Consider the double expansion in $1-z_{a}$ and $1-z_{b}$, illustrated in Fig. 2,

[TABLE]

where $\hat{\sigma}_{ij}^{(m_{a},m_{b})}(z_{a},z_{b})\sim(1-z_{a})^{m_{a}}(1-z_{b})^{m_{b}}$. Expanding around the soft $z=z_{a}z_{b}\to 1$ limit corresponds to counting powers of $(1-z)^{m_{a}+m_{b}+1}$. The leading-power result in Eq. (5) gives the $m_{a}=m_{b}=-1$ term. At the $m$th order, NmLPsoft, we keep all terms with $m_{a}+m_{b}+2\leq m$. At leading power in the generalized expansion, Eq. (II) includes all terms with $\min\{m_{a},m_{b}\}=-1$. Similarly, at the $m$th order, NmLPgen, we keep all terms with $\min\{m_{a},m_{b}\}=m-1$, so the missing corrections at NmLPgen are $\mathcal{O}[(1-z_{a})^{m}(1-z_{b})^{m}]$.

In Fig. 5, we show the deviation from the exact result at various orders in both expansions for Drell-Yan at NLO, where we have full analytic control. Analogous results for $gg\to H$ are provided in sup . The generalized expansion performs significantly better than the soft one for both partonic channels. We expect this to hold in general, since expanding a two-dimensional function along a one-dimensional boundary is superior to expanding it in a single point on that boundary.

In fact, as seen in Fig. 2, each order in the generalized expansion fully contains two orders in the soft expansion, and in particular, the LPgen result Eq. (II) contains the entire NLPsoft contribution. This does not mean it can be used to perform the NLPsoft resummation because the $\mu$ evolution of $B_{i}(t,x)$ does not predict its $x$ dependence. It does, however, show that $H_{ij}$ factorizes for all partonic channels and reduces the problem to deriving the NLPsoft factorization for Eq. (16).

At N3LO, Eq. (II) predicts a highly nontrivial set of terms for any color-singlet process, since all terms $\sim\mathcal{L}_{n}(1-z_{a})$ in $\hat{\mathcal{I}}_{ij}(z_{a},z_{b})$ are known from its $\mu$ evolution, where $\mathcal{L}_{n}(y)\equiv[\ln^{n}(y)/y]_{+}$. To illustrate this, the coefficient of $\alpha_{s}^{3}/(4\pi)^{3}$ in $\hat{\sigma}_{ij}(\mu=Q)$ with $n\geq 3$ is

[TABLE]

The required ingredients are given in sup . The extension down to $\mathcal{L}_{0}(1-z_{a})$ and to the full $z_{a}$ dependence for $z_{b}\to 1$ is straightforward Billis et al. . The $\delta(1-z_{a})$ coefficient requires the still unknown $\mathcal{O}(\alpha_{s}^{3})$ finite terms of the beam function. For Drell-Yan, Eq. (22) significantly extends the current knowledge at $\mathcal{O}(\alpha_{s}^{3})$ Ahmed et al. (2014b), providing the full $z_{b}$ dependence for all partonic channels. For $gg\to H$, the extension to $\mathcal{L}_{0}(1-z_{a})$ would also provide additional information beyond what is currently known Dulat et al. (2019).

V Summary

We introduced a generalized threshold factorization, which is much more powerful than the often used soft threshold, thus opening the door to numerous applications to improve theoretical predictions for collider processes. It describes all kinematic limits in $(x_{a},x_{b})$ or $(Q,Y)$, including in particular $\lvert Y\rvert\to Y_{\mathrm{max}}$ at generic $Q$, which is directly accessible at the LHC. It enables the corresponding threshold resummation where only one PDF is probed at large $x$, which is not captured by the soft limit. At the partonic level, it captures all singularities of $\hat{\sigma}_{ij}(z_{a},z_{b})$, including nondiagonal partonic channels. It can be used to predict a rich set of terms at higher fixed order or to resum them to all orders. It is the weakest known limit in which the process-dependent virtual corrections ($H_{ij}$) factorize.

While we only considered color-singlet processes here, the same methods can be used to generalize the soft threshold factorization in other situations, such as for processes with heavy particles, jets, or identified hadrons in the final state.

Acknowledgements.

Acknowledgments

We thank M. Beneke, G. Billis, A. Broggio, G. Das, M. Diehl, L. Dixon, M. Ebert, S. Forte, A. Kulesza, S. Marzani, B. Mistlberger, F. Ringer, D. Scott, D. Soper, and W. Waalewijn for discussions. We thank M. Stahlhofen for providing the results of Ref. Gaunt and Stahlhofen (2014), and M. Prakash and V. Ravindran for help in comparing with Refs. Mathews et al. (2005); Ravindran et al. (2007). The authors thank each other’s institutions for hospitality. This work was supported in part by the ERC grant ERC-STG-2015-677323 and the D-ITP consortium, a program of the Netherlands Organization for Scientific Research (NWO) that is funded by the Dutch Ministry of Education, Culture and Science (OCW).

Supplemental material

V.1 Factorization theorem

Here we give some more details on the derivation of the factorization theorem in Eq. (II), which underlies all other factorization theorems, and which we repeat here for easy reference,

[TABLE]

As in the main text, we always implicitly sum over repeated flavor indices $i,j,k$. Equation (S1) is valid up to power corrections in $\lambda^{2}$ in the generalized threshold limit

[TABLE]

We require $\lambda_{\mathrm{QCD}}\ll\lambda$ for reasons that will be apparent soon. Without loss of generality we can then consider $\lambda_{\mathrm{QCD}}\sim\lambda^{2}$. This relation is to be interpreted as follows: First, in our context, $\lambda_{\mathrm{QCD}}$ denotes the scale of the PDFs, which is generically allowed to be as large as $\lambda^{2}$ and does not necessarily have to be nonperturbative. If it happens to be a perturbative scale, then the physics below $\lambda_{\mathrm{QCD}}$ is simply described by the PDF evolution. Conversely, it also means that $\lambda^{2}$ is in principle allowed to be as small as $\lambda_{\mathrm{QCD}}$ including being nonperturbative, i.e., it is only relevant that $\lambda\gg\lambda_{\mathrm{QCD}}$ is perturbative.

The key step in deriving Eq. (S1) is to identify the relevant degrees of freedom (modes) in the effective field theory (EFT) that describe the physical situation. They are defined by the relative scaling of their light-cone momentum components and are summarized in Table S1. We note that rather than matching QCD directly onto these modes, one may also perform a multi-stage matching, as was done for endpoint DIS in Ref. Hoang et al. (2016), which yields the same end result.

The $p_{{\bar{n}}}$ modes describe the hadronic final state of the collision. Their scaling is determined by the fact that in the limit of Eq. (S2), there is only $p_{\bar{n}}^{-}\sim\lambda^{2}E_{\mathrm{cm}}\sim\lambda^{2}q^{-}$ minus momentum available. On the other hand, their plus momentum is unconstrained, which means it has generic scaling set by the hard interaction, $p_{\bar{n}}^{+}\sim\xi_{b}E_{\mathrm{cm}}\sim x_{b}E_{\mathrm{cm}}\sim q^{+}$. Since $p_{\bar{n}}^{2}\sim\lambda^{2}q^{+}q^{-}\sim\lambda^{2}Q^{2}\gg\Lambda_{\mathrm{QCD}}^{2}$, the $p_{\bar{n}}$ modes are perturbative. Therefore, they describe the perturbative QCD final state produced in the partonic collision in addition to $L$ (but excluding the beam remnant). The $P_{n}$ and $P_{\bar{n}}$ modes describe the incoming protons, or more precisely, the partons in the proton with the typical momentum fractions required to produce the hard final state. This means their scaling is determined by $P_{n}^{-}\sim q^{-}$ and $P_{\bar{n}}^{+}\sim\xi_{b}E_{\mathrm{cm}}\sim x_{b}E_{\mathrm{cm}}\sim q^{+}$ and $P_{n}^{2}\sim P_{\bar{n}}^{2}\sim\Lambda_{\mathrm{QCD}}^{2}$.

The soft $P_{s}$ modes describe the interactions between the $p_{\bar{n}}$ and $P_{n}$ modes. Their scaling is thus determined by $P_{s}^{-}\sim p_{\bar{n}}^{-}\sim\lambda^{2}Q$ and $P_{s\perp}\sim P_{n\perp}\sim\Lambda_{\mathrm{QCD}}$ or equivalently $P_{s}^{2}\sim P_{n}^{2}\sim\Lambda_{\mathrm{QCD}}^{2}$. Hence, they keep the $p_{\bar{n}}$ modes on shell and have a SCETI-like relation to them. Their interactions with the $p_{\bar{n}}$ modes in the leading-power SCET Lagrangian are decoupled and moved into soft Wilson lines in the SCET current via the BPS field redefinition. At the same time, the $P_{s}$ modes have a SCETII-like relation to the $P_{n}$ modes, i.e., they have the same virtuality but are parametrically separated in rapidity. Hence, their interactions with the $P_{n}$ modes, which take the $P_{n}$ modes off shell, are described by soft Wilson lines in the SCET current that are directly produced during the matching onto SCET. The distinction of the $P_{s}$ modes relies on $\lambda_{\mathrm{QCD}}\ll\lambda$, while for $\lambda_{\mathrm{QCD}}\sim\lambda$, they would become degenerate with the $p_{\bar{n}}$ and $P_{\bar{n}}$ modes.

The power counting and the relations between the modes are easiest in the leptonic frame, which is the frame where the color-singlet final state has total rapidity $\hat{Y}=0$. Boosting from the lab frame to the leptonic frame by $Y$, we have $\hat{q}^{\pm}=q^{\pm}e^{\pm Y}=\sqrt{q^{+}q^{-}}\sim Q$. In the leptonic frame, the $p_{\bar{n}}$ modes are genuinely ${\bar{n}}$-collinear with $p_{\bar{n}}^{-}\sim\lambda^{2}p_{\bar{n}}^{+}$, and the soft modes are homogeneous, $P_{s}\sim\lambda^{2}Q\sim\lambda_{\mathrm{QCD}}Q$. By contrast, in the lab frame we must separately keep track of $q^{+}$ and $q^{-}$, i.e., we cannot count them as $q^{+}\sim q^{-}$, because we want to take the limit of large $q^{-}$ for generic $q^{+}$. As a result, the $p_{\bar{n}}$ modes do not necessarily appear to be ${\bar{n}}$-collinear in the lab frame because $q^{+}\sim\lambda^{2}q^{-}$ is allowed. However, the key requirement for their factorization is that they are collinear relative to the soft modes, which in the lab frame are boosted in the $n$-collinear direction and become $n$-collinear-soft (csoft) modes Bauer et al. (2012); Procura et al. (2015).

Finally, the ultrasoft (usoft) $P_{us}$ and Glauber $P_{G}$ modes describe the interactions between the $P_{n}$ and $P_{\bar{n}}$ modes that are possible without pushing either of them off shell, which requires $P_{us}^{-}\sim P_{\bar{n}}^{-}$ and $P_{us}^{+}\sim P_{n}^{+}$. The $\perp$ component of the usoft modes is fixed by requiring them to be on-shell modes, $P_{us\perp}^{2}\sim P_{us}^{+}P_{us}^{-}$. The corresponding Glauber modes are allowed to be off shell, so their $\perp$ component can be as large as $P_{G\perp}\sim P_{n\perp}\sim P_{{\bar{n}}\perp}\sim\Lambda_{\mathrm{QCD}}$. The effects of the usoft and Glauber modes cancel, so we do not need to consider them further. This directly follows from the collinear factorization theorem Bodwin (1985); Collins et al. (1985, 1988) because in the limit we consider, the measurement is still fully inclusive over perpendicular momenta at the scale $\lambda_{\mathrm{QCD}}Q$. This is another reason why we require $\lambda_{\mathrm{QCD}}\ll\lambda$. Note also that there is only a single collinear sector at the perturbative $\lambda Q$ scale, so there are no perturbative Glauber modes with scaling $Q(\lambda^{2},\lambda^{2},\lambda)$ that could spoil factorization.

Since there are no interactions between the modes in the leading-power SCET Lagrangian, the cross section factorizes into separate forward matrix elements in each sector. The detailed derivation closely follows Ref. Stewart et al. (2010a), with the matrix element of the combined $p_{\bar{n}}$ and $P_{\bar{n}}$ modes giving $B_{j}(t,x)$ and the combined matrix element of the $P_{n}$ and $P_{s}$ modes giving $f_{i}^{\mathrm{thr}}$. The factorization of the threshold PDF into separate collinear and csoft nonperturbative matrix elements is discussed in Refs. Fleming and Labun (2015); Hoang et al. (2016). It is not needed for our purposes. The arguments of the functions and their convolution structure in Eq. (S1) follow from overall momentum conservation among all sectors, which at leading power in $\lambda$ must hold separately for the large (label) and small (residual) momenta carried by the matrix elements. Denoting them as $\omega_{n,{\bar{n}}}\sim Q$, $k_{{\bar{n}}\perp}\sim\lambda Q$, and $k_{s,{\bar{n}}}\sim\lambda^{2}Q$, it takes the form

[TABLE]

The first three $\delta$ functions set the $x_{\pm}$ and $\vec{\kern 0.0ptq}_{T}$ arguments of $f_{i}^{\rm thr}$ and $B_{j}$. The analogous $\delta(k_{s}^{+}-k_{\bar{n}}^{+})$ disappeared by absorbing the $k_{\bar{n}}^{+}$ dependence into $\omega_{\bar{n}}$ Stewart et al. (2010a). By combining the $P_{n}$ and $P_{s}$ modes, the threshold PDF depends on $(\omega_{n}+k_{s}^{-})/P_{a}^{-}=x_{-}(1+k_{s}^{-}/q^{-})$. Here, we do have to keep track of the momentum of the $P_{s}$ modes $k_{s}^{-}\sim\lambda^{2}Q$, as it is much larger than the typical residual minus momentum $k_{n}^{-}\sim\lambda_{\mathrm{QCD}}^{2}Q$ of the $P_{n}$ modes that is absorbed into $\omega_{n}$. The $\delta(k_{s}^{-}-k_{\bar{n}}^{-})$ then yields the convolution in $t=q^{+}k_{\bar{n}}^{-}$ in Eq. (S1). While Eq. (S1) is formally derived in the leptonic frame, it has exactly the same form in the lab frame. This is because the measured observables $x_{\pm}$ and $\vec{\kern 0.0ptq}_{T}$ are boost invariant along the beam axis, and reparametrization invariance (RPI) Manohar et al. (2002) forces all functions to only depend on boost-invariant quantities.

To see that the generalized threshold limit contains the soft threshold limit, note that the limits $x_{-}\to 1$ and $x_{+}\to 1$ commute, so taking one limit after the other is equivalent to taking $x_{-},x_{+}\to 1$ simultaneously. To see this, consider the hierarchy $\lambda_{-}^{2}\sim 1-x_{-}\ll\lambda_{+}^{2}\sim 1-x_{+}\ll 1$, which we can interpret as taking $x_{+}\to 1$ after having already taken $x_{-}\to 1$. In this limit, the $p_{\bar{n}}$ modes factorize into perturbative ${\bar{n}}$-csoft modes $p_{{\bar{n}},cs}\sim Q(\lambda_{+}^{2},\lambda_{-}^{2},\lambda_{+}\lambda_{-})$ and ${\bar{n}}$-csoft modes $P_{{\bar{n}},cs}\sim Q(\lambda_{+}^{2},\lambda_{\mathrm{QCD}}^{2}/\lambda_{+}^{2},\lambda_{\mathrm{QCD}})$. Including $\lambda_{+}$, the condition $\lambda_{\mathrm{QCD}}\ll\lambda$ becomes $\lambda_{\mathrm{QCD}}\ll\lambda_{-}\lambda_{+}$ and without loss of generality we can consider $\lambda_{\mathrm{QCD}}\sim\lambda_{-}^{2}\lambda_{+}^{2}$. This also implies that the $P_{s}$ modes now get boosted in the $n$ direction and become $n$-csoft modes, $P_{s}\equiv P_{n,cs}$. The beam function matching onto PDFs now takes the form

[TABLE]

where the combined $P_{{\bar{n}},cs}$ and $P_{\bar{n}}$ modes yield the threshold PDF, and $\mathcal{S}(k^{-},k^{+},\vec{\kern 0.0ptk}_{T})$ is the matrix element of the $p_{{\bar{n}},cs}$ modes. It has the same Wilson line structure as the soft function appearing in Eq. (4), and thus upon integration over $\vec{\kern 0.0ptk}_{T}$ becomes equal to it to all orders by reparametrization invariance. We now have $k_{T}^{2}/Q^{2}\sim\lambda_{-}^{2}\lambda_{+}^{2}\ll t/Q^{2}\sim\lambda_{-}^{2}$, so $B_{j}$ and $\tilde{B}_{j}$ become the same and integrating Eq. (S4) over $\vec{\kern 0.0ptk}_{T}$ yields Eq. (16) for either of them.

V.2 Plus distribution identities

In Refs. Anastasiou et al. (2003a, b, 2004), the partonic cross section for $(Q,Y)$ is given in terms of partonic variables $(z,y)$, defined as

[TABLE]

where $z_{a,b}$ are defined by Eq. (3), and the integration limits $0\leq z_{a,b}\leq 1$ correspond to $0\leq z\leq 1$ and $0\leq y\leq 1$. In the following, we derive relations between plus distributions in $(z,y)$ and $(z_{a},z_{b})$.

In general, the plus distributions are uniquely defined by their functional form in the bulk, i.e. for $z_{a,b}<1$ away from any singularity, and their integrals (against unit test functions) over arbitrary integration regions that include the singularities. First, the functional form in the bulk is easily obtained simply by plugging in Eq. (V.2). Next, the correct boundary terms at $z_{a}=1$ as a function of $z_{b}$, at $z_{b}=1$ as a function of $z_{a}$, and at $z_{a}=z_{b}=1$ are determined by comparing integrals over the integration region $x_{a}\leq z_{a}(z,y)\leq 1$ and $x_{b}\leq z_{b}(z,y)\leq 1$ for generic $x_{a}$ and $x_{b}$. This integration region is indicated by the gray box in the $(z_{a},z_{b})$ plane in Fig. S1. It is sufficiently general to fix all boundary terms and precisely corresponds to the relevant integration region for the physical cross section in $(Q,Y)$. In terms of $(z,y)$, the integration region is given by

[TABLE]

where the integration bounds in $y$, also illustrated in Fig. S1, are given by

[TABLE]

In Table S2, we collect the distribution identities at leading power in $1-z_{a}$ and arbitrary $z_{b}$ that are required for validating Eq. (14) at NLO, where we denote the plus distributions as

[TABLE]

The relations are derived by integrating each structure in the left column in terms of $(z,y)$ over the region in Eq. (V.2), expanding the result to leading power in $1-x_{a}$, and comparing to the corresponding (straightforward) integral over the same region in terms of $z_{a,b}$ of each structure in the right column.

The last entry in Table S2 has the most intricate structure. Its exact integral without any expansion is given by

[TABLE]

where the imaginary parts from the branch cuts cancel between the different terms. Matching this with the exact distribution in the bulk, we obtain the distributional identity

[TABLE]

Expanding the right-hand side to leading power in $1-z_{a}$ yields the result given in the last line of Table S2. We note that Eq. (V.2) is the two-dimensional analog of a typical distribution identity like

[TABLE]

Namely, it expresses plus distributions of a function $y(z_{a},z_{b})$ in terms of simpler plus distributions of $1-z_{a,b}$ plus regular terms. Moving the regular terms out of the plus distribution incurs additional boundary terms.

A key property of the left-hand side of Eq. (V.2) is that it vanishes when integrated over all of $y$. It is instructive to see how this is reproduced by the right-hand side by projecting onto $z$ via $\int\!\mathrm{d}z_{a}\mathrm{d}z_{b}\,\delta(z-z_{a}z_{b})$. The only nontrivial projection integral that is required is

[TABLE]

These terms precisely cancel the contributions from the other terms upon projecting onto $z$.

V.3 Exact NLO partonic cross sections in terms of $(z_{a},z_{b})$

Since Eq. (II) captures the full singularity structure as $z_{a}\to 1$ and/or $z_{b}\to 1$, the power corrections to it are of relative $\mathcal{O}[(1-z_{a})(1-z_{b})]$ and are fully integrable. Hence, we can construct the exact partonic cross section in terms of $(z_{a},z_{b})$ from the results in terms of $(z,y)$ from Refs. Anastasiou et al. (2003b, a) as

[TABLE]

Here, the term in square brackets is evaluated in the bulk, away from any singularities, so we can simply plug in Eq. (V.2). In the following, we collect the resulting expressions for the Drell-Yan and $gg\to H$ cross sections to NLO written explicitly in terms of $(z_{a},z_{b})$.

V.3.1 Results for Drell-Yan

The Born cross section for Drell-Yan production, $q\bar{q}\to Z/\gamma^{\ast}\to\ell^{+}\ell^{-}$, is given by

[TABLE]

with $N_{c}$ the number of colors, $Q_{q}$ the quark charge in units of $\lvert e\rvert$, $v_{\ell,q}$ and $a_{\ell,q}$ the standard electroweak vector and axial couplings of the leptons and quarks, and $m_{Z}$ and $\Gamma_{Z}$ the mass and width of the $Z$ boson. To restrict to $q\bar{q}\to\gamma^{\ast}\to\ell^{+}\ell^{-}$, only the first term $\propto Q_{q}^{2}$ is kept. The complete LO cross section is given by

[TABLE]

where the sum runs over $q=\{u,d,c,s,b\}$. We expand the partonic cross section for Drell-Yan as

[TABLE]

The LO result corresponding to Eq. (S15) is given by $\hat{\sigma}^{(0)}_{q{\bar{q}}}(z_{a},z_{b},Q,\mu)=\sigma_{B,q}^{\mathrm{DY}}\,\delta(1-z_{a})\,\delta(1-z_{b})$. Writing ${\bar{z}}_{a}\equiv 1-z_{a}$ and ${\bar{z}}_{b}\equiv 1-z_{b}$ for short, the NLO coefficient for the $q{\bar{q}}$ channel is given by

[TABLE]

where $(z_{a}\leftrightarrow z_{b})$ indicates all previous expressions in the curly brackets repeated with $z_{a}$ and $z_{b}$ interchanged. For the $qg$ channel we have

[TABLE]

The $gq$ channel is given by $\sigma_{gq}^{(1)}(Q,\mu,z_{a},z_{b})=\sigma_{qg}^{(1)}(Q,\mu,z_{b},z_{a})$. The results for $q\leftrightarrow{\bar{q}}$ are identical.

V.3.2 Results for gluon-fusion Higgs production

For gluon-fusion Higgs production, $gg\to H$, we use the effective Lagrangian in the limit $m_{H}^{2}\ll 4m_{t}^{2}$,

[TABLE]

As in Ref. Anastasiou et al. (2003a), the Wilson coefficient $C_{t}$ is always perturbatively expanded against other fixed-order contributions. The Born cross section and LO rapidity spectrum are given by

[TABLE]

We write the partonic cross section as

[TABLE]

The NLO coefficient function for the $gg$ channel is given by

[TABLE]

For the $gq$ and $qg$ channels we find, with $\eta_{qg}^{(1)}(z_{a},z_{b},m_{H},\mu)=\eta_{gq}^{(1)}(z_{b},z_{a},m_{H},\mu)$,

[TABLE]

The $q\bar{q}$ channel is fully regular and given by

[TABLE]

The $N_{c}$-dependent prefactor accounts for the different color average compared to the Born cross section.

The above expressions are in full agreement with Refs. Mathews et al. (2005); Ravindran et al. (2007), which in turn agree with the earliest result for the NLO Drell-Yan rapidity spectrum in Ref. Kubar et al. (1980). In Refs. Mathews et al. (2005); Ravindran et al. (2007), the cross section was also parametrized in terms of $x_{a,b}$ and $z_{a,b}$, but all subtractions were written out in full at the level of the hadronic cross section.

In Ref. Anastasiou et al. (2003b), only the sum of the $gq$ and $qg$ coefficient functions for Higgs productions was given. The separation of the singular terms into the two channels is unique because only the $gq$ ($qg$) channel can be singular as $y\to 0$ ($y\to 1$). We determined the separation of the regular terms by comparing to Ref. Ravindran et al. (2007). To the best of our knowledge, this is the first time that the explicit analytic agreement between the independent NLO calculations in terms of $(z_{a},z_{b})$ and $(z,y)$ has been established.

Finally, we have also compared our numerical implementation of Eq. (S16)-(S24) with the rapidity spectra obtained from Vrap 0.9 Anastasiou et al. (2004) for Drell-Yan and from SusHi 1.7.0 Harlander et al. (2013, 2017) for Higgs production, finding excellent agreement. Since Vrap 0.9 implements the $(z,y)$ parametrization, this effectively confirms the distributional identities in Table S2 numerically with physical PDFs as test functions.

V.4 Comment on rapidity-dependent soft threshold results in the literature

As we have discussed, the soft threshold limit is fully contained in the generalized threshold limit. Our results thus provide an independent confirmation that Eq. (4) is the correct soft threshold factorization for the cross section differential in both $Q$ and $Y$, or equivalently Eq. (5) for the two-dimensional partonic cross section in $(z_{a},z_{b})$.

Several results in the literature Bolzoni (2006); Mukherjee and Vogelsang (2006); Becher et al. (2008); Bonvini et al. (2011, 2015) considering the soft threshold factorization differential in rapidity differ from Eq. (4). The difference is manifest already at fixed NLO in the term $\mathcal{L}_{0}(1-z)[\mathcal{L}_{0}(y)+\mathcal{L}_{0}(1-y)]$, which appears in the flavor-diagonal partonic cross sections. The distributional identity in Eq. (V.2) unambiguously shows that this term has a double singularity in the limit $z_{a}\to 1$ and $z_{b}\to 1$, which means it contributes a priori at leading power in the soft limit $z=z_{a}z_{b}\to 1$, i.e., it contributes to the $m_{a}=m_{b}=-1$ term in Eq. (21). This can already be seen just by considering the distribution in the bulk since

[TABLE]

Moreover, Eq. (V.2) shows that it contributes at leading-logarithmic order. The soft function in Eqs. (4) and (5) precisely contains the leading-power contribution of this term, which is given by the first four terms on the right-hand side of Eq. (V.2).

By contrast, this term and analogous ones at higher order are missing in the leading-power resummed results in Refs. Bolzoni (2006); Mukherjee and Vogelsang (2006); Becher et al. (2008); Bonvini et al. (2011, 2015). There, it is effectively argued that the contribution of such terms to the rapidity spectrum is power-suppressed in $1-z$, leading to the incorrect conclusion that the rapidity dependence in the soft threshold limit can be included simply by taking $\hat{\sigma}_{ij}(z_{a},z_{b})$ to be $\hat{\sigma}_{ij}(z)\,[\delta(y)+\delta(1-y)]/2$ or, depending on the reference, $\hat{\sigma}_{ij}(z)\,\delta(y-1/2)$, where $\hat{\sigma}_{ij}(z)$ is the inclusive, rapidity-integrated, partonic cross section. In the following, we give a critical appraisal of the arguments used to support this conclusion and show why they are flawed.

This replacement first appeared in Ref. Laenen and Sterman (1992), where it was conjectured to provide an approximation to the threshold-resummed rapidity spectrum at small $Y$. The phenomenological impact of the correct convolution structure on PDF determinations relying on soft threshold resummation was discussed in Ref. Westmark and Owens (2017). A detailed numerical study of the difference at the level of the resummed Drell-Yan rapidity spectrum was performed in Ref. Banerjee et al. (2018b).

Argument based on PDF momentum fractions

What makes the $\mathcal{L}_{0}(1-z)[\mathcal{L}_{0}(y)+\mathcal{L}_{0}(1-y)]$ term subtle is that it vanishes upon integration over $y$, so it drops out in the inclusive cross section. This fact alone is of course insufficient to argue that it is power suppressed at each point in the spectrum. It simply means that different leading-power terms conspire to cancel upon integration, which is clear in terms of $(z_{a},z_{b})$, as discussed below Eq. (V.2).

The argument in Ref. Becher et al. (2008) rests on the observation that the PDF arguments, $x_{a}/z_{a}(z,y)$ and $x_{b}/z_{b}(z,y)$, in the two-dimensional convolution integral are independent of $y$ at $z=1$, from which it is concluded that the $y$ dependence of the PDF arguments is power suppressed in $1-z$ and can be dropped. If this is done, the $y$ integral becomes unconstrained and can be carried out freely, which eliminates this term. More generally, one could then replace $\hat{\sigma}_{ij}(z,y)=\hat{\sigma}_{ij}(z)\,\delta(y-1/2)[1+\mathcal{O}(1-z)]$ underneath the convolution integral.

However, a closer inspection of the PDF arguments in terms of $(z,y)$ reveals that

[TABLE]

Hence, the $y$ dependence is not power suppressed but multiplies the leading dependence of the PDF arguments on $z$ itself, and so it cannot be dropped. This becomes even clearer in terms of $(z_{a},z_{b})$,

[TABLE]

which are precisely the $\mathcal{O}(\lambda^{2})$ convolution variables $k^{\mp}/(Qe^{\pm Y})$ in Eq. (4).

To illustrate explicitly that the $y$ dependence has a leading-power effect, consider the hadronic soft threshold limit $1-x_{a}\sim 1-x_{b}\ll 1$ and a simple toy PDF with a power-law behavior near the endpoint, with $\alpha>0$,

[TABLE]

Using Eq. (V.2), it is straightforward to show that $\mathcal{L}_{0}(1-z)[\mathcal{L}_{0}(y)+\mathcal{L}_{0}(1-y)]$ gives rise to double logarithms of $1-x_{a,b}$, but performing the integral directly in terms of $(z,y)$ is tedious, essentially as tedious as deriving Eq. (V.2) itself. Instead, to disprove the above argument and show that the $y$ dependence is not power suppressed, it suffices to consider two terms that have the same $y$ integral,

[TABLE]

and show that they give different results at leading power, while the above argument would imply that they do not. Convolving $A(z,y)$ and $B(z,y)$ against the toy PDFs over the domain shown in Fig. S1 yields

[TABLE]

where $H_{\alpha}$ is the harmonic number. To evaluate the integrals it is convenient to already expand at integrand level, e.g. $1-x_{a}/z_{a}=1-x_{a}+y(1-z)$ up to higher powers in $1-z$ and $1-x_{a}$. The maximum in the second case arises because the integration region in $z$ along fixed $y=1/2$ is cut off by the square of the larger of the two momentum fractions. (The order of expanding in $1-x_{a}$ and $1-x_{b}$ also needs to be picked accordingly.) Clearly, the two results only coincide for $Y=0$, where $x_{a}=x_{b}$. Away from $Y=0$, the logarithmic dependence on $x_{a,b}$ and thus on $Y$ differs at leading power.

Fourier-transform argument

An alternative line of argument Bolzoni (2006); Mukherjee and Vogelsang (2006); Bonvini et al. (2011) relies on taking the Fourier transform of the partonic cross section to also argue that the $y$ dependence is trivial. A first step is to change variables from $y$ to

[TABLE]

(Note that the variables $u$ and $y$ are precisely interchanged in the notation used in Ref. Bonvini et al. (2011).) One then considers the Fourier transform of the partonic cross section $C(z,u)$ with respect to $u$,

[TABLE]

The second equality, which is in question, is based on observing that $C(z,u)$ only has support on an interval bounded by $u_{\mathrm{max}}\sim 1-z$, and concluding that the Fourier kernel can be expanded in $u\sim 1-z$ as $e^{\mathrm{i}Mu}\stackrel{{\scriptstyle?}}{{=}}1+\mathcal{O}(1-z)$, and so $\tilde{C}(z,M)$ is independent of $M$ at leading power. Thus, taking the inverse Fourier transform, the partonic cross section may be approximated as

[TABLE]

This argument is flawed because in order to satisfy the Fourier inversion theorem, one must count $M\sim(1-z)^{-1}$ if one wants to count $u\sim 1-z$. In particular, one is not allowed to count $M\sim 1$ when taking the limit $z\to 1$ (or equivalently $N\to\infty$ for the Mellin conjugate $N$ of $z$). This is essential because $C(z,u)$ contains distributional terms in $u$ that cancel the suppression by the integration domain.

To disprove Eq. (S33), it again suffices to consider $A(z,y)$ and $B(z,y)$ defined in Eq. (S29). Changing variables to $u$, we have

[TABLE]

Both terms satisfy the assumptions of the above argument, i.e., they only have support for $\lvert u\rvert\leq u_{\mathrm{max}}$. Changing variables back to $y$, Eq. (S33) would imply that up to power corrections in $1-z$,

[TABLE]

In fact, the overall factor found in Ref. Bonvini et al. (2011) is $\delta(y-1/2)$, while it is $[\delta(y)+\delta(1-y)]/2$ in Ref. Becher et al. (2008), and the above argument was used in Ref. Bonvini et al. (2011) to argue that the two are equivalent. As a distributional identity, this is obviously incorrect. The only thing that is equal between $A(z,y)$ and $B(z,y)$ are their $y$ integrals, and as demonstrated before, this is insufficient because the $y$ dependence of the PDF arguments is a leading-power effect and cannot be neglected.

Summary

For $pp$ production processes in general, to correctly describe the soft threshold limit of differential observables that are sensitive to the total rapidity of the Born system, one must maintain the two-dimensional dependence on $z_{a}$ and $z_{b}$ in the convolutions against the PDFs. Equivalently, in Mellin space one must maintain two Mellin fractions $N_{a}$ and $N_{b}$ as in the original Ref. Catani and Trentadue (1989). In terms of the Mellin conjugate $N$ of $z$ and a Fourier conjugate $M$ of another variable like $u$, one has to keep the dependence on $M\sim\lvert N_{a}-N_{b}\rvert\sim N$. In particular, reducing the two-dimensional convolution structure to one dimension – such that the rapidity dependence is only carried by the luminosity function – amounts to making an additional assumption that is not justified by taking the soft limit.

V.5 Perturbative ingredients

Here, we collect the required perturbative ingredients to evaluate the generalized threshold factorization theorem to two loops as well as the highest few logarithmic terms at three loops.

V.5.1 Hard functions

The hard function for Drell-Yan production, $q\bar{q}\to Z/\gamma^{\ast}\to\ell^{+}\ell^{-}$, is given by

[TABLE]

where the sum runs over $q=\{u,d,c,s,b\}$, $\sigma_{B,q}^{\mathrm{DY}}$ is the Born cross section given in Eq. (S14), and $C_{q{\bar{q}}}^{V}$ is the Wilson coefficient from matching the QCD quark vector current onto SCET. In principle, we also need the matching coefficient $C_{q{\bar{q}}}^{A}$ for the axial-vector current, which differs from $C_{q{\bar{q}}}^{V}$ starting at $\mathcal{O}(\alpha_{s}^{2})$ by small singlet corrections due to the large mass splitting between bottom and top quarks. We neglect these terms, as is done in Vrap 0.9, and use $H_{ij}^{\mathrm{DY}}\propto|C_{q{\bar{q}}}^{V}|^{2}$ throughout. The hard function for gluon-fusion Higgs production in the limit $m_{H}^{2}\ll 4m_{t}^{2}$ reads

[TABLE]

where $\sigma_{B}^{ggH}$ and $C_{t}$ are given in Eqs. (S19) and (S20). The Wilson coefficient $C_{gg}$ arises from matching the $gg\to H$ operator in Eq. (S19) onto SCET. The Wilson coefficients contain the IR-finite virtual corrections to the Born process. They are normalized as $C=1+\mathcal{O}(\alpha_{s})$ and can be found in Ref. Ebert et al. (2017) in our notation.

In the main text, we also refer to the fixed-order expansion of the hard function, where the coefficients $H_{ij}^{(n)}$ include all prefactors that are present at Born level,

[TABLE]

V.5.2 Beam functions

For $t\gg\Lambda_{\mathrm{QCD}}^{2}$, the modified beam function defined by Eq. (15) can be matched onto PDFs as

[TABLE]

which directly follows from the analogous matching relations for the inclusive and double-differential beam functions, $B_{i}(t,x,\mu)$ and $B_{i}(t,\vec{\kern 0.0ptk}_{T},x,\mu)$, with matching coefficients $\mathcal{I}_{ij}(t,z,\mu)$ and $\mathcal{I}_{ij}(t,\vec{\kern 0.0ptk}_{T},z,\mu)$ Stewart et al. (2010a, c); Jain et al. (2012).

The modified beam function satisfies the same RGE as the other beam functions,

[TABLE]

because the RGE does not change the $\vec{\kern 0.0ptk}_{T}$ dependence Jain et al. (2012). The matching coefficients satisfy the RGE Stewart et al. (2010c)

[TABLE]

where $P_{ij}(\alpha_{s},z)$ is the PDF anomalous dimension. We define the perturbative expansion of $\tilde{\mathcal{I}}_{ij}$ as

[TABLE]

Solving Eq. (S41) order by order in $\alpha_{s}$, we obtain a general expression for $\tilde{\mathcal{I}}_{ij}^{(n)}(t,z,\mu)$

[TABLE]

where $\mathcal{L}_{n}(t,\mu^{2})\equiv(1/\mu^{2})\mathcal{L}_{n}(t/\mu^{2})$, $\tilde{I}^{(n)}_{ij}(z)$ is the $\mathcal{O}(\alpha_{s}^{n})$ boundary term that is not predicted by the RGE, and the ellipses in the three-loop result indicate the terms proportional to $\mathcal{L}_{0,1,2}(t,\mu^{2})$. We also introduced the shorthand $(g\otimes h)(z)\equiv\int_{z}^{1}\!\mathrm{d}z^{\prime}/z^{\prime}\,g(z^{\prime})\,h(z/z^{\prime})$ for the Mellin convolution of two functions of $z$. In practice, we use the MT package Höschele et al. (2014) to evaluate them analytically. Expanding the three-loop expression for $\tilde{\mathcal{I}}_{ij}(t,z,\mu)$ against the hard function yields Eq. (22) in the main text.

The anomalous dimension coefficients in Eq. (V.5.2) are as follows. The QCD $\beta$ function coefficients are

[TABLE]

The anomalous dimensions are expanded as

[TABLE]

The coefficients of the cusp anomalous dimension to two loops are Korchemsky and Radyushkin (1987)

[TABLE]

where $C_{q}=C_{F}$ and $C_{g}=C_{A}$. The beam function anomalous dimension coefficients are Stewart et al. (2010c); Berger et al. (2011); Gaunt et al. (2014a, b)

[TABLE]

The one-loop PDF anomalous dimensions are

[TABLE]

with the standard color-stripped splitting functions,

[TABLE]

V.5.3 Calculation of beam function boundary terms

The structure of Eq. (V.5.2) is exactly the same as for the inclusive beam function $\mathcal{I}_{ij}(t,z,\mu)$ Gaunt et al. (2014a, b), except for the different boundary terms $I^{(n)}_{ij}(z)\neq\tilde{I}^{(n)}_{ij}(z)$. The definition in Eq. (15) implies for the matching coefficients

[TABLE]

which we use to calculate $\tilde{I}^{(n)}_{ij}(z)$ to the extent that the double-differential matching coefficients are known, i.e., to one loop for $i=g$ Jain et al. (2012) and two loops for $i=q$ Jain et al. (2012); Gaunt and Stahlhofen (2014). Note that for $i=g$, the integral over all $\vec{\kern 0.0ptk}_{T}$ leaves behind only the polarization-independent piece of the double-differential gluon beam function. We have also verified that the $\mu$-dependent pieces $\propto\mathcal{L}_{n}(t,\mu)$ obtained from Eq. (S49) agree with the RGE prediction Eq. (V.5.2), i.e., we have explicitly checked that the projection and the RGE commute.

At one loop, we find that the following simple relation holds for all partonic channels,

[TABLE]

where $\mathcal{I}_{ij}^{(1)}$ is the one-loop matching coefficient for the inclusive beam function. Explicitly, the one-loop finite terms of the modified beam function are given by

[TABLE]

We decompose the two-loop quark finite terms $\tilde{I}_{ij}^{(2)}(z)$ by their flavor structure as

[TABLE]

As was done for $I^{(2)}_{ij}(z)$ in Refs. Gaunt et al. (2014a, b), we find it convenient to pull common rational factors out of recurring terms with transcendental weight three ($\tilde{S}_{3},\tilde{T}_{3},\tilde{U}_{3},\tilde{V}_{3},\tilde{R}_{\dots}$), and group terms of lower transcendental weight separately by color factor and flavor structure ($\tilde{C}_{\dots}$). We also pull out a conventional factor of four:

[TABLE]

Here, overall factors of $\theta(z)\,\theta(1-z)$ are understood, but omitted for brevity. The $D_{qqV,\dots}(z)$ contain all distributional terms in $1-z$ and are the same as for the standard inclusive beam function Gaunt et al. (2014a),

[TABLE]

All remaining terms in Eq. (V.5.3) are integrable for $z\to 1$. Their full expressions are lengthy and are available from the authors upon request. As an example of the structures that occur, we give

[TABLE]

Here, we have used the recent PolyLogTools package Duhr and Dulat (2019) to convert all polylogarithms of rational functions of $z$ to standard harmonic polylogarithms $H(a_{1},\dots,a_{n};z)$ as well as multiple polylogarithms $G(a_{1},\dots,a_{n};z)$ of $z\in[0,1]$. The latter are as defined in Ref. Duhr and Dulat (2019), and for all $a_{i}=0,\pm 1$ reduce to standard harmonic polylogarithms up to a sign. We find no evidence for a simple generalization of the one-loop relation Eq. (S50) at two loops. Finally, we note that the two-loop $\tilde{I}_{ij}^{(2)}(z)$ for the modified beam function are substantially more complicated than those for the standard inclusive beam function. For example, the latter does not involve polylogarithms with fractional weights, which only arise from the projection integral in Eq. (S49). We expect that similarly the three-loop cross section in terms of $(q^{+},q^{-},\vec{\kern 0.0ptq}_{T})$ will have a much simpler structure than in terms of $(Q,Y,\vec{\kern 0.0ptq}_{T})$.

V.6 Breakdown of NNLO validation into partonic channels

Here, we provide the breakdown of the numerical validation of Eq. (14) at NNLO in Fig. 3 into individual partonic channels. Following Vrap Anastasiou et al. (2004), we take the $ij=q\bar{q}$ channel to include all topologies where $i$ and $j$ are part of the same quark line. The leading-power limit of these diagrams corresponds to the $qqV$ beam function matching coefficient in the decomposition in Eq. (V.5.3). In addition, the $q\bar{q}$ channel also includes purely nonsingular contributions with topologies $q\bar{q}\to g\to q\bar{q}V$. We then take the $qq^{\prime}$ channel to include the remaining quark-initiated processes, which at leading power reduces to the sum of the $qqS$ and $q\bar{q}V$ beam function contributions in Eq. (V.5.3). The $qg$ channel maps onto the $qg$ beam function contribution at leading power, while the $gq$ and $gg$ channels are purely nonsingular.

The results are shown in Fig. S2. In all cases, the prediction of Eq. (14) is in excellent agreement with the singular limit of the full calculation, with their difference vanishing as a power of $1-x_{a}$ as it should. The excellent numerical stability of Vrap for the nondiagonal channels allows us to extend the check down to $1-x_{a}=10^{-5}$, where it becomes limited by MC statistics. For the $q\bar{q}$ channel, we start to see a systematic deviation at the $10^{-4}$ level below $1-x_{a}\lesssim 10^{-4}$. We observe a similar deviation already at NLO, where the partonic cross sections agree analytically. We thus attribute this to a systematic effect in the PDF integrations in Vrap.

V.7 Results for Drell-Yan at $\mu=Q/2$ and for gluon-fusion Higgs production

Here, we provide additional results for the generalized threshold expansion at fixed order. The analogs of Figs. 4 and 5 for the Drell-Yan rapidity spectrum at a different scale choice $\mu=Q/2$ are shown in Fig. S3, and for the $gg\to H$ NLO rapidity spectrum at $\mu=m_{H}$ and $\mu=m_{H}/2$ in Fig. S4.

When performing any threshold expansion for different scale choices, there are several options how to treat the terms in the partonic cross section that are predicted by the running of the PDFs or $\alpha_{s}$. One option is to expand these terms to the working order in the threshold expansion. This ensures that the partonic cross section has homogeneous power counting at any scale, but leaves the running of the PDFs and the coupling uncanceled beyond the working order. Another option is to threshold-expand the partonic cross section at a given reference scale and treat the running exactly. This leads to a privileged scale where the expansion was performed, but ensures the cancellation of $\alpha_{s}$ and PDF running to all powers (up to higher orders in $\alpha_{s}$). In the following results, we choose the first option for definiteness. The difference between the two approaches could serve as a way to estimate the size of power corrections.

For Drell-Yan at $\mu=Q/2$, the generalized threshold expansion performs similarly well as for $\mu=Q$ in the main text, and again much better than the soft expansion. For $gg\to H$, the generalized threshold expansion again performs in a manner clearly superior to the soft one. Here, the increment from the leading-power soft to the leading-power generalized approximation of the $gg$ channel at NLO is roughly comparable to the piece still missing to the full result; either contribution amounts to $\mathcal{O}(20\%)$ in units of the Born cross section. This is consistent with the expectation that for a gluon-induced process, hard central radiation plays a larger role than for Drell-Yan. The shape of the NLO contributions at large $Y$ is well captured by the leading-power generalized approximation for both the $gg$ and $gq+qg$ channel. The leading-power soft approximation for the $gg$ channel (on close inspection) turns out to be off at large $Y$, and in both channels there is barely any convergence beyond leading power in the soft expansion at any $Y$.