# Generalized Threshold Factorization with Full Collinear Dynamics

Gillian Lustermans<EMAIL_ADDRESS>Institute for Theoretical Physics
Amsterdam and Delta Institute for Theoretical Physics, University of
Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands  Nikhef,
Theory Group, Science Park 105, 1098 XG, Amsterdam, The Netherlands  Johannes
K. L. Michel<EMAIL_ADDRESS>Theory Group, Deutsches Elektronen-
Synchrotron (DESY), D-22607 Hamburg, Germany  Frank J. Tackmann
<EMAIL_ADDRESS>Theory Group, Deutsches Elektronen-Synchrotron (DESY),
D-22607 Hamburg, Germany

(August 2, 2019)

###### 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.

††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$,

$\displaystyle p(P_{a}^{\mu})+p(P_{b}^{\mu})\to L(q^{\mu})+X(P_{X}^{\mu})\,,$
(1)

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

$\displaystyle x_{a}=\frac{Q}{E_{\mathrm{cm}}}e^{+Y}\,,\quad
x_{b}=\frac{Q}{E_{\mathrm{cm}}}e^{-Y}\,,\quad\tau=x_{a}x_{b}\,,$ (2)

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)

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}x_{a}\mathrm{d}x_{b}}$
$\displaystyle=\int\\!\frac{\mathrm{d}z_{a}}{z_{a}}\frac{\mathrm{d}z_{b}}{z_{b}}\,\hat{\sigma}_{ij}(z_{a},z_{b})\,f_{i}\Bigl{(}\frac{x_{a}}{z_{a}}\Bigr{)}\,f_{j}\Bigl{(}\frac{x_{b}}{z_{b}}\Bigr{)}\,,$
(3)

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),

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}x_{a}\mathrm{d}x_{b}}$
$\displaystyle=H_{ij}(Q^{2})\int\\!\mathrm{d}k^{-}\,\mathrm{d}k^{+}S(k^{-},k^{+})$
(4) $\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{]}.$

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$,

$\displaystyle\hat{\sigma}_{ij}(z_{a},z_{b})$
$\displaystyle=H_{ij}\,\hat{S}(z_{a},z_{b})\,,$ (5)

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

$\displaystyle
x_{\mp}\equiv\frac{q^{\mp}}{P_{a,b}^{\mp}}=\frac{\sqrt{Q^{2}+q_{T}^{2}}}{E_{\mathrm{cm}}}\,e^{\pm
Y}\,.$ (6)

We consider the generalized threshold limit

$\lambda_{\mathrm{QCD}}^{2}\ll\lambda^{2}\sim 1-x_{-}\ll 1\quad\text{for
generic $x_{+}$}\,,$ (7)

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

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

Figure 1: Illustration of Drell-Yan at large dilepton rapidity.

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

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}x_{-}\mathrm{d}x_{+}}$
$\displaystyle=H_{ij}(q^{+}q^{-})\int\\!\mathrm{d}t\,f_{i}^{\mathrm{thr}}\Bigl{[}x_{-}\Bigl{(}1+\frac{t}{q^{+}q^{-}}\Bigr{)}\Bigr{]}$
$\displaystyle\quad\times B_{j}(t,x_{+})\,.$ (9)

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)

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

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,

$\displaystyle p_{\bar{n}}$ $\displaystyle\sim
Q\,\bigl{(}1,\lambda^{2},\lambda\bigr{)}\,,\quad P_{\bar{n}}\sim
Q\,(1,\lambda_{\mathrm{QCD}}^{2},\lambda_{\mathrm{QCD}})\,,$ (11)
$\displaystyle P_{n}$ $\displaystyle\sim
Q\,(\lambda_{\mathrm{QCD}}^{2},1,\lambda_{\mathrm{QCD}})\,,\quad P_{s}\sim
Q\,\Bigl{(}\frac{\lambda_{\mathrm{QCD}}^{2}}{\lambda^{2}},\lambda^{2},\lambda_{\mathrm{QCD}}\Bigr{)}\,.$

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

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}x_{-}\mathrm{d}x_{+}\,\mathrm{d}\vec{\kern
0.0ptq}_{T}}$
$\displaystyle=H_{ij}(q^{+}q^{-})\int\\!\mathrm{d}t\,f_{i}^{\mathrm{thr}}\Bigl{[}x_{-}\Bigl{(}1+\frac{t}{q^{+}q^{-}}\Bigr{)}\Bigr{]}$
$\displaystyle\quad\times B_{j}(t,\vec{\kern 0.0ptq}_{T},x_{+})\,.$ (12)

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 $0$-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

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}x_{a}\mathrm{d}x_{b}\,\mathrm{d}\vec{\kern
0.0ptq}_{T}}$
$\displaystyle=H_{ij}(Q^{2})\int\\!\mathrm{d}t\,f_{i}^{\mathrm{thr}}\Bigl{[}x_{a}\Bigl{(}1+\frac{q_{T}^{2}}{2Q^{2}}+\frac{t}{Q^{2}}\Bigr{)}\Bigr{]}$
$\displaystyle\quad\times B_{j}(t,\vec{\kern 0.0ptq}_{T},x_{b})\,.$ (13)

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

$\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})\,.$
(14)

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

$\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{)}\,.$ (15)

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

$\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{]}$
(16)

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

$\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{]}\,,$ (17)

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

$\displaystyle\hat{\sigma}_{ij}(z_{a},z_{b})$
$\displaystyle=H_{k\ell}\bigl{[}\delta_{ki}\,\hat{\mathcal{I}}_{\ell
j}(z_{a},z_{b})+\hat{\mathcal{I}}_{ki}(z_{b},z_{a})\,\delta_{\ell j}$
$\displaystyle\quad\quad\quad-\delta_{ki}\,\delta_{\ell
j}\,\hat{S}(z_{a},z_{b})\bigr{]}\,,$ (18)

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

$\displaystyle\hat{\mathcal{I}}_{ij}(z_{a},z_{b})$ $\displaystyle\equiv
Q^{2}\,\tilde{\mathcal{I}}_{ij}\bigl{[}Q^{2}(1-z_{a}),z_{b}\bigr{]}\,,$
$\displaystyle\hat{S}(z_{a},z_{b})$ $\displaystyle\equiv
Q^{2}\,S\bigl{[}Qe^{Y}(1-z_{a}),Qe^{-Y}(1-z_{b})\bigr{]}\,.$ (19)

Figure 2: Terms in the partonic cross section $\hat{\sigma}(z_{a},z_{b})$
captured by the soft and generalized threshold expansions.

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

Figure 3: Validation of the $\mathcal{O}(\alpha_{s}^{2})$ contribution to
$\mathrm{d}\sigma/\mathrm{d}x_{a}\mathrm{d}x_{b}$ predicted for $x_{a}\to 1$
by Eq. (14) (blue) against the full result from Vrap (red). Their difference
(green) vanishes like a power as $1-x_{a}\to 0$, as it must. The error bars
indicate the integration uncertainties.

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

$\displaystyle\hat{\mathcal{I}}_{ij}(z_{a},z_{b})$
$\displaystyle=\delta_{ij}\,\hat{S}(z_{a},z_{b})\,\bigl{[}1+\mathcal{O}(1-z_{b})\bigr{]}\,,$
(20)

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.

Figure 4: The $\mathcal{O}(\alpha_{s})$ (top) and
$\mathcal{O}(\alpha_{s}^{2})$ (bottom) contributions to
$\sigma\equiv\mathrm{d}\sigma/\mathrm{d}Q\mathrm{d}Y$ normalized to the LO
result. Shown are the full result (red), the generalized threshold
approximation (blue), and the soft threshold approximation (gray).

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,

$\hat{\sigma}_{ij}(z_{a},z_{b})=\sum_{m_{a},m_{b}}\hat{\sigma}_{ij}^{(m_{a},m_{b})}(z_{a},z_{b})\,,$
(21)

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}]$.

Figure 5: Convergence of the generalized (blue) and soft (gray) threshold
expansions. Shown are the deviations from the full NLO result for the
$q\bar{q}$ channel (top) and the $qg+gq$ channel (bottom) normalized to the LO
result.

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

$\displaystyle\hat{\sigma}_{ij}^{(3)}$
$\displaystyle=\mathcal{L}_{5}(1-z_{a})H^{(0)}_{ij}\delta(1-z_{b})\,\frac{(\Gamma_{0}^{i})^{3}}{8}$
(22)
$\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\,.$

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).

## References

  * Bodwin (1985) G. T. Bodwin, Phys. Rev. D31, 2616 (1985), [Erratum: Phys. Rev.D34,3932(1986)].
  * Collins  _et al._ (1985) J. C. Collins, D. E. Soper, and G. F. Sterman, Nucl. Phys. B261, 104 (1985).
  * Collins  _et al._ (1988) J. C. Collins, D. E. Soper, and G. F. Sterman, Nucl. Phys. B308, 833 (1988).
  * Sterman (1987) G. F. Sterman, Nucl. Phys. B281, 310 (1987).
  * Catani and Trentadue (1989) S. Catani and L. Trentadue, Nucl. Phys. B327, 323 (1989).
  * Ravindran  _et al._ (2007) V. Ravindran, J. Smith, and W. L. van Neerven, Nucl. Phys. B767, 100 (2007), arXiv:hep-ph/0608308 .
  * Westmark and Owens (2017) D. Westmark and J. F. Owens, Phys. Rev. D95, 056024 (2017), arXiv:1701.06716 [hep-ph] .
  * Banerjee  _et al._ (2018a) P. Banerjee, G. Das, P. K. Dhani, and V. Ravindran, Phys. Rev. D97, 054024 (2018a), arXiv:1708.05706 [hep-ph] .
  * Banerjee  _et al._ (2018b) P. Banerjee, G. Das, P. K. Dhani, and V. Ravindran, Phys. Rev. D98, 054018 (2018b), arXiv:1805.01186 [hep-ph] .
  * Appell  _et al._ (1988) D. Appell, G. F. Sterman, and P. B. Mackenzie, Nucl. Phys. B309, 259 (1988).
  * Magnea (1991) L. Magnea, Nucl. Phys. B349, 703 (1991).
  * Korchemsky and Marchesini (1993) G. P. Korchemsky and G. Marchesini, Nucl. Phys. B406, 225 (1993), arXiv:hep-ph/9210281 .
  * Contopanagos  _et al._ (1997) H. Contopanagos, E. Laenen, and G. F. Sterman, Nucl. Phys. B484, 303 (1997), arXiv:hep-ph/9604313 .
  * Catani  _et al._ (1996) S. Catani, M. L. Mangano, P. Nason, and L. Trentadue, Nucl. Phys. B478, 273 (1996), arXiv:hep-ph/9604351 .
  * Belitsky (1998) A. V. Belitsky, Phys. Lett. B442, 307 (1998), arXiv:hep-ph/9808389 .
  * Moch and Vogt (2005) S. Moch and A. Vogt, Phys. Lett. B631, 48 (2005), arXiv:hep-ph/0508265 .
  * Laenen and Magnea (2006) E. Laenen and L. Magnea, Phys. Lett. B632, 270 (2006), arXiv:hep-ph/0508284 .
  * Idilbi  _et al._ (2006) A. Idilbi, X.-d. Ji, and F. Yuan, Nucl. Phys. B753, 42 (2006), arXiv:hep-ph/0605068 .
  * Mukherjee and Vogelsang (2006) A. Mukherjee and W. Vogelsang, Phys. Rev. D73, 074005 (2006), arXiv:hep-ph/0601162 .
  * Bolzoni (2006) P. Bolzoni, Phys. Lett. B643, 325 (2006), arXiv:hep-ph/0609073 .
  * Becher  _et al._ (2008) T. Becher, M. Neubert, and G. Xu, JHEP 07, 030 (2008), arXiv:0710.0680 [hep-ph] .
  * Bonvini  _et al._ (2011) M. Bonvini, S. Forte, and G. Ridolfi, Nucl. Phys. B847, 93 (2011), arXiv:1009.5691 [hep-ph] .
  * Bonvini  _et al._ (2015) M. Bonvini, S. Marzani, J. Rojo, L. Rottoli, M. Ubiali, R. D. Ball, V. Bertone, S. Carrazza, and N. P. Hartland, JHEP 09, 191 (2015), arXiv:1507.01006 [hep-ph] .
  * Fuks  _et al._ (2013) B. Fuks, M. Klasen, D. R. Lamprea, and M. Rothering, Eur. Phys. J. C73, 2480 (2013), arXiv:1304.0790 [hep-ph] .
  * Bonvini and Marzani (2014) M. Bonvini and S. Marzani, JHEP 09, 007 (2014), arXiv:1405.3654 [hep-ph] .
  * Schmidt and Spira (2016) T. Schmidt and M. Spira, Phys. Rev. D93, 014022 (2016), arXiv:1509.00195 [hep-ph] .
  * A H  _et al._ (2019) A. A H, A. Chakraborty, G. Das, P. Mukherjee, and V. Ravindran, (2019), arXiv:1905.03771 [hep-ph] .
  * Bonocore  _et al._ (2016) D. Bonocore, E. Laenen, L. Magnea, L. Vernazza, and C. D. White, JHEP 12, 121 (2016), arXiv:1610.06842 [hep-ph] .
  * Del Duca  _et al._ (2017) V. Del Duca, E. Laenen, L. Magnea, L. Vernazza, and C. D. White, JHEP 11, 057 (2017), arXiv:1706.04018 [hep-ph] .
  * Beneke  _et al._ (2019) M. Beneke, A. Broggio, M. Garny, S. Jaskiewicz, R. Szafron, L. Vernazza, and J. Wang, JHEP 03, 043 (2019), arXiv:1809.10631 [hep-ph] .
  * Anastasiou  _et al._ (2014) C. Anastasiou, C. Duhr, F. Dulat, E. Furlan, T. Gehrmann, F. Herzog, and B. Mistlberger, Phys. Lett. B737, 325 (2014), arXiv:1403.4616 [hep-ph] .
  * Ahmed  _et al._ (2014a) T. Ahmed, M. Mahakhud, N. Rana, and V. Ravindran, Phys. Rev. Lett. 113, 112002 (2014a), arXiv:1404.0366 [hep-ph] .
  * Ahmed  _et al._ (2014b) T. Ahmed, M. K. Mandal, N. Rana, and V. Ravindran, Phys. Rev. Lett. 113, 212003 (2014b), arXiv:1404.6504 [hep-ph] .
  * de Florian  _et al._ (2014) D. de Florian, J. Mazzitelli, S. Moch, and A. Vogt, JHEP 10, 176 (2014), arXiv:1408.6277 [hep-ph] .
  * Li  _et al._ (2015) Y. Li, A. von Manteuffel, R. M. Schabinger, and H. X. Zhu, Phys. Rev. D91, 036008 (2015), arXiv:1412.2771 [hep-ph] .
  * Anastasiou  _et al._ (2016) C. Anastasiou, C. Duhr, F. Dulat, E. Furlan, T. Gehrmann, F. Herzog, A. Lazopoulos, and B. Mistlberger, JHEP 05, 058 (2016), arXiv:1602.00695 [hep-ph] .
  * Dulat  _et al._ (2018) F. Dulat, B. Mistlberger, and A. Pelloni, JHEP 01, 145 (2018), arXiv:1710.03016 [hep-ph] .
  * Dulat  _et al._ (2019) F. Dulat, B. Mistlberger, and A. Pelloni, Phys. Rev. D99, 034004 (2019), arXiv:1810.09462 [hep-ph] .
  * Stewart  _et al._ (2010a) I. W. Stewart, F. J. Tackmann, and W. J. Waalewijn, Phys. Rev. D81, 094035 (2010a), arXiv:0910.0467 [hep-ph] .
  * Stewart  _et al._ (2010b) I. W. Stewart, F. J. Tackmann, and W. J. Waalewijn, Phys. Rev. Lett. 105, 092002 (2010b), arXiv:1004.2489 [hep-ph] .
  * Stewart  _et al._ (2010c) I. W. Stewart, F. J. Tackmann, and W. J. Waalewijn, JHEP 09, 005 (2010c), arXiv:1002.2213 [hep-ph] .
  * Berger  _et al._ (2011) C. F. Berger, C. Marcantonini, I. W. Stewart, F. J. Tackmann, and W. J. Waalewijn, JHEP 04, 092 (2011), arXiv:1012.4480 [hep-ph] .
  * Gaunt  _et al._ (2014a) J. R. Gaunt, M. Stahlhofen, and F. J. Tackmann, JHEP 04, 113 (2014a), arXiv:1401.5478 [hep-ph] .
  * Gaunt  _et al._ (2014b) J. Gaunt, M. Stahlhofen, and F. J. Tackmann, JHEP 08, 020 (2014b), arXiv:1405.1044 [hep-ph] .
  * Melnikov  _et al._ (2019a) K. Melnikov, R. Rietkerk, L. Tancredi, and C. Wever, JHEP 02, 159 (2019a), arXiv:1809.06300 [hep-ph] .
  * Melnikov  _et al._ (2019b) K. Melnikov, R. Rietkerk, L. Tancredi, and C. Wever, JHEP 06, 033 (2019b), arXiv:1904.02433 [hep-ph] .
  * Bauer  _et al._ (2000) C. W. Bauer, S. Fleming, and M. E. Luke, Phys. Rev. D63, 014006 (2000), arXiv:hep-ph/0005275 .
  * Bauer  _et al._ (2001) C. W. Bauer, S. Fleming, D. Pirjol, and I. W. Stewart, Phys. Rev. D63, 114020 (2001), arXiv:hep-ph/0011336 .
  * Bauer and Stewart (2001) C. W. Bauer and I. W. Stewart, Phys. Lett. B516, 134 (2001), arXiv:hep-ph/0107001 .
  * Bauer  _et al._ (2002a) C. W. Bauer, D. Pirjol, and I. W. Stewart, Phys. Rev. D65, 054022 (2002a), arXiv:hep-ph/0109045 .
  * Bauer  _et al._ (2002b) C. W. Bauer, S. Fleming, D. Pirjol, I. Z. Rothstein, and I. W. Stewart, Phys. Rev. D66, 014017 (2002b), arXiv:hep-ph/0202088 .
  * (52) .
  * Fleming and Labun (2015) S. Fleming and O. Z. Labun, Phys. Rev. D91, 094011 (2015), arXiv:1210.1508 [hep-ph] .
  * Hoang  _et al._ (2016) A. H. Hoang, P. Pietrulewicz, and D. Samitz, Phys. Rev. D93, 034034 (2016), arXiv:1508.04323 [hep-ph] .
  * Jain  _et al._ (2012) A. Jain, M. Procura, and W. J. Waalewijn, JHEP 04, 132 (2012), arXiv:1110.0839 [hep-ph] .
  * Gaunt and Stahlhofen (2014) J. R. Gaunt and M. Stahlhofen, JHEP 12, 146 (2014), arXiv:1409.8281 [hep-ph] .
  * Procura  _et al._ (2015) M. Procura, W. J. Waalewijn, and L. Zeune, JHEP 02, 117 (2015), arXiv:1410.6483 [hep-ph] .
  * Lustermans  _et al._ (2019) G. Lustermans, J. K. L. Michel, F. J. Tackmann, and W. J. Waalewijn, JHEP 03, 124 (2019), arXiv:1901.03331 [hep-ph] .
  * Manohar (2003) A. V. Manohar, Phys. Rev. D68, 114019 (2003), arXiv:hep-ph/0309176 .
  * Becher  _et al._ (2007) T. Becher, M. Neubert, and B. D. Pecjak, JHEP 01, 076 (2007), arXiv:hep-ph/0607228 .
  * Chay and Kim (2013) J. Chay and C. Kim, JHEP 09, 126 (2013), arXiv:1303.1637 [hep-ph] .
  * Kang  _et al._ (2013) D. Kang, C. Lee, and I. W. Stewart, Phys. Rev. D88, 054004 (2013), arXiv:1303.6952 [hep-ph] .
  * Anastasiou  _et al._ (2003a) C. Anastasiou, L. J. Dixon, and K. Melnikov, Nucl. Phys. Proc. Suppl. 116, 193 (2003a), hep-ph/0211141 .
  * Anastasiou  _et al._ (2003b) C. Anastasiou, L. J. Dixon, K. Melnikov, and F. Petriello, Phys. Rev. Lett. 91, 182002 (2003b), arXiv:hep-ph/0306192 .
  * Kubar  _et al._ (1980) J. Kubar, M. Le Bellac, J. L. Meunier, and G. Plaut, Nucl. Phys. B175, 251 (1980).
  * Mathews  _et al._ (2005) P. Mathews, V. Ravindran, K. Sridhar, and W. L. van Neerven, Nucl. Phys. B713, 333 (2005), arXiv:hep-ph/0411018 .
  * Ebert  _et al._ (2018) M. A. Ebert, J. K. L. Michel, F. J. Tackmann, _et al._ , DESY-17-099 (2018), webpage: http://scetlib.desy.de.
  * Anastasiou  _et al._ (2004) C. Anastasiou, L. J. Dixon, K. Melnikov, and F. Petriello, Phys. Rev. D69, 094008 (2004), arXiv:hep-ph/0312266 .
  * Harland-Lang  _et al._ (2015) L. A. Harland-Lang, A. D. Martin, P. Motylinski, and R. S. Thorne, Eur. Phys. J. C75, 204 (2015), arXiv:1412.3989 [hep-ph] .
  * (70) G. Billis, M. A. Ebert, J. K. L. Michel, and F. J. Tackmann, (in preparation).
  * Bauer  _et al._ (2012) C. W. Bauer, F. J. Tackmann, J. R. Walsh, and S. Zuberi, Phys. Rev. D85, 074006 (2012), arXiv:1106.6047 [hep-ph] .
  * Manohar  _et al._ (2002) A. V. Manohar, T. Mehen, D. Pirjol, and I. W. Stewart, Phys. Lett. B539, 59 (2002), arXiv:hep-ph/0204229 .
  * Harlander  _et al._ (2013) R. V. Harlander, S. Liebler, and H. Mantler, Comput. Phys. Commun. 184, 1605 (2013), arXiv:1212.3249 [hep-ph] .
  * Harlander  _et al._ (2017) R. V. Harlander, S. Liebler, and H. Mantler, Comput. Phys. Commun. 212, 239 (2017), arXiv:1605.03190 [hep-ph] .
  * Laenen and Sterman (1992) E. Laenen and G. F. Sterman, in  _The Fermilab Meeting DPF 92. Proceedings, 7th Meeting of the American Physical Society, Division of Particles and Fields, Batavia, USA, November 10-14, 1992. Vol. 1, 2_ (1992) pp. 987–989.
  * Ebert  _et al._ (2017) M. A. Ebert, J. K. L. Michel, and F. J. Tackmann, JHEP 05, 088 (2017), arXiv:1702.00794 [hep-ph] .
  * Höschele  _et al._ (2014) M. Höschele, J. Hoff, A. Pak, M. Steinhauser, and T. Ueda, Comput. Phys. Commun. 185, 528 (2014), arXiv:1307.6925 [hep-ph] .
  * Korchemsky and Radyushkin (1987) G. P. Korchemsky and A. V. Radyushkin, Nucl. Phys. B283, 342 (1987).
  * Duhr and Dulat (2019) C. Duhr and F. Dulat, (2019), arXiv:1904.07279 [hep-th] .

## Supplemental material

###  V.1 Factorization theorem

Mode | Lab frame | Leptonic ($\hat{Y}=0$) frame
---|---|---
| $(+,-,\perp)$ | $(+,-,\perp)$
$p_{\bar{n}}$ | $(q^{+},\lambda^{2}q^{-},\lambda\sqrt{q^{-}q^{+}})$ | $Q\,(1,\lambda^{2},\lambda)$
$P_{\bar{n}}$ | $\displaystyle\Bigl{(}q^{+},\frac{\Lambda_{\mathrm{QCD}}^{2}}{q^{+}},\Lambda_{\mathrm{QCD}}\Bigr{)}$ | $Q\,(1,\lambda_{\mathrm{QCD}}^{2},\lambda_{\mathrm{QCD}})$
$P_{n}$ | $\displaystyle\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}^{2}}{q^{-}},q^{-},\Lambda_{\mathrm{QCD}}\Bigr{)}$ | $Q\,(\lambda_{\mathrm{QCD}}^{2},1,\lambda_{\mathrm{QCD}})$
$P_{s}$ | $\displaystyle\Bigl{(}\frac{1}{\lambda^{2}}\,\frac{\Lambda_{\mathrm{QCD}}^{2}}{q^{-}},\lambda^{2}q^{-},\Lambda_{\mathrm{QCD}}\Bigr{)}$ | $\displaystyle Q\,\Bigl{(}\frac{\lambda_{\mathrm{QCD}}^{2}}{\lambda^{2}},\lambda^{2},\lambda_{\mathrm{QCD}}\Bigr{)}$
$P_{us}$ | $\displaystyle\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}^{2}}{q^{-}},\frac{\Lambda_{\mathrm{QCD}}^{2}}{q^{+}},\frac{\Lambda_{\mathrm{QCD}}^{2}}{\sqrt{q^{+}q^{-}}}\Bigr{)}$ | $Q\,(\lambda_{\mathrm{QCD}}^{2},\lambda_{\mathrm{QCD}}^{2},\lambda_{\mathrm{QCD}}^{2})$
$P_{G}$ | $\displaystyle\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}^{2}}{q^{-}},\frac{\Lambda_{\mathrm{QCD}}^{2}}{q^{+}},\Lambda_{\mathrm{QCD}}\Bigr{)}$ | $Q\,(\lambda_{\mathrm{QCD}}^{2},\lambda_{\mathrm{QCD}}^{2},\lambda_{\mathrm{QCD}})$
Table S1: Relevant EFT modes in the limit
$\lambda_{\mathrm{QCD}}\sim\lambda^{2}\sim 1-x^{-}\ll 1$ in the lab (hadronic
center-of-mass) frame and the leptonic frame where $\hat{Y}=0$. In the right
column we used that in the leptonic frame,
$q^{\pm}\to\hat{q}^{\pm}=\sqrt{q^{+}q^{-}}\sim Q$.

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,

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}x_{-}\mathrm{d}x_{+}\,\mathrm{d}\vec{\kern
0.0ptq}_{T}}$
$\displaystyle=H_{ij}(q^{+}q^{-})\int\\!\mathrm{d}t\,f_{i}^{\mathrm{thr}}\Bigl{[}x_{-}\Bigl{(}1+\frac{t}{q^{+}q^{-}}\Bigr{)}\Bigr{]}B_{j}(t,\vec{\kern
0.0ptq}_{T},x_{+})\,.$ (S1)

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

$\lambda^{2}\sim 1-x_{-}\ll 1\qquad\text{for generic
}x_{+}\,,\qquad\lambda_{\mathrm{QCD}}\equiv\Lambda_{\mathrm{QCD}}/Q\sim\lambda^{2}\ll\lambda\,.$
(S2)

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

$\delta(\omega_{\bar{n}}-q^{+})\,\delta(\omega_{n}-q^{-})\,\delta(\vec{\kern
0.0ptk}_{{\bar{n}}\perp}+\vec{\kern
0.0ptq}_{T})\,\delta(k_{s}^{-}-k_{\bar{n}}^{-})\,.$ (S3)

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

$B_{j}(\omega_{\bar{n}}k^{-},\vec{\kern
0.0ptk}_{T},x_{+})=\int\\!\frac{\mathrm{d}k^{+}}{\omega_{\bar{n}}}\,\mathcal{S}(k^{-},k^{+},\vec{\kern
0.0ptk}_{T})\,f_{j}^{\mathrm{thr}}\Bigl{[}x_{+}\Bigl{(}1+\frac{k^{+}}{\omega_{\bar{n}}}\Bigr{)}\Bigr{]}\,\Bigl{[}1+\mathcal{O}(\lambda_{+})\Bigr{]}\,,$
(S4)

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

Figure S1: The $(z_{a},z_{b})$ plane as parametrized by $(z,y)$. The gray area
shows the integration region $x_{a}\leq z_{a}\leq 1$ and $x_{b}\leq z_{b}\leq
1$ (for the case of $x_{a}>x_{b}$) used to derive the distribution identities.
The dotted black lines indicate integration paths over $y$ at representative
fixed values of $z$. The solid red lines indicate the edge cases in $z$. The
dashed blue lines are lines of constant $y$.

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

$\displaystyle z$ $\displaystyle=z_{a}z_{b}\,,\qquad$ $\displaystyle y$
$\displaystyle=\frac{z_{b}(1-z_{a}^{2})}{(1-z_{a}z_{b})(z_{a}+z_{b})}\,,\qquad$
$\displaystyle 1-y$
$\displaystyle=\frac{z_{a}(1-z_{b}^{2})}{(1-z_{a}z_{b})(z_{a}+z_{b})}=y|_{a\leftrightarrow
b}\,,$ $\displaystyle z_{a}$
$\displaystyle=\sqrt{\frac{z[1-y(1-z)]}{z+y(1-z)}}\,,\qquad$ $\displaystyle
z_{b}$ $\displaystyle=z_{a}|_{y\leftrightarrow 1-y}\,,\qquad$
$\displaystyle\frac{\mathrm{d}z\,\mathrm{d}y}{\mathrm{d}z_{a}\mathrm{d}z_{b}}$
$\displaystyle=\frac{2[1-y(1-z)][1-(1-y)(1-z)]}{1-z^{2}}\,,$ (S5)

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

$\displaystyle\theta\bigl{[}z_{a}(z,y)\geq
x_{a}\bigr{]}\,\theta\bigl{[}z_{b}(z,y)\geq x_{b}\bigr{]}$
$\displaystyle=\theta\bigl{[}x_{a}x_{b}\leq
z<\min\\{x_{a},x_{b}\\}\bigr{]}\,\theta\bigl{[}y_{\mathrm{low}}(z)\leq y\leq
y_{\mathrm{high}}(z)\bigr{]}+\theta\bigl{[}\max\\{x_{a},x_{b}\\}\leq
z\bigr{]}$ $\displaystyle\quad+\theta\bigl{[}x_{b}\leq
z<x_{a}\bigr{]}\,\theta\bigl{[}y\leq
y_{\mathrm{high}}(z)\bigr{]}+\theta\bigl{(}x_{a}\leq
z<x_{b})\,\theta\bigl{[}y_{\mathrm{low}}(z)\leq y\bigr{]}\,,$ (S6)

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

$\displaystyle
y_{\mathrm{high}}(z)=\frac{z(1-x_{a}^{2})}{(1-z)(x_{a}^{2}+z)}\,,\qquad
y_{\mathrm{low}}(z)=\frac{x_{b}^{2}-z^{2}}{(1-z)(z+x_{b}^{2})}\,.$ (S7)

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

$\displaystyle\mathcal{L}_{n}(x)=\Bigl{[}\frac{\theta(x)\ln^{n}(x)}{x}\Bigr{]}_{+}=\lim_{\epsilon\to
0}\frac{\mathrm{d}}{\mathrm{d}x}\Bigl{[}\theta(x-\epsilon)\frac{\ln^{n+1}x}{n+1}\Bigr{]}\quad\text{with}\quad\mathcal{L}_{n}(x>0)=\frac{\ln^{n}(x)}{x}\,,\qquad\int^{1}\\!\mathrm{d}x\,\mathcal{L}_{n}(x)=0\,.$
(S8)

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.

$\mathrm{d}z\,\mathrm{d}y\,\text{LHS}$ | $=\mathrm{d}z_{a}\mathrm{d}z_{b}\,\text{RHS}$
---|---
$f(z)\,\delta(y)$ | $f(z_{b})\,\delta(1-z_{a})$
$r(z)\,\mathcal{L}_{0}(y)$ | $\displaystyle r(z_{b})\Bigl{[}\mathcal{L}_{0}(1-z_{a})+\delta(1-z_{a})\,\ln\frac{2z_{b}}{(1+z_{b})(1-z_{b})}\Bigr{]}+\mathcal{O}(1)$
$f(z)\,\delta(1-y)$ | $f(z_{a})\delta(1-z_{b})$
$r(z)\,\mathcal{L}_{0}(1-y)$ | $\mathcal{O}(1)$
$r(z)\,\mathcal{O}[y^{0}(1-y)^{0}]$ | $\mathcal{O}(1)$
$\mathcal{L}_{0}(1-z)\,\bigl{[}\mathcal{L}_{0}(y)+\mathcal{L}_{0}(1-y)\bigr{]}$ | $\displaystyle-\mathcal{L}_{1}(1-z_{a})\,\delta(1-z_{b})+\mathcal{L}_{0}(1-z_{a})\,\mathcal{L}_{0}(1-z_{b})-\delta(1-z_{a})\,\mathcal{L}_{1}(1-z_{b})$
|
$\displaystyle+\frac{\pi^{2}}{6}\delta(1-z_{a})\,\delta(1-z_{b})+\delta(1-z_{a})\,\frac{1}{1-z_{b}}\ln\frac{2z_{b}}{1+z_{b}}+\mathcal{O}(1)$
Table S2: Translation identities of two-dimensional plus distributions between
the $(z,y)$ and $(z_{a},z_{b})$ parametrizations. Here, $f(x)$ is an arbitrary
function of $x$, potentially distribution-valued for $x\to 1$, while
$r(x)=\mathcal{O}[(1-x)^{0}]$ has at most an integrable singularity for $x\to
1$. When indicated, the relations receive power corrections in $1-z_{a}$
starting at $\mathcal{O}(1)\equiv\mathcal{O}[(1-z_{a})^{0}]$. Overall factors
of $\theta(1-z)\,\theta(y)\,\theta(1-y)=\theta(1-z_{a})\,\theta(1-z_{b})$ are
understood on both sides.

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

$\displaystyle\int\\!\mathrm{d}z\mathrm{d}y\,\theta\bigl{[}z_{a}(z,y)-x_{a}\bigr{]}\,\theta\bigl{[}z_{b}(z,y)-x_{b}\bigr{]}\,\mathcal{L}_{0}(1-z)\Bigl{[}\mathcal{L}_{0}(y)+\mathcal{L}_{0}(1-y)\Bigr{]}$
$\displaystyle\qquad=F_{a}(x_{a})+F_{b}(x_{b})-F_{a}(x_{a}x_{b})-F_{b}(x_{a}x_{b})\,,$
$\displaystyle\text{with}\qquad F_{a,b}(z)$
$\displaystyle=-\ln(1-z)\ln(-z+\mathrm{i}0)-\text{Li}_{2}\Big{(}\frac{1-z}{1+x_{a,b}}\Big{)}-\text{Li}_{2}\Big{(}\frac{1-z}{1-x_{a,b}}-\mathrm{i}0\Big{)}-\text{Li}_{2}(z)\,,$
(S9)

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

$\displaystyle\mathrm{d}z\,\mathrm{d}y\,\mathcal{L}_{0}(1-z)\,\bigl{[}\mathcal{L}_{0}(y)+\mathcal{L}_{0}(1-y)\bigr{]}$
$\displaystyle\qquad=\mathrm{d}z_{a}\mathrm{d}z_{b}\biggl{[}\frac{\pi^{2}}{6}\,\delta(1-z_{a})\,\delta(1-z_{b})-\mathcal{L}_{1}(1-z_{a})\,\delta(1-z_{b})+\mathcal{L}_{0}(1-z_{a})\,\mathcal{L}_{0}(1-z_{b})-\delta(1-z_{a})\,\mathcal{L}_{1}(1-z_{b})$
$\displaystyle\qquad\qquad\qquad\quad+\delta(1-z_{a})\,\frac{1}{1-z_{b}}\ln\frac{2z_{b}}{1+z_{b}}+\delta(1-z_{b})\,\frac{1}{1-z_{a}}\ln\frac{2z_{a}}{1+z_{a}}+\frac{1}{(1+z_{a})(1+z_{b})}\biggr{]}\,.$
(S10)

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

$\displaystyle\Bigl{[}\frac{1+z^{2}}{1-z}\Bigr{]}_{+}=2\mathcal{L}_{0}(1-z)+\frac{3}{2}\delta(1-z)-(1+z)\,.$
(S11)

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

$\int\\!\mathrm{d}z_{a}\mathrm{d}z_{b}\,\delta(z-z_{a}z_{b})\,\mathcal{L}_{0}(1-z_{a})\mathcal{L}_{0}(1-z_{b})=2\mathcal{L}_{1}(1-z)-\frac{\pi^{2}}{6}\delta(1-z)-\frac{\ln
z}{1-z}\,.$ (S12)

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

$\displaystyle\hat{\sigma}_{ij}(z_{a},z_{b})$
$\displaystyle=H_{ik}\hat{\mathcal{I}}_{kj}(z_{a},z_{b})+H_{kj}\hat{\mathcal{I}}_{ki}(z_{b},z_{a})-H_{ij}\hat{S}(z_{a},z_{b})$
$\displaystyle\quad+\biggl{[}\frac{\mathrm{d}z\,\mathrm{d}y}{\mathrm{d}z_{a}\mathrm{d}z_{b}}\,\hat{\sigma}_{ij}\bigl{[}z(z_{a},z_{b}),y(z_{a},z_{b})\bigr{]}-H_{ik}\hat{\mathcal{I}}_{kj}(z_{a},z_{b})-H_{kj}\hat{\mathcal{I}}_{ki}(z_{b},z_{a})+H_{ij}\hat{S}(z_{a},z_{b})\biggr{]}_{z_{a,b}<1}\,.$
(S13)

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

$\displaystyle\sigma_{B,q}^{\mathrm{DY}}$
$\displaystyle=\frac{4\pi\alpha_{\text{em}}^{2}}{3N_{c}Q^{2}}\biggl{[}Q_{q}^{2}+\frac{(v_{q}^{2}+a_{q}^{2})(v_{\ell}^{2}+a_{\ell}^{2})-2Q_{q}v_{q}v_{\ell}(1-m_{Z}^{2}/Q^{2})}{(1-m_{Z}^{2}/Q^{2})^{2}+m_{Z}^{2}\Gamma_{Z}^{2}/Q^{4}}\biggr{]}\,,$
(S14)

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

$\displaystyle\frac{\mathrm{d}\sigma^{\mathrm{DY}}_{\text{LO}}}{\mathrm{d}Q\,\mathrm{d}Y}=\frac{2Q}{E_{\mathrm{cm}}^{2}}\,\frac{\mathrm{d}\sigma^{\mathrm{DY}}_{\text{LO}}}{\mathrm{d}x_{a}\mathrm{d}x_{b}}=\frac{2Q}{E_{\mathrm{cm}}^{2}}\,\sum_{q}\sigma_{B,q}^{\mathrm{DY}}\,\bigl{[}f_{q}(x_{a})\,f_{\bar{q}}(x_{b})+f_{\bar{q}}(x_{a})\,f_{q}(x_{b})\bigr{]}\,,$
(S15)

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

$\displaystyle\hat{\sigma}_{ij}(z_{a},z_{b},Q,\mu)=\sum_{n=0}^{\infty}\Bigl{[}\frac{\alpha_{s}(\mu)}{4\pi}\Bigr{]}^{n}\,\sigma_{ij}^{(n)}(z_{a},z_{b},Q,\mu)\,.$
(S16)

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

$\displaystyle\frac{1}{\sigma_{B,q}^{\mathrm{DY}}}\frac{\hat{\sigma}^{(1)}_{q{\bar{q}}}(z_{a},z_{b},Q,\mu)}{C_{F}}$
$\displaystyle=\delta({\bar{z}}_{a})\delta({\bar{z}}_{b})(2\pi^{2}-16)+4\mathcal{L}_{1}({\bar{z}}_{a})\delta({\bar{z}}_{b})+4\mathcal{L}_{0}({\bar{z}}_{a})\mathcal{L}_{0}({\bar{z}}_{b})+4\delta({\bar{z}}_{a})\mathcal{L}_{1}({\bar{z}}_{b})$
$\displaystyle\quad+\biggl{\\{}-2(1+z_{b})\mathcal{L}_{0}({\bar{z}}_{a})+\delta({\bar{z}}_{a})\Bigl{[}2{\bar{z}}_{b}-4(1+z_{b})\ln{\bar{z}}_{b}-\frac{2(1+z_{b}^{2})\ln
z_{b}}{{\bar{z}}_{b}}$
$\displaystyle\qquad\quad+\frac{4}{{\bar{z}}_{b}}\ln\frac{2z_{b}}{1+z_{b}}+2(1+z_{b})\ln\frac{1-z_{b}^{2}}{2z_{b}}\Bigr{]}-4\ln\frac{\mu}{Q}\,\delta({\bar{z}}_{a})\Bigl{[}2\mathcal{L}_{0}({\bar{z}}_{b})+\frac{3}{2}\delta({\bar{z}}_{b})-(1+z_{b})\Bigr{]}$
$\displaystyle\qquad\quad+\frac{2(z_{a}^{2}+z_{b}^{2})[(1+z_{a})^{2}+z_{a}z_{b}(3+2z_{a}+z_{a}z_{b})]}{(1+z_{a})(1+z_{b})(z_{a}+z_{b})^{2}}+(z_{a}\leftrightarrow
z_{b})\biggr{\\}}\,,$ (S17)

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

$\displaystyle\frac{1}{\sigma_{B,q}^{\mathrm{DY}}}\frac{\hat{\sigma}_{qg}^{(1)}(z_{a},z_{b},Q,\mu)}{T_{F}}$
$\displaystyle=2(z_{b}^{2}+{\bar{z}}_{b}^{2})\mathcal{L}_{0}({\bar{z}}_{a})+\delta({\bar{z}}_{a})\Bigl{[}2(z_{b}^{2}+{\bar{z}}_{b}^{2})\ln\frac{2{\bar{z}}_{b}}{1+z_{b}}+4z_{b}{\bar{z}}_{b}\Bigr{]}-4\ln\frac{\mu}{Q}\,\delta({\bar{z}}_{a})\bigl{(}z_{b}^{2}+{\bar{z}}_{b}^{2}\bigr{)}$
(S18)
$\displaystyle\quad+\frac{1}{(1+z_{a})(z_{a}+z_{b})^{3}}\Bigl{[}-4z_{a}^{5}z_{b}^{3}-4z_{a}^{4}z_{b}^{2}(-1+z_{b}+2z_{b}^{2})+2z_{a}^{3}(1+4z_{b}^{2}+2z_{b}^{3}-4z_{b}^{4}-4z_{b}^{5})$
$\displaystyle\qquad+2z_{a}^{2}z_{b}(1+4z_{b}+8z_{b}^{2}-8z_{b}^{3}-4z_{b}^{4})+2z_{a}z_{b}^{2}(1+4z_{b}-2z_{b}^{2}-4z_{b}^{3})-2z_{b}^{3}(1-2z_{b}+2z_{b}^{2})\Bigr{]}\,.$

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}$,

$\mathcal{L}_{\mathrm{eff}}(m_{H})=-\frac{C_{t}}{12\pi
v}\,\alpha_{s}\,G_{\mu\nu}^{a}G^{a,\mu\nu}H\,,\qquad
C_{t}=1+\frac{\alpha_{s}}{4\pi}(5C_{A}-3C_{F})+\mathcal{O}(\alpha_{s}^{2})\,.$
(S19)

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

$\displaystyle\sigma_{B}^{ggH}=\frac{1}{72\pi
v^{2}(N_{c}^{2}-1)}\,,\qquad\frac{\mathrm{d}\sigma^{ggH}_{\text{LO}}}{\mathrm{d}Y}=x_{a}x_{b}\,\sigma_{B}^{ggH}\,\alpha_{s}^{2}\,f_{g}(x_{a})\,f_{g}(x_{b})\,.$
(S20)

We write the partonic cross section as

$\displaystyle\hat{\sigma}_{ij}(z_{a},z_{b},m_{t},m_{H},\mu)\equiv\sigma_{B}^{ggH}\,\alpha_{s}^{2}(\mu)\,|C_{t}(m_{t},\mu)|^{2}\,\sum_{n=0}^{\infty}\Bigl{[}\frac{\alpha_{s}(\mu)}{4\pi}\Bigr{]}^{n}\,\hat{\eta}_{ij}^{(n)}(z_{a},z_{b},m_{H},\mu)\,.$
(S21)

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

$\displaystyle\frac{\hat{\eta}_{gg}^{(1)}(z_{a},z_{b},m_{H},\mu)}{C_{A}}$
$\displaystyle=2\pi^{2}\delta({\bar{z}}_{a})\delta({\bar{z}}_{b})+4\mathcal{L}_{1}({\bar{z}}_{a})\delta({\bar{z}}_{b})+4\delta({\bar{z}}_{a})\mathcal{L}_{1}({\bar{z}}_{b})+4\mathcal{L}_{0}({\bar{z}}_{a})\mathcal{L}_{0}({\bar{z}}_{b})$
(S22)
$\displaystyle\quad+\biggl{\\{}4\mathcal{L}_{0}({\bar{z}}_{a})\Bigl{[}\frac{1}{z_{b}}-2+z_{b}-z_{b}^{2}\Bigr{]}+4\delta({\bar{z}}_{a})\frac{1}{{\bar{z}}_{b}}\Bigl{[}\frac{\ln(2{\bar{z}}_{b})}{z_{b}}-3\ln{\bar{z}}_{b}+2\ln\frac{1+z_{b}}{2}-\frac{\ln(1+z_{b})}{z_{b}}$
$\displaystyle\qquad\quad-
z_{b}(3-2z_{b}+z_{b}^{2})\ln\frac{1+z_{b}}{2{\bar{z}}_{b}}\Bigr{]}-4\ln\frac{\mu}{m_{H}}\delta({\bar{z}}_{a})\Bigl{[}2\mathcal{L}_{0}({\bar{z}}_{b})\frac{(1-z_{b}+z_{b}^{2})^{2}}{z_{b}}\Bigr{]}$
$\displaystyle\qquad\quad+\frac{4z_{b}}{z_{a}(1+z_{a})(1+z_{b})(z_{a}+z_{b})^{4}}\Bigl{[}2z_{b}^{2}+z_{b}^{3}+3z_{a}^{6}z_{b}^{4}+2z_{a}^{5}z_{b}^{3}(5+5z_{b}+2z_{b}^{2})$
$\displaystyle\qquad\qquad+z_{a}^{4}z_{b}^{2}(16+17z_{b}+12z_{b}^{2}+6z_{b}^{3}+2z_{b}^{4})+z_{a}^{3}z_{b}(5+22z_{b}+12z_{b}^{2}+8z_{b}^{3}+8z_{b}^{4}+2z_{b}^{5})$
$\displaystyle\qquad\qquad+z_{a}^{2}(3+2z_{b}^{2}+7z_{b}^{3}+2z_{b}^{4}+4z_{b}^{5}+2z_{b}^{6})+z_{a}z_{b}(4+z_{b}+z_{b}^{2}+z_{b}^{3}+z_{b}^{5})\Bigr{]}+(z_{a}\leftrightarrow
z_{b})\biggr{\\}}\,.$

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)$,

$\displaystyle\frac{\hat{\eta}_{gq}^{(1)}(z_{a},z_{b},m_{H},\mu)}{C_{F}}$
$\displaystyle=2\mathcal{L}_{0}({\bar{z}}_{a})\,\frac{2-2z_{b}+z_{b}^{2}}{z_{b}}+2\delta({\bar{z}}_{a})\Bigl{[}z_{b}+\frac{2-2z_{b}+z_{b}^{2}}{z_{b}}\Bigl{(}\ln\frac{2{\bar{z}}_{b}}{1+z_{b}}-2\ln\frac{\mu}{m_{H}}\Bigr{)}\Bigr{]}$
$\displaystyle\quad+\frac{2}{(1+z_{a})z_{b}(z_{a}+z_{b})^{3}}\Bigl{[}z_{a}^{3}(2-2z_{b}+z_{b}^{2})+z_{a}^{2}(4-2z_{b}-4z_{b}^{2}+7z_{b}^{3}-2z_{b}^{5})$
$\displaystyle\qquad+z_{a}z_{b}(4-4z_{b}+4z_{b}^{2}+z_{b}^{3}-2z_{b}^{4})-z_{b}^{2}(-2+2z_{b}-2z_{b}^{2}+z_{b}^{3})\Bigr{]}\,.$
(S23)

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

$\displaystyle\frac{\hat{\eta}_{q\bar{q}}^{(1)}(z_{a},z_{b},m_{H},\mu)}{C_{F}}=\frac{N_{c}^{2}-1}{N_{c}}\,\frac{4(1+z_{a}z_{b})(z_{a}^{4}z_{b}^{2}+z_{a}^{2}z_{b}^{4}-4z_{a}^{2}z_{b}^{2}+z_{a}^{2}+z_{b}^{2})}{(z_{a}+z_{b})^{4}}\,.$
(S24)

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

$\mathrm{d}z\,\mathrm{d}y\,\frac{1}{1-z}\Bigl{(}\frac{1}{y}+\frac{1}{1-y}\Bigr{)}=\mathrm{d}z_{a}\mathrm{d}z_{b}\,\frac{1}{(1-z_{a})(1-z_{b})}\bigl{[}1+\mathcal{O}(1-z_{a},1-z_{b})\bigr{]}\,.$
(S25)

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

$\frac{x_{a}}{z_{a}(z,y)}=x_{a}\Bigl{\\{}1+y(1-z)+\mathcal{O}[(1-z)^{2}]\Bigr{\\}}\,,\qquad\frac{x_{b}}{z_{b}(z,y)}=x_{b}\Bigl{\\{}1+(1-y)(1-z)+\mathcal{O}[(1-z)^{2}]\Bigr{\\}}\,.$
(S26)

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})$,

$y(1-z)=1-z_{a}+\mathcal{O}[(1-z)^{2}]\,,\qquad(1-y)(1-z)=1-z_{b}+\mathcal{O}[(1-z)^{2}]\,,$
(S27)

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$,

$f(x)\equiv\theta(1-x)(1-x)^{\alpha}\,.$ (S28)

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,

$\displaystyle
A(z,y)\equiv\mathcal{L}_{0}(1-z)\,\frac{\delta(y)+\delta(1-y)}{2}\,,\qquad
B(z,y)\equiv\mathcal{L}_{0}(1-z)\,\delta\Bigl{(}y-\frac{1}{2}\Bigr{)}\,,$
(S29)

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

$\displaystyle\int\\!\mathrm{d}z\,\mathrm{d}y\,A(z,y)\,f\Bigl{[}\frac{x_{a}}{z_{a}(z,y)}\Bigr{]}f\Bigl{[}\frac{x_{b}}{z_{b}(z,y)}\Bigr{]}$
$\displaystyle=f(x_{a})f(x_{b})\Bigl{[}\frac{1}{2}\ln\bigl{(}1-x_{a}\bigr{)}+\frac{1}{2}\ln\bigl{(}1-x_{b}\bigr{)}-H_{\alpha}+\mathcal{O}(1\\!-\\!x_{a},1\\!-\\!x_{b})\Bigr{]}\,,$
$\displaystyle\int\\!\mathrm{d}z\,\mathrm{d}y\,B(z,y)\,f\Bigl{[}\frac{x_{a}}{z_{a}(z,y)}\Bigr{]}f\Bigl{[}\frac{x_{b}}{z_{b}(z,y)}\Bigr{]}$
$\displaystyle=f(x_{a})f(x_{b})\Bigl{[}\ln\bigl{(}1-\max\\{x_{a},x_{b}\\}\bigr{)}-H_{\alpha}+\mathcal{O}(1\\!-\\!x_{a},1\\!-\\!x_{b})\Bigr{]}\,,$
(S30)

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

$u\equiv\frac{1}{2}\ln\frac{z_{a}}{z_{b}}\,,\qquad-u_{\mathrm{max}}\leq u\leq
u_{\mathrm{max}}\,,\qquad u_{\mathrm{max}}\equiv\ln\frac{1}{\sqrt{z}}\,.$
(S31)

(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$,

$\tilde{C}(z,M)\equiv\int\\!\mathrm{d}u\,e^{\mathrm{i}Mu}\,C(z,u)\stackrel{{\scriptstyle?}}{{=}}\int\\!\mathrm{d}u\,C(z,u)\bigl{[}1+\mathcal{O}(1-z)\bigr{]}\,.$
(S32)

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

$C(z,u)=\int\\!\frac{\mathrm{d}M}{2\pi}\,e^{-\mathrm{i}Mu}\,\tilde{C}(z,M)\stackrel{{\scriptstyle?}}{{=}}\delta(u)\int\\!\mathrm{d}u^{\prime}\,C(z,u^{\prime})\bigl{[}1+\mathcal{O}(1-z)\bigr{]}\,.$
(S33)

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

$\mathrm{d}y\,A(z,y)=\mathrm{d}u\,\mathcal{L}_{0}(1-z)\,\frac{\delta(u+u_{\mathrm{max}})+\delta(u-u_{\mathrm{max}})}{2}\,,\qquad\mathrm{d}y\,B(z,y)=\mathrm{d}u\,\mathcal{L}_{0}(1-z)\,\delta(u)\,.$
(S34)

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$,

$A(z,y)=\mathcal{L}_{0}(1-z)\,\frac{\delta(y)+\delta(1-y)}{2}\stackrel{{\scriptstyle?}}{{=}}\mathcal{L}_{0}(1-z)\,\delta\Bigl{(}y-\frac{1}{2}\Bigr{)}=B(z,y)\,.$
(S35)

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

$\displaystyle H_{ij}^{\mathrm{DY}}(Q^{2},\mu)$
$\displaystyle=\sum_{q}\sigma_{B,q}^{\mathrm{DY}}\,\bigl{(}\delta_{iq}\delta_{j{\bar{q}}}+\delta_{i{\bar{q}}}\delta_{jq}\bigr{)}\,\lvert
C_{q{\bar{q}}}^{V}(Q^{2},\mu)\rvert^{2}\,,$ (S36)

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

$\displaystyle H_{ij}^{ggH}(m_{t},m_{H},\mu)$
$\displaystyle=\sigma_{B}^{ggH}\,\bigl{|}\alpha_{s}(\mu)\,C_{t}(m_{t},\mu)\bigr{|}^{2}\,\delta_{ig}\delta_{jg}\,\lvert
C_{gg}(m_{H},\mu)\rvert^{2}\,,$ (S37)

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,

$H_{ij}(Q^{2},\mu)=\sum_{n=0}^{\infty}\Bigl{[}\frac{\alpha_{s}(\mu)}{4\pi}\Bigr{]}^{n}H_{ij}^{(n)}(Q^{2},\mu)\,.$
(S38)

####  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

$\displaystyle\tilde{B}_{i}(t,x,\mu)=\int\\!\frac{\mathrm{d}z}{z}\,\tilde{\mathcal{I}}_{ij}(t,z,\mu)\,f_{j}\Bigl{(}\frac{x}{z},\mu\Bigr{)}\,\Bigl{[}1+\mathcal{O}\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}^{2}}{t}\Bigr{)}\Bigr{]}\,,$
(S39)

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,

$\displaystyle\mu\frac{\mathrm{d}}{\mathrm{d}\mu}\tilde{B}_{i}(t,x,\mu)=\int\\!\mathrm{d}t^{\prime}\,\gamma_{B}^{i}(t-t^{\prime},\mu)\,\tilde{B}_{i}(t^{\prime},x,\mu)\,,\quad\gamma_{B}^{i}(t,\mu)=-2\Gamma_{\mathrm{cusp}}^{i}\bigl{[}\alpha_{s}(\mu)\bigr{]}\,\mathcal{L}_{0}(t,\mu^{2})+\gamma_{B}^{i}\bigl{[}\alpha_{s}(\mu)\bigr{]}\delta(t)\,,$
(S40)

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)

$\displaystyle\mu\frac{\mathrm{d}}{\mathrm{d}\mu}\tilde{\mathcal{I}}_{ij}(t,z,\mu)$
$\displaystyle=\int\\!\mathrm{d}t^{\prime}\,\int\\!\frac{\mathrm{d}z^{\prime}}{z^{\prime}}\,\tilde{\mathcal{I}}_{ik}\Bigl{(}t-t^{\prime},\frac{z}{z^{\prime}},\mu\Bigr{)}\,\Bigl{\\{}\delta_{kj}\delta(1-z^{\prime})\,\gamma_{B}^{i}(t^{\prime},\mu)-2\delta(t^{\prime})\,P_{kj}\bigl{[}\alpha_{s}(\mu),z^{\prime}\bigr{]}\Bigr{\\}}\,,$
(S41)

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

$\displaystyle\tilde{\mathcal{I}}_{ij}(t,z,\mu)=\sum_{n=0}^{\infty}\Bigl{[}\frac{\alpha_{s}(\mu)}{4\pi}\Bigr{]}^{n}\tilde{\mathcal{I}}_{ij}^{(n)}(t,z,\mu)\,.$
(S42)

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

$\displaystyle\tilde{\mathcal{I}}_{ij}^{(0)}(t,z,\mu)$
$\displaystyle=\delta(t)\,\delta_{ij}\delta(1-z)\,,$
$\displaystyle\tilde{\mathcal{I}}_{ij}^{(1)}(t,z,\mu)$
$\displaystyle=\mathcal{L}_{1}(t,\mu^{2})\,\Gamma_{0}^{i}\,\delta_{ij}\delta(1-z)+\mathcal{L}_{0}(t,\mu^{2})\Bigl{[}P^{(0)}_{ij}(z)-\frac{\gamma_{B\,0}^{i}}{2}\,\delta_{ij}\delta(1-z)\Bigr{]}+\delta(t)\,\tilde{I}^{(1)}_{ij}(z)\,,$
$\displaystyle\tilde{\mathcal{I}}_{ij}^{(2)}(t,z,\mu)$
$\displaystyle=\mathcal{L}_{3}(t,\mu^{2})\,\frac{(\Gamma_{0}^{i})^{2}}{2}\,\delta_{ij}\delta(1-z)$
$\displaystyle\quad+\mathcal{L}_{2}(t,\mu^{2})\,\frac{\Gamma_{0}^{i}}{2}\biggl{[}-\Bigl{(}\beta_{0}+\frac{3}{2}\gamma_{B\,0}^{i}\Bigr{)}\delta_{ij}\delta(1-z)+3P^{(0)}_{ij}(z)\biggr{]}$
$\displaystyle\quad+\mathcal{L}_{1}(t,\mu^{2})\biggl{[}\Bigl{(}-\frac{\pi^{2}}{6}(\Gamma_{0}^{i})^{2}+\frac{\beta_{0}}{2}\gamma_{B\,0}^{i}+\frac{(\gamma_{B\,0}^{i})^{2}}{4}+\Gamma_{1}^{i}\Bigr{)}\delta_{ij}\delta(1-z)-\bigl{(}\beta_{0}+\gamma_{B\,0}^{i}\bigr{)}P^{(0)}_{ij}(z)$
$\displaystyle\qquad\qquad\qquad+\bigl{(}P^{(0)}_{ik}\otimes
P^{(0)}_{kj}\bigr{)}(z)+\Gamma_{0}^{i}\tilde{I}^{(1)}_{ij}(z)\biggr{]}$
$\displaystyle\quad+\mathcal{L}_{0}(t,\mu^{2})\biggl{[}\Bigl{(}\zeta_{3}(\Gamma_{0}^{i})^{2}+\frac{\pi^{2}}{12}\Gamma_{0}^{i}\gamma_{B\,0}^{i}-\frac{\gamma_{B\,1}^{i}}{2}\Bigr{)}\delta_{ij}\delta(1-z)-\frac{\pi^{2}}{6}\Gamma_{0}^{i}P^{(0)}_{ij}(z)+P^{(1)}_{ij}(z)$
$\displaystyle\qquad\qquad\qquad-\Bigl{(}\beta_{0}+\frac{\gamma_{B\,0}^{i}}{2}\Bigr{)}\tilde{I}^{(1)}_{ij}(z)+\bigl{(}\tilde{I}^{(1)}_{ik}\otimes
P^{(0)}_{kj}\bigr{)}(z)\biggr{]}$
$\displaystyle\quad+\delta(t)\,\tilde{I}^{(2)}_{ij}(z)\,,$
$\displaystyle\tilde{\mathcal{I}}_{ij}^{(3)}(t,z,\mu)$
$\displaystyle=\mathcal{L}_{5}(t,\mu^{2})\,\frac{(\Gamma_{0}^{i})^{3}}{8}\delta_{ij}\delta(1-z)$
$\displaystyle\quad+\mathcal{L}_{4}(t,\mu^{2})\,\frac{5}{8}(\Gamma_{0}^{i})^{2}\biggl{[}-\Bigl{(}\frac{2}{3}\beta_{0}+\frac{\gamma_{B\,0}^{i}}{2}\Bigr{)}\delta_{ij}\delta(1-z)+P^{(0)}_{ij}(z)\biggr{]}$
$\displaystyle\quad+\mathcal{L}_{3}(t,\mu^{2})\,\Gamma_{0}^{i}\biggl{[}\Bigl{(}-\frac{\pi^{2}}{6}(\Gamma_{0}^{i})^{2}+\frac{\beta_{0}^{2}}{3}+\frac{5}{6}\beta_{0}\gamma_{B\,0}^{i}+\frac{(\gamma_{B\,0}^{i})^{2}}{4}+\Gamma_{1}^{i}\Bigr{)}\delta_{ij}\delta(1-z)-\Bigl{(}\frac{5}{3}\beta_{0}+\gamma_{B\,0}^{i}\Bigr{)}P^{(0)}_{ij}(z)$
$\displaystyle\qquad\qquad\qquad\quad+\bigl{(}P^{(0)}_{ik}\otimes
P^{(0)}_{kj}\bigr{)}(z)+\frac{\Gamma_{0}^{i}}{2}\,\tilde{I}^{(1)}_{ij}(z)\biggr{]}$
$\displaystyle\quad+\dotsb+\delta(t)\,\tilde{I}^{(3)}_{ij}(z)\,,$ (S43)

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

$\displaystyle\beta(\alpha_{s})=-2\alpha_{s}\sum_{n=0}^{\infty}\beta_{n}\left(\frac{\alpha_{s}}{4\pi}\right)^{n+1}\,,\qquad\beta_{0}=\frac{11}{3}\,C_{A}-\frac{4}{3}\,T_{F}\,n_{f}\,,\qquad\beta_{1}=\frac{34}{3}\,C_{A}^{2}-2T_{F}\,n_{f}\Bigl{(}\frac{10}{3}\,C_{A}+2C_{F}\Bigr{)}\,.$
(S44)

The anomalous dimensions are expanded as

$\Gamma^{i}_{\mathrm{cusp}}(\alpha_{s})=\sum_{n=0}^{\infty}\Gamma^{i}_{n}\Bigl{(}\frac{\alpha_{s}}{4\pi}\Bigr{)}^{n+1}\,,\qquad\gamma_{B}^{i}(\alpha_{s})=\sum_{n=0}^{\infty}\gamma^{i}_{B\,n}\Bigl{(}\frac{\alpha_{s}}{4\pi}\Bigr{)}^{n+1}\,,\qquad
P_{ij}(\alpha_{s},z)=\sum_{n=0}^{\infty}P^{(n)}_{ij}(z)\left(\frac{\alpha_{s}}{4\pi}\right)^{n+1}\,.$
(S45)

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

$\displaystyle\Gamma_{0}^{i}=4C_{i}\,,\qquad\Gamma_{1}^{i}=4C_{i}\Bigl{[}C_{A}\Bigl{(}\frac{67}{9}-\frac{\pi^{2}}{3}\Bigr{)}-\frac{20}{9}\,T_{F}\,n_{f}\Bigr{]}=\frac{4}{3}C_{i}\bigl{[}(4-\pi^{2})C_{A}+5\beta_{0}\bigr{]}\,,$

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)

$\displaystyle\gamma_{B\,0}^{q}$ $\displaystyle=6C_{F}\,,\qquad$
$\displaystyle\gamma_{B\,1}^{q}$
$\displaystyle=C_{F}\Bigl{[}\Bigl{(}\frac{146}{9}-80\zeta_{3}\Bigr{)}C_{A}+(3-4\pi^{2}+48\zeta_{3})C_{F}+\Bigl{(}\frac{121}{9}+\frac{2\pi^{2}}{3}\Bigr{)}\beta_{0}\Bigr{]}\,,$
$\displaystyle\gamma_{B\,0}^{g}$ $\displaystyle=2\beta_{0}\,,\qquad$
$\displaystyle\gamma_{B\,1}^{g}$
$\displaystyle=\Bigl{(}\frac{182}{9}-32\zeta_{3}\Bigr{)}C_{A}^{2}+\Bigl{(}\frac{94}{9}-\frac{2\pi^{2}}{3}\Bigr{)}C_{A}\,\beta_{0}+2\beta_{1}.$
(S46)

The one-loop PDF anomalous dimensions are

$\displaystyle P^{(0)}_{q_{i}q_{j}}(z)$
$\displaystyle=P^{(0)}_{\bar{q}_{i}\bar{q}_{j}}(z)=2C_{F}\,\delta_{ij}\,\theta(z)P_{qq}(z)\,,\quad$
$\displaystyle P^{(0)}_{q_{i}g}(z)$
$\displaystyle=P^{(0)}_{\bar{q}_{i}g}(z)=2T_{F}\,\theta(z)P_{qg}(z)\,,$
$\displaystyle P^{(0)}_{gg}(z)$
$\displaystyle=2C_{A}\,\theta(z)P_{gg}(z)+\beta_{0}\,\delta(1-z)\,,\quad$
$\displaystyle P^{(0)}_{gq_{i}}(z)$
$\displaystyle=P^{(0)}_{g\bar{q}_{i}}(z)=2C_{F}\,\theta(z)P_{gq}(z)\,,$ (S47)

with the standard color-stripped splitting functions,

$\displaystyle P_{qq}(z)$
$\displaystyle=\mathcal{L}_{0}(1-z)(1+z^{2})+\frac{3}{2}\delta(1-z)\,,\qquad$
$\displaystyle P_{qg}(z)$
$\displaystyle=\theta(1-z)\bigl{[}1-2z(1-z)\bigr{]}\,,$ $\displaystyle
P_{gg}(z)$
$\displaystyle=2\mathcal{L}_{0}(1-z)\frac{(1-z+z^{2})^{2}}{z}\,,\qquad$
$\displaystyle P_{gq}(z)$ $\displaystyle=\theta(1-z)\frac{1+(1-z)^{2}}{z}\,.$
(S48)

####  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

$\displaystyle\tilde{\mathcal{I}}_{ij}(t,z,\mu)=\int\\!\mathrm{d}^{2}\vec{\kern
0.0ptk}_{T}\,\mathcal{I}_{ij}\Bigl{(}t-\frac{k_{T}^{2}}{2},\vec{\kern
0.0ptk}_{T},z,\mu\Bigr{)}\,,$ (S49)

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,

$\displaystyle\tilde{\mathcal{I}}^{(1)}_{ij}(t,z,\mu)=\mathcal{I}_{ij}^{(1)}(t,z,\mu)+\delta(t)\,P^{(0)}_{ij}(z)\ln\frac{2z}{1+z}\,,$
(S50)

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

$\displaystyle\tilde{I}^{(1)}_{q_{i}q_{j}}=\tilde{I}^{(1)}_{\bar{q}_{i}\bar{q}_{j}}\equiv\delta_{ij}\tilde{I}^{(1)}_{qq}(z)$
$\displaystyle=2C_{F}\,\delta_{ij}\,\theta(z)\Bigl{[}\mathcal{L}_{1}(1-z)(1+z^{2})-\frac{\pi^{2}}{6}\delta(1-z)+\theta(1-z)(1-z)+P_{qq}(z)\ln\frac{2}{1+z}\Bigr{]}\,,$
$\displaystyle\tilde{I}^{(1)}_{q_{i}g}=\tilde{I}^{(1)}_{\bar{q}_{i}g}=\tilde{I}^{(1)}_{qg}(z)$
$\displaystyle=2T_{F}\,\theta(z)\Bigl{[}P_{qg}(z)\ln\frac{2(1-z)}{1+z}+\theta(1-z)\,2z(1-z)\Bigr{]}\,,$
$\displaystyle\tilde{I}_{gg}^{(1)}(z)$
$\displaystyle=2C_{A}\,\theta(z)\Bigl{[}\mathcal{L}_{1}(1-z)\frac{2(1-z+z^{2})^{2}}{z}-\frac{\pi^{2}}{6}\delta(1-z)+P_{gg}(z)\ln\frac{2}{1+z}\Bigr{]}\,,$
$\displaystyle\tilde{I}^{(1)}_{gq_{i}}=\tilde{I}^{(1)}_{g\bar{q}_{j}}=\tilde{I}_{gq}^{(1)}(z)$
$\displaystyle=2C_{F}\,\theta(z)\Bigl{[}P_{gq}(z)\ln\frac{2(1-z)}{1+z}+\theta(1-z)\,z\Bigr{]}\,.$
(S51)

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

$\displaystyle\tilde{I}_{q_{i}q_{j}}^{(2)}(z)=\tilde{I}_{{\bar{q}}_{i}{\bar{q}}_{j}}^{(2)}(z)$
$\displaystyle=\delta_{ij}\tilde{I}_{qqV}^{(2)}(z)+\tilde{I}_{qqS}^{(2)}(z)\,,\qquad\tilde{I}_{q_{i}{\bar{q}}_{j}}^{(2)}(z)=\tilde{I}_{{\bar{q}}_{i}q_{j}}^{(2)}(z)=\delta_{ij}\tilde{I}_{q{\bar{q}}V}^{(2)}(z)+\tilde{I}_{qqS}^{(2)}(z)\,,$
$\displaystyle\tilde{I}_{q_{i}g}^{(2)}(z)=\tilde{I}_{{\bar{q}}_{i}g}^{(2)}(z)$
$\displaystyle=\tilde{I}_{qg}^{(2)}(z)\,.$ (S52)

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:

$\displaystyle\tilde{I}^{(2)}_{qqV}(z)$
$\displaystyle=4C_{F}^{2}\Bigl{\\{}D_{qqV,C_{F}}(z)-\frac{2}{1-z}\tilde{T}_{3}(z)+\frac{1+z^{2}}{1-z}\bigl{[}\tilde{V}_{3}(z)-2\tilde{U}_{3}(z)\bigr{]}+\tilde{C}_{qqV,C_{F}}(z)\Bigr{\\}}$
$\displaystyle\quad+4C_{F}C_{A}\Bigl{\\{}D_{qqV,C_{A}}(z)+\frac{1+z^{2}}{1-z}\bigl{[}\tilde{U}_{3}(z)+\tilde{R}_{qqV}(z)\bigr{]}+\tilde{C}_{qqV,C_{A}}(z)\Bigr{\\}}+4C_{F}\beta_{0}\Bigl{[}D_{qqV,\beta_{0}}(z)+\tilde{C}_{qqV,\beta_{0}}(z)\Bigr{]}\,,$
$\displaystyle\tilde{I}^{(2)}_{q\bar{q}V}(z)$
$\displaystyle=4C_{F}(2C_{F}-C_{A})\Bigl{[}\frac{1+z^{2}}{1+z}\tilde{S}_{3}(z)+\tilde{C}_{q\bar{q}V}(z)\Bigr{]}\,,$
$\displaystyle\tilde{I}^{(2)}_{qqS}(z)$
$\displaystyle=4C_{F}T_{F}\Bigl{[}-2(1+z)\,\tilde{T}_{3}(z)+\tilde{C}_{qqS}(z)\Bigr{]}\,,$
$\displaystyle\tilde{I}^{(2)}_{qg}(z)$
$\displaystyle=4T_{F}C_{F}\Bigl{\\{}-2(1-z)^{2}\,\tilde{T}_{3}(z)+P_{qg}(z)\bigl{[}\tilde{V}_{3}(z)+\tilde{R}_{qg,C_{F}}(z)\bigr{]}+\tilde{C}_{qg,C_{F}}(z)\Bigr{\\}}$
$\displaystyle\quad+4T_{F}C_{A}\Bigl{\\{}-2(1+4z)\,\tilde{T}_{3}(z)-P_{qg}(z)\bigl{[}\tilde{U}_{3}(z)+\tilde{R}_{qg,C_{A}}(z)\bigr{]}+P_{qg}(-z)\,\tilde{S}_{3}(z)+\tilde{C}_{qg,C_{A}}(z)\Bigr{\\}}\,.$
(S53)

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),

$\displaystyle D_{qqV,C_{F}}(z)$
$\displaystyle=(1+z^{2})\Bigl{[}\mathcal{L}_{3}(1-z)-\frac{5\pi^{2}}{6}\mathcal{L}_{1}(1-z)+4\zeta_{3}\mathcal{L}_{0}(1-z)\Bigr{]}+\frac{7\pi^{4}}{120}\delta(1-z)\,,$
$\displaystyle D_{qqV,C_{A}}(z)$
$\displaystyle=(1+z^{2})\Bigl{[}\Bigl{(}\frac{2}{3}-\frac{\pi^{2}}{6}\Bigr{)}\mathcal{L}_{1}(1-z)+\Bigl{(}-\frac{8}{9}+\frac{7\zeta_{3}}{2}\Bigr{)}\mathcal{L}_{0}(1-z)\Bigr{]}+\Bigl{(}\frac{52}{27}-\frac{\pi^{2}}{6}-\frac{\pi^{4}}{36}\Bigr{)}\delta(1-z)\,,$
$\displaystyle D_{qqV,\beta_{0}}(z)$
$\displaystyle=(1+z^{2})\Bigl{[}-\frac{1}{4}\mathcal{L}_{2}(1-z)+\frac{5}{6}\mathcal{L}_{1}(1-z)+\Bigl{(}-\frac{7}{9}+\frac{\pi^{2}}{12}\Bigr{)}\mathcal{L}_{0}(1-z)\Bigr{]}+\Bigl{(}\frac{41}{27}-\frac{5\pi^{2}}{24}-\frac{5\zeta_{3}}{6}\Bigr{)}\delta(1-z)\,.$
(S54)

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

$\displaystyle\tilde{S}_{3}(z)$
$\displaystyle=2G\bigl{(}-1,-1,-\tfrac{1}{2};z\bigr{)}-3G\bigl{(}-1,0,-\tfrac{1}{2};z\bigr{)}-3G\bigl{(}0,-1,-\tfrac{1}{2};z\bigr{)}+4G\bigl{(}0,0,-\tfrac{1}{2};z\bigr{)}$
$\displaystyle\quad-2H(-1,-1,0;z)+2H(-1,0,-1;z)-3H(-1,0,0;z)-2H(-1,1,0;z)+H(0,-1,0;z)$
$\displaystyle\quad-H(0,0,-1;z)-H(0,0,1;z)-2H(1,-1,0;z)-\frac{\ln^{3}(z)}{12}-2H(-1,0;z)\ln(1-z)$
$\displaystyle\quad+\ln(2)\Bigl{[}-\frac{\pi^{2}}{6}+G\bigl{(}-1,-\tfrac{1}{2};z\bigr{)}-G\bigl{(}0,-\tfrac{1}{2};z\bigr{)}-2H(-1,0;z)+\frac{\ln^{2}(z)}{2}\Bigr{]}$
$\displaystyle\quad+\ln(1+z)\Bigl{[}-\frac{5\pi^{2}}{12}-G\bigl{(}-1,-\tfrac{1}{2};z\bigr{)}+G\bigl{(}0,-\tfrac{1}{2};z\bigr{)}\Bigr{]}$
$\displaystyle\quad+\ln(z)\Bigl{[}\frac{\pi^{2}}{4}+2G\bigl{(}-1,-\tfrac{1}{2};z\bigr{)}-2G\bigl{(}0,-\tfrac{1}{2};z\bigr{)}+2H(-1,0,z)\Bigr{]}+\frac{9\zeta_{3}}{4}\,.$
(S55)

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.

Figure S2: Breakdown of Fig. 3 into partonic channels. Shown are the
$\mathcal{O}(\alpha_{s}^{2})$ contribution to
$\mathrm{d}\sigma/\mathrm{d}x_{a}\mathrm{d}x_{b}$ predicted for $x_{a}\to 1$
by Eq. (14) (blue), the full result from Vrap (red). In all cases, their
difference (green) vanishes like a power as $1-x_{a}\to 0$, as it must. The
$gg$ channel is power suppressed, so its full result by itself vanishes like a
power. The error bars indicate the integration uncertainties. Dashing in the
blue line indicates a negative result.

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

Figure S3: Top row: Generalized threshold approximation of the
$\mathcal{O}(\alpha_{s})$ (top left) and $\mathcal{O}(\alpha_{s}^{2})$
contribution (top right) to the Drell-Yan rapidity spectrum
$\sigma\equiv\mathrm{d}\sigma/(\mathrm{d}Q\mathrm{d}Y)$ at $\mu=Q/2$
normalized to the LO result. This is the analog of Fig. 4.
Bottom row: Convergence of the generalized and soft threshold expansions for
the $q{\bar{q}}$ (bottom left) and $qg+gq$ channels (bottom right) for
$\mu=Q/2$. This is the analog of Fig. 5.

Figure S4: Top row: Generalized threshold approximation of the
$\mathcal{O}(\alpha_{s})$ contribution to the $gg\to H$ rapidity spectrum
$\sigma\equiv\mathrm{d}\sigma/\mathrm{d}Y$ at $\mu=m_{H}$ (top left) and
$\mu=m_{H}/2$ (top right) normalized to the LO result. This is the analog of
Fig. 4.
Bottom row: Convergence of the generalized and soft threshold expansions for
the $gg$ (bottom left) and $gq+qg$ channels (bottom right) for $\mu=m_{H}/2$.
This is the analog of Fig. 5.

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$.