# Zero-bin subtraction and the $q_{\mathrm{T}}$ spectrum beyond leading power

Giancarlo Ferrera(a,b), Wan-Li Ju(b), Marek Schönherr(c)

###### Abstract

In this paper, we present an algorithm to construct the $q_{\mathrm{T}}$
distribution at NLO accuracy to arbitrary power precision, including the
assembly of suitable zero-bin subtrahends, in a mathematically well-defined
way for a generic choice of rapidity-divergence regularisation prescription.
In its derivation, we divide the phase space into two sectors, the interior of
the integration domain as well as the integration boundary, which we include
here for the first time. To demonstrate the applicability and usefulness of
our algorithm, we calculate the N2LP corrections for Higgs hadroproduction for
the first time. We observe that our approximate N2LP-accurate $q_{\mathrm{T}}$
spectra replicate the asymptotic behaviour of the full QCD calculation to a
much better degree than the previously available results, both within the
$q_{\mathrm{T}}\to 0$ limit and in the large-$q_{\mathrm{T}}$ domain for all
the involved partonic processes. While playing a minor role at larger
transverse momenta, we show that the newly incorporated boundary contribution
plays a vital role in the $q_{\mathrm{T}}\to 0$ limit, where any subleading
power accuracy would be lost without them. In particular, our N2LP-accurate
$q_{\mathrm{T}}$ expansion can approximate the exact $q_{\mathrm{T}}$
distribution up to $q_{\mathrm{T}}\approx 30\,\text{GeV}$ at the percent level
for rapidities $|Y_{H}|\lesssim 3$.

###### Contents

  1. 1 Introduction
  2. 2 Theoretical framework
    1. 2.1 Analysis of the NLO fixed order results
    2. 2.2 Prerequisites for the power expansion
    3. 2.3 Power corrections on the boundary
    4. 2.4 Power corrections over the interior domain
      1. 2.4.1 Asymptotic expansion in momentum space
      2. 2.4.2 Rapidity regularisation and zero-bin subtraction
      3. 2.4.3 Comparison with existing results
    5. 2.5 Discussion and extrapolation
  3. 3 Implementation up to N2LP
    1. 3.1 Power expansion with momentum cutoffs
    2. 3.2 Power expansion with the pure rapidity regulator
    3. 3.3 Power expansion with the exponential rapidity regulator
  4. 4 Numerical results
  5. 5 Conclusions
  6. A Application and adaptation of the dissipative regulators
    1. A.1 The $\Delta$-regulator
    2. A.2 The $\eta$-regulator
  7. B Comparison of the NLP results with the previous literature
  8. C Small $q_{\mathrm{T}}$ expansion for multi-boson hadroproduction

##  1 Introduction

Differential observables of the Drell-Yan processes and Higgs hadroproduction
play paramount roles in precisely determining the Standard Model (SM) input
parameters and probing New Physics scenarios. They have therefore drawn
extensive experimental and theoretical attention in the recent decades. The
latest measurements of the Drell-Yan processes have been carried out by the
ATLAS [1, 2, 3, 4], CMS [5, 6, 7, 8, 9] and LHCb [10] collaborations at the
LHC at colliding energies of $\sqrt{s}=7$, $8$, and $13\,\text{TeV}$, as well
as the CDF and DØ collaborations at $\sqrt{s}=1.96~{}$TeV at the Tevatron [11,
12, 13, 14]. Inclusive and differential fiducial cross section measurements
for Higgs boson production, on the other hand, have been presented in [15, 16,
17, 18, 19, 20, 21, 22, 23, 24]. Along with the progress in experimental
precision, strides were also made in the fixed-order calculations. The QCD
corrections are known up to third order for both the Drell-Yan processes [25,
26, 27, 28, 29, 30, 31, 32, 33, 34] and Higgs production [35, 36, 37, 38, 39,
40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53]. Electroweak
corrections have mostly been studied for Drell-Yan production [54, 55, 56, 57,
58, 59, 60], but are known for Higgs production as well [61, 62].

Of the differential observables, the transverse momentum distribution of the
colourless final state is of the particular concern in this work. Even though
fixed order calculations are able to produce reliable predictions in the
majority of the phase space, substantial corrections can emerge from higher
perturbative order within the low $q_{\mathrm{T}}$ regime, as a result of soft
and collinear radiation, where the bulk of the cross section resides. Hence,
in order to arrive at a sufficiently convergent result in perturbation theory,
a variety of the resummation techniques that are capable of exponentiating the
most singular, so-called leading power (LP), behaviour of the $q_{\mathrm{T}}$
distribution at every order have been developed in the past decades. Examples
are the CSS formalism [63, 64, 65, 66, 67, 68], the momentum space resummation
[69, 70, 71, 72, 73], and the SCET-based analyses [74, 75, 76, 77, 78, 79,
80]. The subsequent resummation improved $q_{\mathrm{T}}$ distributions have
been evaluated up to approximate N4LL accuracy [74, 81, 82, 83, 84, 85, 86,
87, 73, 88, 72, 89, 90, 91, 92, 93, 94, 32, 95] for the Drell-Yan process and
N3LL′ in the case of Higgs production [96, 67, 89, 97, 98, 71, 99, 73, 100,
101, 50, 53] as to the Higgs production. Recently, to further improve the
precision as well as the phase space coverage of the $q_{\mathrm{T}}$
resummation, the factorisation pattern of the subleading power contribution
has been investigated within both the non-perturbative region
$\Lambda_{H}\lesssim q_{\mathrm{T}}\ll Q_{H}$ [102, 103, 104, 105, 106, 107,
108, 109, 110, 111, 112] and the perturbative zone $\Lambda_{H}\ll
q_{\mathrm{T}}\ll Q_{H}$ [113, 114, 115], where $\Lambda_{H}$ and $Q_{H}$
characterise the hadronic and hard scales in the process under consideration,
respectively. At present, a resummation of those subleading power corrections
has not been accomplished yet beyond LL.

Similarly, subleading power corrections also play a role in the limit
$q_{\mathrm{T}}\to 0$. Here, the LP approximation can only recover the most
singular behaviour of the exact distribution at any given order in a fixed-
order expansion, leaving behind integrable (but numerically problematic)
singularities and constant terms. These non-vanishing remainders are of
particular concern in $q_{\mathrm{T}}$-based subtraction and slicing methods
for higher-order calculations [73, 99, 92]. Related results on power
suppressed contribution can be found from [116, 117, 118, 90, 119] concerning
the slicing subtraction method, the endpoint singularity in the Higgs
production and decay [120, 121, 122, 123, 124, 125, 126, 127], the event
shapes in the leptonic and hadron colliders [128, 129, 130, 131, 132, 133,
134, 135], and the threshold resummation [136, 137, 138, 139, 140, 141, 142,
143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154].

To derive the factorisation formulae beyond the leading power approximation,
one of the prerequisites is to consistently combine the contributions from the
constituent dynamic regions, generally comprising the hard, collinear, and
soft modes [102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114,
115].111 The irrelevance of the Glauber vertices is presumed in the subleading
power factorisation of [102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112,
113, 114, 115]. Up to now, a rigorous proof of the cancellation of the Glauber
contributions beyond the leading power approximation is still absent, to our
best knowledge. The leading power discussion on inclusive observables can be
found in [155, 156, 157, 158].  Although these dynamic regions are well
defined in given segments of the phase space or the loop integrals, in
practical calculations one often extrapolates their contributions from their
respective intrinsic domains to the entire integration range, necessitating a
robust systematic subtraction process to remove their overlap. This procedure
is frequently referred to as the _zero-bin subtraction_ [159, 160, 161, 162],
the soft subtraction [163, 164, 160, 161, 165], or the overlap reduction [166,
167, 114]. At leading power, zero-bin subtraction proceeds with deducting the
soft contribution from the collinear sectors, and has been demonstrated to be
valid up to N3LO [168, 169, 79, 80, 170, 171, 172, 173]. Similar conclusions
can also be found in the analysis of the subleading angular coefficients
(suppressed by a factor of $q_{\mathrm{T}}/Q_{H}$ w.r.t. leading singular
terms) in the Drell-Yan process and the semi-inclusive deep-inelastic
scatterings [106, 107, 108, 110]. However, for the subleading inclusive
observables of the Drell-Yan process the zero-bin subtraction becomes more
involved, as it comprises mixed contributions in the overlapping area from the
lowest power accuracy up to NLP [167, 114].

Even though such recipes to construct zero-bin subtractions can establish a
factorisation at the NLP level, a robust algorithm that is capable of
generating zero-bin subtrahends at an arbitrary power and choice of a
rapidity-divergence regularisation prescription is still missing.
Investigations towards this goal were pioneered by [174], in which a
mathematically well-defined derivation was delivered for the power series of
the one-loop integrals in a variety of kinematical limits. Therein, the
analysis of the massive Sudakov form factor is intimately related to the
asymptotic behaviour of the $q_{\mathrm{T}}$ distribution. Their power
expansion begins with the introduction of a group of auxiliary cutoff scales
along the integration path, from which the relevant scales of the process are
well separated and the defining dynamic regions can be prescribed accordingly.
Within these regions, the expansions are straightforward and always well-
defined. This is then followed by a recombination procedure to lift all
dependences on the auxiliary scales, thereby extrapolating the expanded
integrands from their intrinsic domains to the entire integration range.
During this procedure, however, a set of constituents emerge, consisting of
doubly and triply expanded integrands in line with their respective distinct
scaling laws which account for the overlapping contributions amongst different
modes and turn out to be in part the zero-bin subtrahends. While integrals
with a single definite scaling rule for every integration variable are
subjected only to a single expansion, multiple expansions can induce ambiguous
interpretations of the resulting integral as the integration variables in
general observe distinct scaling rules for each dynamic mode. Subsequently,
this ambiguity limits the applicability of the formalism in [174] to only a
special range of rapidity regularisation schemes, such as the analytic
regulator [74, 76] and the pure rapidity regulator [113, 175]. When using more
generic regularisation schemes, e.g. the exponential regulator [79, 80],
doubly and triply expanded integrals instead call for an unambiguous counting
rule at each power accuracy to yield results that are independent of the
scheme. As a matter of fact, any reliance on the choice of rapidity regulators
in carrying out the power expansion will hinder the establishment of the
rapidity renormalisation group [77], which, akin to the renormalisation group
governing the virtuality divergence, requires the equivalence of the power
series resulting from various rapidity-divergence regularisation approaches.

In the following, we will make use of the NLO $q_{\mathrm{T}}$ distribution of
the process $pp\to H+X$ as a demonstrative example to accomplish this goal.
Since our derivation mainly concerns the pattern of the denominators of the
squared transition amplitudes, the conclusion here can also be generalised to
the analysis of the Drell-Yan processes, $pp\to V+X$ with $V=W,Z$. We commence
by following the strategy outlined in [174] to categorise the phase space so
as to detach the emerging scales and thereby carry out the expansion in
momentum space. Afterwards, instead of extrapolating the expanded integrands
on the cumulative level, in this work we shift the auxiliary boundaries at a
particular power accuracy only while maintaining the equality to the power
series derived from the momentum space. As a result, a set of multiple
expanded integrals emerge to balance the overextended collinear sectors power
by power. We use these to construct the zero-bin subtrahends and thereby
retrieve the scalings for the occurring integration variables. We will
demonstrate that the zero-bin subtrahends proposed in this paper are
straightforwardly applicable to all rapidity regulators that preserve their
expressions before and after asymptotic expansions, including the analytic
regulator [74, 76], the exponential regulator [79, 80], and the pure rapidity
regulator [113, 175]. In addition, with a few adaptations, our algorithm
presented here can also be generalised to use the $\Delta$-regulator [162,
165, 169, 168, 176] and the $\eta$-regulator [77, 78]. This accommodates the
majority of the rapidity-divergence regularisation schemes in use at this
time.

Our paper is organised as follows. In Sec. 2 we develop our framework to
compute, in a general way, all power corrections to NLO accuracy. To this end,
we will use Sec. 2.1 and 2.2, taking the example of Higgs hadroproduction at
LHC, to review the fixed order calculation on the $q_{\mathrm{T}}$ spectrum
and thereby categorise the phase space into two sectors according to the
origins of the power corrections. The first sector, associated with the
asymptotic expansion of the boundary condition, will be investigated in Sec.
2.3. The other one takes in the bulk of phase space and thus forms the main
concern of this paper. We will devote Sec. 2.4 to elaborate on its expansion
procedure, in particular highlighting the structural pattern of the emerging
zero-bin subtrahends. At last, in Sec. 2.5, we combine the relevant
ingredients and present the power series of the $q_{\mathrm{T}}$ spectrum of
our example process $pp\to H+X$. Sec. 3 then is dedicated to the application
of the framework derived in Sec. 2, where the power corrections will be
appraised up to N2LP for the first time. We calculate it in three different
approaches, including the momentum cutoff scales, the homogenous
regularisation scheme, and the inhomogeneous regulator, to regulate the
rapidity divergences. We will demonstrate that after appropriate combination,
all three methods result in the same analytic expressions at least up to N2LP,
echoing the rapidity renormalisation group [77]. Finally, we will scrutinise
our N2LP results in Sec. 4 by comparing to the full QCD calculation.
Eventually, we summarise our findings in Sec. 5.

##  2 Theoretical framework

In this section, we introduce the formalism we derive the next-to-leading and
next-to-next-to-leading power expansions, dubbed NLP and N2LP respectively. We
will use the case of Higgs production at hadron colliders, $p+p\to H+X$, as an
illustrative example. The results we derive, however, are similarly applicable
to the production of any other colour-neutral final state.

###  2.1 Analysis of the NLO fixed order results

To introduce our notation, we start our analysis by recalling the general
expression for the differential cross section of the process $p+p\to H+X$ in
proton-proton collisions. According to the QCD factorization theorem [177],
the differential $q_{\mathrm{T}}$ spectra at NLO can be calculated as follows,

$\displaystyle\frac{{\rm d}\sigma_{H}}{{\rm d}Y_{H}\,{\rm
d}q_{\mathrm{T}}^{2}}\;=\;$
$\displaystyle\frac{1}{16\pi\,s^{2}}\sum_{i,j}\int^{k_{+}^{\mathrm{max}}}_{0}\\!\\!\\!\\!{\rm
d}k_{+}\,\int^{k_{-}^{\mathrm{max}}}_{0}\\!\\!\\!\\!{\rm
d}k_{-}\,\delta(k_{+}k_{-}-q_{\mathrm{T}}^{2})\;\frac{f_{i/n}(\xi_{n})}{\xi_{n}}\,\frac{f_{j/{\bar{n}}}(\xi_{{\bar{n}}})}{\xi_{{\bar{n}}}}\;\overline{\sum_{\mathrm{col},\mathrm{pol}}}\big{|}\mathcal{M}(i+j\to
H+k)\big{|}^{2}\,$ $\displaystyle\;=\;$
$\displaystyle\frac{1}{16\pi\,s^{2}}\sum_{i,j}\int^{k_{+}^{\mathrm{max}}}_{k_{+}^{\mathrm{min}}}\,\frac{{\rm
d}k_{+}}{k_{+}}\frac{f_{i/n}(\xi_{n})}{\xi_{n}}\frac{f_{j/{\bar{n}}}(\xi_{{\bar{n}}})}{\xi_{{\bar{n}}}}\;\overline{\sum_{\mathrm{col},\mathrm{pol}}}\big{|}\mathcal{M}(i+j\to
H+k)\big{|}^{2}\,,$ (2.1)

where $\sqrt{s}$ denotes the collider centre-of-momentum energy of the
incoming protons. $Y_{H}$ and $q_{\mathrm{T}}$ stand for the rapidity and
transverse momentum of the Higgs boson measured in the laboratory reference
frame, respectively. The boundaries of the phase space integral are given by
the variables $k_{+}^{\mathrm{max}}$ and $k_{+}^{\mathrm{min}}$, with

$\displaystyle k_{+}^{\mathrm{max}}=\sqrt{s}-m_{\mathrm{T}}e^{-Y_{H}}\,,\qquad
k_{+}^{\mathrm{min}}=\frac{q_{\mathrm{T}}^{2}}{k_{-}^{\mathrm{max}}}=\frac{q_{\mathrm{T}}^{2}}{\sqrt{s}-m_{\mathrm{T}}e^{+Y_{H}}}\;.$
(2.2)

Therein, $m_{\mathrm{T}}^{2}\equiv m_{H}^{2}+q_{\mathrm{T}}^{2}$ is the
transverse mass of the Higgs boson, $m_{H}$ being the mass of the Higgs boson.

In the integrand of Eq. (2.1), $f_{i/n}$ and $f_{j/{\bar{n}}}$ represent the
parton distribution functions (PDFs) for the partons $i$ and $j$ in the
protons $n$ and ${\bar{n}}$, respectively, with the argument
$\xi_{n({\bar{n}})}$ encoding the momentum fractions carried by the incident
particle $i(j)$ with respect to its mother proton $n({\bar{n}})$. In this
paper, we will use the 5 flavour scheme, and correspondingly, the partons
$\\{i,j,k\\}\in\\{u,\bar{u},d,\bar{d},c,\bar{c},s,\bar{s},b,\bar{b},g\\}$ will
be considered massless and taken fully into account in the following analysis.

$|\mathcal{M}|^{2}$ denotes the squared amplitudes for the contributing
partonic channel $i+j\to H+k$. Since the main focus of this paper is on the
region where $q_{\mathrm{T}}$ approaches, but is not equal to, zero, only tree
level contributions will participate in Eq. (2.1). In the large top mass
limit, the tree level amplitudes for the processes $i+j\to H+k$ can be
calculated in the low energy Higgs effective field theory (HEFT) [178, 179,
180, 181, 182, 183, 184]. Having summed and averaged, as appropriate, over all
the colour and polarisation configurations, the results read [185, 186],

$\begin{split}\overline{\sum_{\mathrm{col},\mathrm{pol}}}\big{|}\mathcal{M}(g+g\to
H+g)\big{|}^{2}&=\,\frac{3\lambda_{t}}{32}\,\frac{m_{H}^{8}+s_{ij}^{4}+s_{ik}^{4}+s_{jk}^{4}}{s_{ij}s_{ik}s_{jk}}\,,\\\
\overline{\sum_{\mathrm{col},\mathrm{pol}}}\big{|}\mathcal{M}(g+q(\bar{q})\to
H+q(\bar{q}))\big{|}^{2}&=\,-\frac{\lambda_{t}}{24}\,\frac{s^{2}_{ij}+s^{2}_{ik}}{s_{jk}}\,,\\\
\overline{\sum_{\mathrm{col},\mathrm{pol}}}\big{|}\mathcal{M}(q+\bar{q}\to
H+g)\big{|}^{2}&=\,\frac{\lambda_{t}}{9}\,\frac{s^{2}_{ik}+s^{2}_{jk}}{s_{ij}}\,,\\\
\end{split}$ (2.3)

where the coupling parameter
$\lambda_{t}\equiv{4\alpha_{s}^{3}C_{t}^{2}}/{(9\pi v^{2})}$ collects the
strong coupling constant $\alpha_{s}$ as well as the vacuum expectation value
$v$, and the Wilson coefficient $C_{t}$ of the HEFT governing the Higgs-gluon
interaction. Up to now, $C_{t}$ has been calculated up to the four-loop order
in the literature [182, 187, 188, 184, 189, 190, 191]. In this work, we only
take into account its LO contribution in accordance with the perturbative
accuracy of the amplitudes of Eq. (2.3). The missing expressions for
$q(\bar{q})g$ and $\bar{q}q$ initial states can be derived from the above by
swapping the roles of $i$ and $j$.

Further, the expressions in Eq. (2.3) introduce the scalar products of the
momenta of the initial and final particles,

$\begin{split}s_{ij}=&\phantom{+}2p_{i}\cdot
p_{j}=\left(k_{+}+m_{\mathrm{T}}e^{-Y_{H}}\right)\,\left(k_{-}+m_{\mathrm{T}}e^{+Y_{H}}\right)\,,\\\
s_{ik}=&-2p_{i}\cdot
p_{k}=-\left(k_{+}+m_{\mathrm{T}}e^{-Y_{H}}\right)k_{-}\,,\\\
s_{jk}=&-2p_{j}\cdot
p_{k}=-\left(k_{-}+m_{\mathrm{T}}e^{+Y_{H}}\right)k_{+}\,.\end{split}$ (2.4)

Therein, $k_{\pm}$ denotes the light cone component of the momentum of the
emitted parton $k$. To be precise, with the help of the light-like reference
vectors $n^{\mu}$ and ${\bar{n}}^{\mu}$, aligned with incoming beams and
satisfying $n\cdot{\bar{n}}=2$, it is defined as

$\displaystyle p_{k}^{\mu}=\frac{p_{k}\cdot
n}{2}{\bar{n}}^{\mu}+\frac{p_{k}\cdot{\bar{n}}}{2}{n}^{\mu}+p_{k,\bot}^{\mu}\equiv\frac{k_{-}}{2}{\bar{n}}^{\mu}+\frac{k_{+}}{2}{n}^{\mu}+k_{\bot}^{\mu}\equiv\big{[}k_{-},k_{+},k_{\bot}\big{]}\,.$
(2.5)

In this notation, the momentum fractions $\xi_{n}$ and $\xi_{{\bar{n}}}$ of
Eq. (2.1) can be written as

$\displaystyle\xi_{n}=\frac{k_{+}+m_{\mathrm{T}}e^{-Y_{H}}}{\sqrt{s}}\qquad\text{and}\qquad\xi_{{\bar{n}}}=\frac{k_{-}+m_{\mathrm{T}}e^{+Y_{H}}}{\sqrt{s}}\;.$
(2.6)

In the following sections, we will investigate the asymptotic properties of
Eq. (2.1) in the vicinity of $q_{\mathrm{T}}=0$, adhering to the following
definition of the expansion accuracy, counting powers of $q_{\mathrm{T}}^{2}$,
throughout,

$\begin{split}\frac{{\rm d}\sigma_{H}}{{\rm d}Y_{H}\,{\rm
d}q_{\mathrm{T}}^{2}}=&\,\underbrace{\sum_{m}\,\frac{\Delta_{\text{LP}}^{(m)}\,(L_{H})^{m}}{q_{\mathrm{T}}^{2}}}_{\text{LP}}\,+\,\underbrace{\sum_{m}\,\Delta_{\text{NLP}}^{(m)}\,(L_{H})^{m}}_{\text{NLP}}\,+\,\underbrace{\sum_{m}\,q_{\mathrm{T}}^{2}\,\Delta_{\text{N${}^{2}$LP}}^{(m)}\,(L_{H})^{m}}_{\text{N${}^{2}$LP}}\,+\,\dots\,.\\\
\end{split}$ (2.7)

Therein, the $q_{\mathrm{T}}$ differential distribution in Eq. (2.1) is
expanded in the low $q_{\mathrm{T}}$ domain. At a given power precision, e.g.
Nω+1LP, all non-logarithmic $q_{\mathrm{T}}$ dependences have been collected
in the pre-factor $q_{\mathrm{T}}^{2\omega}$ with $(\omega\geq-1)$. They leave
behind the logarithmic terms $L_{H}\equiv\ln[q_{\mathrm{T}}^{2}/m_{H}^{2}]$
and the coefficient functions
$\Delta_{\mathrm{N}^{\omega+1}\mathrm{LP}}^{(m)}\equiv\Delta_{\mathrm{N}^{\omega+1}\mathrm{LP}}^{(m)}[m_{H},\,s,\,Y_{H};\;f_{i/n},\,f_{j/{\bar{n}}}]$
which accounts for the convolution of the PDFs $f_{i/n}$ and $f_{j/{\bar{n}}}$
as well as the functions of the hard scales $s$ and $m_{H}$. According to this
definition, the LP contribution is expected to accommodate the most singular
behaviour in the low $q_{\mathrm{T}}$ regime, while the higher power
corrections experience progressively stronger suppression factor
$q_{\mathrm{T}}^{2(\omega-1)}$. Correspondingly, $\Delta_{\text{LP}}$ is known
extremely precisely, up to approximate fourth order [192, 193, 194, 195, 173,
172, 196, 197, 198, 176, 170, 171, 100, 188, 184, 199, 200, 201, 202, 203,
204, 205, 190, 189, 206, 207], while the $\Delta_{\text{NLP}}$ has only been
computed recently to NLO accuracy [113]. All further higher-power corrections
are hitherto unknown, and we will use the remainder of this section to
introduce a mathematically well-defined approach to develop a framework to
compute the $\Delta_{\mathrm{N}^{\omega+1}\mathrm{LP}}$ for all $\omega\geq-1$
at NLO accuracy before deriving the respective NLO-correct expression for
$\Delta_{\text{N${}^{2}$LP}}$ explicitly.

In the following subsections we will be evaluating the coefficients
$\Delta_{\mathrm{N}^{\omega+1}\mathrm{LP}}^{(m)}$. In Sec. 2.2, we will derive
additional power corrections originating from the kinematics of the process
and thereby categorise the phase space integral of Eq. (2.1) into two sectors.
The first sector encompasses the domains in the vicinity of
$k^{\mathrm{max}}_{+}$ and $k^{\mathrm{max}}_{-}$, accounting for the power
suppressed contributions induced by the boundary conditions of the phase space
integral, see Eq. (2.2). We will make use of Sec. 2.3 to discuss the
asymptotic properties from this region. The second sector consists of the
remaining regions, entailing a variety of scales along the integration path
and is therefore the main concern of this paper. We will elaborate in Sec. 2.4
on the power expansion of this contribution. Eventually, with the results of
the two, we are capable of deriving the power series for the $q_{\mathrm{T}}$
spectrum in Sec. 2.5.

###  2.2 Prerequisites for the power expansion

With the light-cone coordinates defined in Eq. (2.5), one can recast the
$q_{\mathrm{T}}$ spectrum of eq. (2.1) into a following form,

$\begin{split}\frac{{\rm d}\sigma_{H}}{{\rm d}Y_{H}\,{\rm
d}q_{\mathrm{T}}^{2}}\,=&\;\lambda^{2}_{t}\,\sum_{\kappa}\,\sum_{\\{\beta\\}}\,\sum_{\rho,\sigma}\,\mathcal{H}^{\rho,\sigma}_{[\kappa],\\{\beta\\}}\left(m_{\mathrm{T}},Y_{H},s\right)\,\\\
&\times\;\int^{k_{+}^{\mathrm{max}}}_{k_{+}^{\mathrm{min}}}\frac{{\rm
d}k_{+}}{k_{+}}\;(k_{-})^{\rho}\,(k_{+})^{\sigma}\,{F}^{(0)}_{i/n,\beta_{n}}\left(k_{+}+m_{\mathrm{T}}e^{-Y_{H}}\right)\,{F}^{(0)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(k_{-}+m_{\mathrm{T}}e^{+Y_{H}}\right)\,,\end{split}$
(2.8)

where $\mathcal{H}^{\rho,\sigma}_{[\kappa],\\{\beta\\}}$ is a function of the
hard scales $m_{\mathrm{T}}$ and $\sqrt{s}$. The contributing partonic
channels are denoted by the subscript
$\kappa=\\{gg,gq(\bar{q}),q(\bar{q})g,q\bar{q},\bar{q}q\\}$. To be precise,

$\begin{split}\overline{\sum_{\mathrm{col},\mathrm{pol}}}\big{|}\mathcal{M}_{[\kappa]}\big{|}^{2}\,\equiv&\;{16\lambda_{t}^{2}\pi
s}\,\sum_{\\{\beta\\}}\sum_{\rho,\sigma}\,\frac{(k_{+})^{\sigma}}{(k_{+}+m_{\mathrm{T}}e^{-Y_{H}})^{\beta_{n}-1}}\,\frac{(k_{-})^{\rho}}{(k_{-}+m_{\mathrm{T}}e^{Y_{H}})^{\beta_{\bar{n}}-1}}\,\mathcal{H}^{\rho,\sigma}_{[\kappa],\\{\beta\\}}\left(m_{\mathrm{T}},Y_{H},s\right)\,.\end{split}$
(2.9)

Further, the integrand of Eq. (2.8) contains the derivatives of the scaled
PDFs of the light-cone components $(k_{\pm}+m_{\mathrm{T}}e^{\mp Y_{H}})$,
generically defined as

$\displaystyle{F}^{(\alpha_{n})}_{i/n,\beta_{n}}\left(Q_{n}\right)$
$\displaystyle\equiv\frac{\partial^{\alpha_{n}}}{\partial
Q_{n}^{\alpha_{n}}}\left[\frac{f_{i/n}\left({Q_{n}}/{\sqrt{s}}\right)}{Q_{n}^{\beta_{n}}}\right]\,,\quad\quad{F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}\left(Q_{\bar{n}}\right)\equiv\frac{\partial^{\alpha_{\bar{n}}}}{\partial
Q_{\bar{n}}^{\alpha_{\bar{n}}}}\left[\frac{f_{j/{\bar{n}}}\left(Q_{\bar{n}}/{\sqrt{s}}\right)}{Q_{\bar{n}}^{\beta_{\bar{n}}}}\right]\,,$
(2.10)

where ${Q}_{n}$ and $Q_{\bar{n}}$ represent the two kinematic variables,
$\alpha_{n}$ and $\alpha_{\bar{n}}$ determine the rank of the derivative, and
$\beta_{n}$ and $\beta_{\bar{n}}$ are the powers to which they scale their
respective PDFs. Collectively, we denote the latter $\\{\alpha\\}$ and
$\\{\beta\\}$.

For the results derived in full QCD, see Eq. (2.1), only the zero-rank scaled
PDFs are involved. Hence, only the $\alpha_{n}=\alpha_{\bar{n}}=0$ case is
present in Eq. (2.8). The $\alpha_{n},\alpha_{\bar{n}}>0$ expressions,
however, will enter in the higher-power $q_{\mathrm{T}}$ expansion.

In the following, we will concern ourselves with the analytic properties of
the functions ${F}^{(\alpha_{n})}_{i/n,\beta_{n}}\left(Q_{n}\right)$ and
${F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}\left(Q_{\bar{n}}\right)$.
However, due to the fact that the PDFs $f_{i/n}$ and $f_{j/{\bar{n}}}$ have
their origin in the non-perturbative regime of QCD and, thus, their analytic
expressions are not fully established, we make the following assumptions:

  1. (a)

For a given $\alpha_{n},\alpha_{\bar{n}}>0$ the
${F}^{(\alpha_{n})}_{i/n,\beta_{n}}\left(Q_{n}\right)$ and
${F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}\left(Q_{\bar{n}}\right)$,
as well as their arbitrary order derivatives, exist and are bounded over the
domain $\Lambda_{\mathrm{H}}\leq{Q}_{n({\bar{n}})}\leq\sqrt{s}$ for all
$\beta_{n}$ and $\beta_{\bar{n}}$, where $\Lambda_{\mathrm{H}}$ represent a
hadronisation scale .

  2. (b)

${F}^{(\alpha_{n})}_{i/n,\beta_{n}}\left(Q_{n}+\delta Q_{n}\right)$ and
${F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}\left(Q_{\bar{n}}+\delta
Q_{\bar{n}}\right)$ have convergence radii
$|\delta{Q}_{n({\bar{n}})}|<{Q}_{n({\bar{n}})}$ as long as
$\left({Q}_{n({\bar{n}})}+\delta{Q}_{n({\bar{n}})}\right)$ and
${Q}_{n({\bar{n}})}$ are both well within the interval
$\left[\Lambda_{\mathrm{H}},\sqrt{s}\right]$ .

Please note, if the PDFs are fitted in a basis of Chebyshev polynomials, as
proposed in [208, 209, 210], and thereby $f_{i/n}$ and $f_{j/{\bar{n}}}$ can
be then expressed in terms of $\ln^{k}(\xi_{n})/\xi_{n}^{h}$ and
$\ln^{k}(\xi_{\bar{n}})/\xi_{\bar{n}}^{h}$, respectively, the second
assumption is immediately fulfilled.

Equipped with the above ansatz, we are able to expand Eq. (2.8) in the
vicinity of $q_{\mathrm{T}}=0$. To this end, we start with the expansion of
the kinematic variable $m_{\mathrm{T}}$ in the hard coefficient function
$\mathcal{H}^{\rho,\sigma}_{[\kappa],\\{\beta\\}}$. In the small
$q_{\mathrm{T}}$ domain where $q_{\mathrm{T}}\ll m_{H}$, we can perform an
expansion of the transverse mass $m_{\mathrm{T}}$ of the Higgs boson around
its invariant mass $m_{H}$,

$\displaystyle
m_{\mathrm{T}}\,=\,m_{H}\,\cdot\,\sum_{h=0}^{\infty}\,\frac{1}{h!}\,\frac{\Gamma[\tfrac{3}{2}]}{\Gamma[\frac{3}{2}-h]}\,\left(\frac{q_{\mathrm{T}}^{2}}{m_{H}^{2}}\right)^{h}\,.$
(2.11)

We can use this result to expand the hard coefficient
$\mathcal{H}^{\rho,\sigma}_{[\kappa],\\{\beta\\}}$ around $q_{\mathrm{T}}=0$,
giving

$\begin{split}\mathcal{H}_{[\kappa],\\{\beta\\}}^{\rho,\sigma}&(m_{\mathrm{T}},Y_{H},s)\\\
=&\,\sum_{h=0}^{\infty}\frac{(m_{\mathrm{T}}-m_{H})^{h}}{h!}\,\left\\{\frac{\partial^{h}}{\partial
m_{\mathrm{T}}^{h}}\,\mathcal{H}_{[\kappa],\\{\beta\\}}^{\rho,\sigma}(m_{\mathrm{T}},Y_{H},s)\bigg{|}_{m_{\mathrm{T}}\to
m_{H}}\right\\}\,\\\
=&\,\sum_{h,l=0}^{\infty}\,\sum_{g=0}^{h}\,\frac{(-1)^{h-g}}{g!\,l!\,(h-g)!}\,\frac{\Gamma[\frac{g}{2}+1]}{\Gamma[\frac{g}{2}-l+1]}\,\left(\frac{q_{\mathrm{T}}^{2}}{m^{2}_{H}}\right)^{l}\,m_{H}^{h}\,\mathcal{H}_{[\kappa],\\{\beta\\}}^{(h),\rho,\sigma}(m_{H},Y_{H},s)\,,\end{split}$
(2.12)

where we have introduced
$\mathcal{H}_{[\kappa],\\{\beta\\}}^{(h),\rho,\sigma}(m_{H},Y_{H},s)\equiv\frac{\partial^{h}}{\partial
m_{\mathrm{T}}^{h}}\,\mathcal{H}_{[\kappa],\\{\beta\\}}^{\rho,\sigma}(m_{\mathrm{T}},Y_{H},s)\bigg{|}_{m_{\mathrm{T}}\to
m_{H}}$. An analogous expansion can also be applied to the derivatives of the
scaled PDFs ${F}^{(\alpha_{n})}_{i/n,\beta_{n}}$ and
${F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}$. To this end, we re-
express their arguments below in order to isolate all $q_{\mathrm{T}}$
dependences,

$\displaystyle k_{\pm}+m_{\mathrm{T}}e^{\mp Y_{H}}\to\left[k_{\pm}+m_{H}e^{\mp
Y_{H}}\right]+\left(m_{\mathrm{T}}-m_{H}\right)\,e^{\mp Y_{H}}.$ (2.13)

In proximity to $q_{\mathrm{T}}=0$, the second term on the right handed side
(r.h.s) is always smaller in magnitude with respect to the first one.
Therefore, we can expand
${F}^{(\alpha_{n})}_{i/n,\beta_{n}}(k_{+}+m_{\mathrm{T}}e^{-Y_{H}})$ and
${F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}(k_{-}+m_{\mathrm{T}}e^{+Y_{H}})$
around $(k_{+}+m_{H}e^{-Y_{H}})$ and $(k_{-}+m_{H}e^{+Y_{H}})$, respectively,

$\begin{split}F^{(0)}_{i/n,\beta_{n}}&\left(k_{+}+m_{\mathrm{T}}e^{-Y_{H}}\right)\\\
&=\,\sum_{h=0}^{\infty}\,\frac{(m_{\mathrm{T}}-m_{H})^{h}}{h!}\,(e^{-Y_{H}})^{h}\,F^{(h)}_{i/n,\beta_{n}}\left(k_{+}+m_{H}e^{-Y_{H}}\right)\\\
&=\,\sum_{h,l=0}^{\infty}\sum_{g=0}^{h}\,\frac{(-1)^{h-g}}{g!\,l!\,(h-g)!}\,\frac{\Gamma[\frac{g}{2}+1]}{\Gamma[\frac{g}{2}-l+1]}\,\left(\frac{q_{\mathrm{T}}^{2}}{m^{2}_{H}}\right)^{l}\,\left(m_{H}e^{-Y_{H}}\right)^{h}\,F^{(h)}_{i/n,\beta_{n}}\left(k_{+}+m_{H}e^{-Y_{H}}\right)\,,\\\
F^{(0)}_{j/{\bar{n}},\beta_{\bar{n}}}&\left(k_{-}+m_{\mathrm{T}}e^{+Y_{H}}\right)\\\
&=\,\sum_{h=0}^{\infty}\,\frac{(m_{\mathrm{T}}-m_{H})^{h}}{h!}\,(e^{+Y_{H}})^{h}\,F^{(h)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(k_{-}+m_{H}e^{+Y_{H}}\right)\\\
&=\,\sum_{h,l=0}^{\infty}\sum_{g=0}^{h}\,\frac{(-1)^{h-g}}{g!\,l!\,(h-g)!}\,\frac{\Gamma[\frac{g}{2}+1]}{\Gamma[\frac{g}{2}-l+1]}\,\left(\frac{q_{\mathrm{T}}^{2}}{m^{2}_{H}}\right)^{l}\,\left(m_{H}e^{+Y_{H}}\right)^{h}\,F^{(h)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(k_{-}+m_{H}e^{+Y_{H}}\right)\,.\end{split}$
(2.14)

Inserting the expressions of Eqs. (2.12)-(2.14) into Eq. (2.8), the
$q_{\mathrm{T}}$ spectrum now takes the following form,

$\displaystyle\frac{{\rm d}\sigma_{H}}{{\rm d}Y_{H}\,{\rm
d}q_{\mathrm{T}}^{2}}\,=\,\lambda_{t}^{2}\,\sum_{\omega}\,\left({q_{\mathrm{T}}^{2}}\right)^{\omega}\,\sum_{[\kappa]}\,\sum_{\\{\alpha,\beta\\}}\,\sum_{\rho,\sigma}\,\widetilde{\mathcal{H}}^{(\omega),\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}\left(m_{H},Y_{H},s\right)\,\bigg{\\{}\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}\,+\,\Delta\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}\bigg{\\}}\;.$
(2.15)

Therein,
$\widetilde{\mathcal{H}}^{(\omega),\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$
is introduced to absorb the factorials, gamma functions, as well as powers of
the hard scales $m_{H}$ and $\sqrt{s}$ emerging from Eqs. (2.12)-(2.14),

$\begin{split}\widetilde{\mathcal{H}}^{(\omega),\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}\,=&\;\left(\frac{1}{m^{2}_{H}}\right)^{\omega}\,\sum_{l_{n}=0}^{\omega}\,\sum_{l_{\bar{n}}=0}^{\omega-
l_{n}}\,\sum_{\alpha_{h}=0}^{\infty}\,\prod_{x=\\{h,n,{\bar{n}}\\}}\left\\{\sum_{g_{x}=0}^{\alpha_{x}}\frac{(-1)^{\alpha_{x}-g_{x}}}{g_{x}!\,l_{x}!\,(\alpha_{x}-g_{x})!}\,\frac{\Gamma[\frac{g_{x}}{2}+1]}{\Gamma[\frac{g_{x}}{2}-l_{x}+1]}\,\right\\}\\\
&\times\left(m_{H}e^{-Y_{H}}\right)^{\alpha_{n}}\,\left(m_{H}e^{+Y_{H}}\right)^{\alpha_{\bar{n}}}\,\left(m_{H}\right)^{\alpha_{h}}\,{\mathcal{H}}^{(\alpha_{h}),\rho,\sigma}_{[\kappa],\\{\beta\\}}\,,\end{split}$
(2.16)

where $l_{h}=\omega-l_{n}-l_{\bar{n}}\geq 0$. The $q_{\mathrm{T}}$ dependences
have been factored out from
$\widetilde{\mathcal{H}}^{(\omega),\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$
and the arguments of ${F}^{(\alpha_{n})}_{i/n,\beta_{n}}$ and
${F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}$. The remaining ones
participate through boundary conditions as defined in Eq. (2.2). To address
them as well, we divide the phase space integral of Eq. (2.15) into two
sectors, $\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$
and
$\Delta\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$,
with

$\displaystyle\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}=$
$\displaystyle\,\int^{\tilde{k}^{\mathrm{max}}_{+}}_{\tilde{k}^{\mathrm{min}}_{+}}\,\frac{{\rm
d}k_{+}}{k_{+}}(k_{-})^{\rho}\,(k_{+})^{\sigma}\,{F}^{(\alpha_{n})}_{i/n,\beta_{n}}\left(k_{+}+m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}\left(k_{-}+m_{H}e^{+Y_{H}}\right)\,,$
(2.17)
$\displaystyle\Delta\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}=$
$\displaystyle\,\left(\int_{\tilde{k}^{\mathrm{max}}_{+}}^{{k}^{\mathrm{max}}_{+}}\\!\\!\\!\\!+\int^{\tilde{k}^{\mathrm{min}}_{+}}_{{k}^{\mathrm{min}}_{+}}\right)\frac{{\rm
d}k_{+}}{k_{+}}(k_{-})^{\rho}\,(k_{+})^{\sigma}\,{F}^{(\alpha_{n})}_{i/n,\beta_{n}}\\!\\!\left(k_{+}+m_{H}e^{-Y_{H}}\right){F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}\\!\\!\left(k_{-}+m_{H}e^{+Y_{H}}\right)\,,\\!\\!\\!\\!\\!\\!$
(2.18)

by means of the following boundaries,

$\displaystyle\tilde{k}_{+}^{\mathrm{max}}=\sqrt{s}-m_{H}e^{-Y_{H}}\,,\quad\quad\tilde{k}_{+}^{\mathrm{min}}={q_{\mathrm{T}}^{2}}/(\sqrt{s}-m_{H}e^{+Y_{H}})\,.$
(2.19)

Consequently,
$\Delta\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$
contains all $m_{\mathrm{T}}$ dependences. In analysing the scaling behaviour,
we find that within the central region $e^{\pm Y_{H}}\sim\mathcal{O}(1)$, the
integration variables $k_{\pm}$ exhibit unambiguous scaling in each
integration region. Therefore, as will be illustrated in Sec. 2.3, the power
expansion of
$\Delta\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$
follows an analogous pattern to those in Eqs. (2.12)-(2.14). Nevertheless,
$\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$ behaves
differently. From its expression in Eq. (2.17), a number of different scales
have been enclosed in the integration range of
$\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$,
starting from the ultrasoft scale
${k^{\mathrm{min}}_{+}}\sim\mathcal{O}(q_{\mathrm{T}}^{2}/m_{H})$, through the
soft one $k_{\pm}\sim\mathcal{O}(q_{\mathrm{T}})$, and ending up with the hard
scale ${k^{\mathrm{max}}_{+}}\sim\mathcal{O}(m_{H})$. They all correspond to a
distinct way of expanding the integrands. To this end, in Sec. 2.4, we will
make use of a set of momentum cutoffs to unambiguously treat them.

###  2.3 Power corrections on the boundary

We have mentioned above that besides the integrand itself, the integration
boundaries are other sources for power corrections. In the following, we will
thus expand the integral
$\Delta\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$
within the central region $e^{\pm Y_{H}}\sim\mathcal{O}(1)$. As illustrated in
Eq. (2.18),
$\Delta\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$
consists of the integral of $k_{+}$ over two intervals,
${k}_{+}\in[\tilde{k}_{+}^{\mathrm{max}},{k}_{+}^{\mathrm{max}}]$ and
${k}_{+}\in[{k}_{+}^{\mathrm{min}},\tilde{k}_{+}^{\mathrm{min}}]$. In each of
them, the momentum $k^{\mu}$ of the emitted particle is subject to an
unambiguous scaling,

$\displaystyle{k}_{+}\in[\tilde{k}_{+}^{\mathrm{max}},{k}_{+}^{\mathrm{max}}]\;:\quad
k_{+}\sim\mathcal{O}(m_{H})\,,\quad
k_{-}\sim\mathcal{O}(q_{\mathrm{T}}^{2}/m_{H})\,,\quad\vec{k}_{\mathrm{T}}=-\vec{q}_{\mathrm{T}}\,,$
(2.20)
$\displaystyle{k}_{+}\in[{k}_{+}^{\mathrm{min}},\tilde{k}_{+}^{\mathrm{min}}]\;:\quad
k_{+}\sim\mathcal{O}(q_{\mathrm{T}}^{2}/m_{H})\,,\quad
k_{-}\sim\mathcal{O}(m_{H})\,,\quad\vec{k}_{\mathrm{T}}=-\vec{q}_{\mathrm{T}}\,,$
(2.21)

where the transverse recoil $\vec{k}_{\mathrm{T}}$ is determined by the
momentum conservation. Considering that the momenta $k_{\pm}$ in Eqs.
(2.20)-(2.21) are nearly corresponding to the hardest emission along the beam
direction that are kinematically allowed (in a logarithmic sense), henceforth,
we dub them the $n$-ultra-collinear and ${\bar{n}}$-ultra-collinear modes,
respectively. In the following, we will make use of those two scaling laws to
expand
$\Delta\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$
within the low $q_{\mathrm{T}}$ regime.

We first consider the $n$-ultra-collinear contribution, applying Eq. (2.20) to
Eq. (2.18) yields

$\begin{split}\Delta\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{uc}\,=\,&\int_{\tilde{k}^{\mathrm{max}}_{+}}^{{k}^{\mathrm{max}}_{+}}\,\frac{{\rm
d}k_{+}}{k_{+}}\,(k_{-})^{\rho}\,(k_{+})^{\sigma}\,{F}^{(\alpha_{n})}_{i/n,\beta_{n}}\left(k_{+}+m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}\left(k_{-}+m_{H}e^{+Y_{H}}\right)\\\
=\,&(q_{\mathrm{T}}^{2})^{\rho}\,\int_{0}^{(m_{H}-m_{\mathrm{T}})e^{-Y_{H}}}\,{\rm
d}\hat{k}_{+}\,(\hat{k}_{+}+\tilde{k}^{\mathrm{max}}_{+})^{\sigma-\rho-1}\\\
&\hskip
56.9055pt\times\,{F}^{(\alpha_{n})}_{i/n,\beta_{n}}\left(\hat{k}_{+}+\sqrt{s}\right)\,{F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}\left(\frac{q_{\mathrm{T}}^{2}}{\hat{k}_{+}+\tilde{k}^{\mathrm{max}}_{+}}+m_{H}e^{+Y_{H}}\right)\\\
=\,&\sum_{\tilde{\alpha}_{n},\tilde{\alpha}_{\bar{n}}=0}^{\infty}\,\frac{(q_{\mathrm{T}}^{2})^{\rho+\tilde{\alpha}_{\bar{n}}}}{\tilde{\alpha}_{n}!\,\tilde{\alpha}_{\bar{n}}!}\,{F}^{(\alpha_{n}+\tilde{\alpha}_{n})}_{i/n,\beta_{n}}\left(\sqrt{s}\right)\,{F}^{(\alpha_{\bar{n}}+\tilde{\alpha}_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)\,\mathcal{B}^{\rho,\sigma}_{+,\\{\alpha,\beta\\}}\,,\end{split}$
(2.22)

where in the second step the integration variable is transformed such that the
power suppressed residual terms can be extracted, i.e.

$k_{+}\,\to\,\underbrace{\quad\hat{k}_{+}\quad}_{\sim\mathcal{O}(q_{\mathrm{T}}^{2}/m_{H})}\,+\,\underbrace{\quad\tilde{k}_{+}^{\mathrm{max}}\quad}_{\sim\mathcal{O}(m_{H})}\,.$
(2.23)

Based on this scaling behaviour, the functions
${F}^{(\alpha_{n})}_{i/n,\beta_{n}}$ and
${F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}$ can be expanded
around $\sqrt{s}$ and $m_{H}e^{+Y_{H}}$, respectively, in the last step of Eq.
(2.22), thereby removed from the $\hat{k}_{+}$-integral. The remaining
$\hat{k}_{+}$ dependences are collected in the function
$\mathcal{B}^{\rho,\sigma}_{+,\\{\alpha,\beta\\}}$, defined through

$\begin{split}\mathcal{B}^{\rho,\sigma}_{+,\\{\alpha,\beta\\}}\equiv&\;\int_{0}^{(m_{H}-m_{\mathrm{T}})e^{-Y_{H}}}{\rm
d}\hat{k}_{+}\;(\hat{k}_{+}+\tilde{k}^{\mathrm{max}}_{+})^{\sigma-\rho-\tilde{\alpha}_{\bar{n}}-1}(\hat{k}_{+})^{\tilde{\alpha}_{n}}\\\
=&\;\int_{0}^{(m_{H}-m_{\mathrm{T}})e^{-Y_{H}}}{\rm
d}\hat{k}_{+}\;(\hat{k}_{+})^{\tilde{\alpha}_{n}}\,\left(\tilde{k}^{\mathrm{max}}_{+}\right)^{\sigma-\rho-\tilde{\alpha}_{\bar{n}}-1}\\\
&\hskip
28.45274pt\times\Bigg{\\{}\theta(\sigma-\rho-\tilde{\alpha}_{\bar{n}}-1)\,\sum_{\gamma=0}^{\sigma-\rho-\tilde{\alpha}_{\bar{n}}-1}\,\frac{1}{\gamma!}\,\left(\frac{\hat{k}_{+}}{\tilde{k}^{\mathrm{max}}_{+}}\right)^{\gamma}\,\frac{\Gamma[\sigma-\rho-\tilde{\alpha}_{\bar{n}}]}{\Gamma[\sigma-\rho-\tilde{\alpha}_{\bar{n}}-\gamma]}\\\
&\hskip
48.36967pt{}+\,\bar{\theta}(\sigma-\rho-\tilde{\alpha}_{\bar{n}}-1)\,\sum_{\gamma=0}^{\infty}\,\frac{(-1)^{\gamma}}{\gamma!}\,\frac{\Gamma[\gamma-\sigma+\rho+\tilde{\alpha}_{\bar{n}}+1]}{\Gamma[-\sigma+\rho+\tilde{\alpha}_{\bar{n}}+1]}\,\left(\frac{\hat{k}_{+}}{\tilde{k}^{\mathrm{max}}_{+}}\right)^{\gamma}\,\Bigg{\\}}\,.\end{split}$
(2.24)

Herein, we use the Heaviside theta function $\theta(x)$, with $\theta(x)=1$
for $x\geq 0$ and $0$ otherwise, and its counter part $\bar{\theta}(x)\equiv
1-\theta(x)$ to account for the expansion of the integrand. Even though within
the central rapidity region $e^{\pm Y_{H}}\sim\mathcal{O}(1)$ and thus
$\tilde{k}^{\mathrm{max}}_{+}\sim\mathcal{O}(m_{H}e^{\pm Y_{H}})$, the
integrand of Eq. (2.24) has already admitted a formally asymptotic series in
$({\hat{k}_{+}}/{\tilde{k}^{\mathrm{max}}_{+}})$, the participation of
$m_{\mathrm{T}}$ into the upper boundary condition can invoke inhomogeneous
behaviour after completing the phase space integral. To this end, a further
expansion is necessary with the aid of Eq. (2.11). It follows that

$\begin{split}\int_{0}^{(m_{H}-m_{\mathrm{T}})e^{-Y_{H}}}\,{\rm
d}\hat{k}_{+}\,\hat{k}_{+}^{\tilde{\alpha}_{n}+\gamma}\,=&\,\sum_{\gamma_{2}=0}^{\infty}\,\sum_{\gamma_{1}=0}^{1+\gamma+\tilde{\alpha}_{n}}\,\frac{(-1)^{\gamma_{1}}}{\gamma_{1}!\,\gamma_{2}!}\,\left(\frac{q_{\mathrm{T}}^{2}}{m_{H}^{2}}\right)^{\gamma_{2}}\,(e^{-Y_{H}}m_{H})^{1+\tilde{\alpha}_{n}+\gamma}\\\
&\hskip
56.9055pt\times\,\frac{\Gamma[1+\tfrac{\gamma_{1}}{2}]}{\Gamma[1-\gamma_{2}+\tfrac{\gamma_{1}}{2}]}\,\frac{\Gamma[\tilde{\alpha}_{n}+\gamma+1]}{\Gamma[\tilde{\alpha}_{n}+\gamma-\gamma_{1}+2]}\,.\end{split}$
(2.25)

Combining Eqs. (2.24)-(2.25) with Eq. (2.22), we then accomplish the expansion
of the $n$-ultra-collinear sector in the small parameter
$(q_{\mathrm{T}}/m_{H})$. An analogous procedure can be similarly applied onto
the ${\bar{n}}$-ultra-collinear case, with the exception of the conversion of
the integration variable $k_{+}$ to $k_{-}$ by means of the on-shell condition
$k_{+}k_{-}=q_{\mathrm{T}}^{2}$ and the substitution of $e^{+Y_{H}}$ for
$e^{-Y_{H}}$ in Eq. (2.25), as appropriate.

From those results, we observe that the boundary corrections here are always
associated with PDFs at the opposite end and therefore are a priori expected
to play an negligible role in comparison to
$\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$ in Eq.
(2.17). While this is expected to hold for the first few terms in the power
expansion, it is possible that the higher power corrections from
$\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$ become
of a similar size as the boundary corrections from
$\Delta\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$,
such that both parts of contributions are essential to reproduce the desired
asymptotic behaviour of the $q_{\mathrm{T}}$ spectrum. In Sec. 4, we will
deliver a quantitative assessment thereof.

###  2.4 Power corrections over the interior domain

####  2.4.1 Asymptotic expansion in momentum space

We can now perform the power expansion on the integral
$\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$ in the
low $q_{\mathrm{T}}$ area within the central rapidity region $e^{\pm
Y_{H}}\sim\mathcal{O}(1)$. As exhibited in Eq. (2.17),
$\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$
encompasses a variety of scales along the integration path, including the
ultra-soft scale $k_{+}\sim\mathcal{O}(q_{\mathrm{T}}^{2}/m_{H})$, the soft
one $k_{+}\sim\mathcal{O}(q_{\mathrm{T}})$, and the hard one
$k_{+}\sim\mathcal{O}(m_{H})$. Every one of those momentum modes in practice
prompts a distinct expansion of
$\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$. Hence,
in the spirit of the expansion by regions method [211, 212, 213, 214, 174] and
also according to the pattern of the resulting power series, we categorise the
phase space as follows 222It is interesting to note that the analogous
strategy in refining the dynamic regions has also been utilised previously in
the literature [215, 216, 105, 217, 218, 219] in probing the $k_{\mathrm{T}}$
factorisation in the small-$x$ area.,

$\begin{split}n\text{-collinear~{}mode}:\;\;\;&\nu_{n}<k_{+}\leq\tilde{k}^{\mathrm{max}}_{+};\,\vphantom{\frac{q_{\mathrm{T}}^{2}}{\nu_{{\bar{n}}}}}\\\
\text{transitional
range}:\;\;\;&\frac{q_{\mathrm{T}}^{2}}{\nu_{{\bar{n}}}}<k_{+}\leq\nu_{n};\,\\\
{\bar{n}}\text{-collinear~{}mode}:\;\;\;&\tilde{k}^{\mathrm{min}}_{+}\leq
k_{+}\leq\frac{q_{\mathrm{T}}^{2}}{\nu_{{\bar{n}}}}\;.\end{split}$ (2.26)

Here, a pair of auxiliary scales $\\{\nu_{n},\nu_{{\bar{n}}}\\}$ are
introduced to separate the $n$\- and ${\bar{n}}$-collinear regions from the
moderate range in between. Throughout this paper,
$\\{\nu_{n},\nu_{{\bar{n}}}\\}$ are chosen to be of similar magnitudes to the
hard scales $m_{H}e^{\pm Y_{H}}$ but always small enough that in the
intermediate domain of Eq. (2.26) the light-cone components $k_{+}$ and
$k_{-}$ are both well contained within the convergence radii of
${F}^{(\alpha_{n})}_{i/n,\beta_{n}}$ and
${F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}$, as stipulated in
Sec. 2.2, more specifically,

$\displaystyle m_{H}e^{-Y_{H}}\gtrsim\nu_{n}\sim\mathcal{O}(m_{H})\gg
q_{\mathrm{T}}\,,\quad\quad
m_{H}e^{Y_{H}}\gtrsim\nu_{{\bar{n}}}\sim\mathcal{O}(m_{H})\gg
q_{\mathrm{T}}\,.$ (2.27)

Alternative choices of those two auxiliary boundaries can be used to examine
the asymptotic behaviour of
$\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$. For
example, the scales proposed in [219] may impact the expressions of the
individual sectors defined in Eq. (2.26) but will leave the resulting power
series unchanged upon summation of all ingredients. Having chosen
$\\{\nu_{n},\nu_{{\bar{n}}}\\}$ according to Eq. (2.27), the momentum
$k^{\mu}$ obeys the unambiguous scaling law $k_{+}\sim\mathcal{O}(m_{H})\gg
k_{-}\sim\mathcal{O}(q_{\mathrm{T}}^{2}/m_{H})$ in the $n$-collinear mode,
from which we can expand
${F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}$ around the intrinsic
scale $m_{H}e^{+Y_{H}}$,

$\displaystyle\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{c}\,=\,\sum_{\omega=\rho}^{\infty}\,(q_{\mathrm{T}}^{2})^{\omega}\,\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{c}\,,$
(2.28)

where

$\displaystyle\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{c}\,\equiv\,\int^{\tilde{k}^{\mathrm{max}}_{+}}_{\nu_{n}}\,\frac{{\rm
d}k_{+}}{k_{+}}\,\frac{(k_{+})^{\sigma-\omega}}{(\omega-\rho)!}\,{F}^{(\alpha_{n})}_{i/n,\beta_{n}}\left(k_{+}+m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}}+\omega-\rho)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)\,.$
(2.29)

Here, all $q_{\mathrm{T}}$ dependences have been extracted into the prefactor
of Eq. (2.28) to express the manifest power accuracy of the contribution.
$\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}$
collects all hard scales and thus is always of $\mathcal{O}(1)$. It merits
noting that in the central rapidity region $e^{\pm Y_{H}}\sim\mathcal{O}(1)$,
the lower boundary $\nu_{n}$ in Eq. (2.29) is well separated from the upper
one and thus there is no need to utilise a further expansion here, in contrast
to Eq. (2.22).

The same method can also be applied to the ${\bar{n}}$-collinear mode, after
transforming the integration variable from $k_{+}$ to $k_{-}$ through the on-
shell condition $k_{+}k_{-}=q_{\mathrm{T}}^{2}$. It evaluates to,

$\displaystyle\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}}\,=\,\sum_{\omega=\sigma}^{\infty}\,(q_{\mathrm{T}}^{2})^{\omega}\,\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}}\,,$
(2.30)

with

$\displaystyle\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}}\,\equiv\,\int^{\frac{q_{\mathrm{T}}^{2}}{\tilde{k}^{\mathrm{min}}_{+}}}_{\nu_{{\bar{n}}}}\,\frac{{\rm
d}k_{-}}{k_{-}}\,\frac{(k_{-})^{\rho-\omega}}{(\omega-\sigma)!}\,{F}^{(\alpha_{n}+\omega-\sigma)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}\left(k_{-}+m_{H}e^{+Y_{H}}\right)\,.$
(2.31)

This leaves the transitional domain, as defined in Eq. (2.26), to be
calculated. In light of the scaling rules of Eq. (2.27), we can expand both
${F}^{(\alpha_{n})}_{i/n,\beta_{n}}$ and
${F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}$ around the scales
$m_{H}e^{\pm Y_{H}}$, more specifically,

$\begin{split}\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{t}\,=\,&\sum_{\lambda,\eta=0}^{\infty}\,\int^{\nu_{n}}_{\frac{q_{\mathrm{T}}^{2}}{\nu_{{\bar{n}}}}}\,\frac{{\rm
d}k_{+}}{k_{+}}\,\frac{(k_{-})^{\rho+\lambda}\,(k_{+})^{\sigma+\eta}}{\lambda!\,\eta!}\,{F}^{(\alpha_{n}+\eta)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}}+\lambda)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)\,\\\
=\,&\sum_{\lambda,\eta=0}^{\infty}\,\frac{{F}^{(\alpha_{n}+\eta)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}}+\lambda)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)}{\lambda!\,\eta!}\,\\\
&\hskip
22.76219pt\times\,\Bigg{\\{}\underbrace{\frac{(q_{\mathrm{T}}^{2})^{\rho+\lambda}\,\nu_{n}^{\sigma+\eta-\rho-\lambda}}{\sigma+\eta-\rho-\lambda}\,\bar{\delta}^{\sigma+\eta}_{\rho+\lambda}}_{n\text{-col}}+\underbrace{\frac{(q_{\mathrm{T}}^{2})^{\sigma+\eta}\,\nu_{{\bar{n}}}^{\rho+\lambda-\sigma-\eta}}{\rho+\lambda-\sigma-\eta}\,\bar{\delta}^{\sigma+\eta}_{\rho+\lambda}}_{{\bar{n}}\text{-col}}+\underbrace{(q_{\mathrm{T}}^{2})^{\sigma+\eta}\,\ln\left[\frac{\nu_{n}\nu_{{\bar{n}}}}{q_{\mathrm{T}}^{2}}\right]\,{\delta}^{\sigma+\eta}_{\rho+\lambda}}_{n\text{-col}\,,~{}{\bar{n}}\text{-col}\,,~{}\text{s}}\,\Bigg{\\}}\,,\end{split}$
(2.32)

where we have already grouped the expression according to their powers in
$\nu_{n}$ and $\nu_{{\bar{n}}}$. $\delta^{\alpha}_{\beta}$ is the Kronecker
Delta, $\delta^{\alpha}_{\beta}=1$ if $\alpha=\beta$ and zero otherwise,
$\bar{\delta}^{\alpha}_{\beta}\equiv 1-\delta^{\alpha}_{\beta}$ is its
complement. We note that even for fixed values of $\lambda$ and $\eta$, the
phase space integral results in terms with a variety of power accuracies. For
instance, in the first term on r.h.s. the power precision is determined by the
exponent of $q_{\mathrm{T}}$, which is identical to that of $k_{-}$ in the
integrand. Hence, the scaling of the emitted momentum $k^{\mu}$ here is
determined by the $n$-collinear mode
$k_{-}\sim\mathcal{O}(q_{\mathrm{T}}^{2}/m_{H})$ and
$k_{+}\sim\mathcal{O}(m_{H})$. Similarly, applying this scaling analysis to
the second term, its power precision is given by the exponent of $k_{+}$ and
the ${\bar{n}}$-collinear mode with $k_{-}\sim\mathcal{O}(m_{H})$ and
$k_{+}\sim\mathcal{O}(q_{\mathrm{T}}^{2}/m_{H})$. The interpretation of the
scaling of the integration variable of the last term of Eq. (2.32) is more
flexible due to the presence of the constraint $\rho+\lambda=\sigma+\eta$.
Hence, the emitted parton can be constituted of either the
$n({\bar{n}})$-collinear particle or the soft one $k_{-}\sim
k_{+}\sim\mathcal{O}(q_{\mathrm{T}})$. This inhomogeneity stems from our
choice of the momentum cutoffs $\nu_{n}$ and $\nu_{{\bar{n}}}$ in Eq. (2.27),
from which multiple dynamic modes are allowed to participate in the
transitional area, such as the $n$-collinear one, ${\bar{n}}$-collinear one,
or the soft one. In principle, reducing $\nu_{n}$ and $\nu_{{\bar{n}}}$ down
to $\mathcal{O}(q_{\mathrm{T}})$ will filter out the collinear modes in the
transitional domain rendering it homogeneous again, at the price of additional
soft scales in both collinear sectors through the boundary condition,
necessitating a second expansion in the small parameters $\nu_{n}/m_{H}$ and
$\nu_{\bar{n}}/m_{H}$ within there. It remains to be noted, while alternative
choices of $\nu_{n}$ and $\nu_{{\bar{n}}}$ may modify the dynamics in each
sector, their combined power series is invariant under this choice. Hence, in
this work we choose the prescription of Eq. (2.27), i.e. homogenous collinear
sectors, for simplicity.

By collecting and refactoring terms of a common power in $q_{\mathrm{T}}$ we
can recast Eq. (2.32) as

$\displaystyle\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{t}\,=\,\left(\sum_{\omega=\rho}^{\infty}\,(q_{\mathrm{T}}^{2})^{\omega}\,\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{cs}\,\right)+\left(\sum_{\omega=\sigma}^{\infty}\,(q_{\mathrm{T}}^{2})^{\omega}\,\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}s}\,\right)\,,$
(2.33)

where

$\begin{split}&\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{cs}=\sum_{\eta=0}^{\infty}\frac{{F}^{(\eta+\alpha_{n})}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right){F}^{(\alpha_{\bar{n}}+\omega-\rho)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)\,}{(\omega-\rho)!\,\eta!}\Bigg{\\{}\frac{\nu_{n}^{\sigma+\eta-\omega}}{\sigma+\eta-\omega}\bar{\delta}^{\sigma+\eta}_{\omega}+\ln\left[\frac{\nu_{n}}{q_{\mathrm{T}}}\right]{\delta}^{\sigma+\eta}_{\omega}\,\Bigg{\\}}\,,\\\
&\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}s}=\sum_{\lambda=0}^{\infty}\frac{{F}^{(\alpha_{n}+\omega-\sigma)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right){F}^{(\alpha_{\bar{n}}+\lambda)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)}{(\omega-\sigma)!\,\lambda!}\Bigg{\\{}\frac{\nu_{{\bar{n}}}^{\rho+\lambda-\omega}}{\rho+\lambda-\omega}\bar{\delta}^{\rho+\lambda}_{\omega}+\ln\left[\frac{\nu_{{\bar{n}}}}{q_{\mathrm{T}}}\right]{\delta}^{\rho+\lambda}_{\omega}\Bigg{\\}}\,.\end{split}$
(2.34)

Here we have re-organized the logarithmic terms in Eq. (2.32) into the
$n$-collinear-soft ($cs$) and ${\bar{n}}$-collinear-soft ($\bar{c}s$) sectors
based on their scale dependences. At variance from the results in Eq. (2.29)
and Eq. (2.31) in the $n$\- and ${\bar{n}}$-collinear sectors, where the
coefficients at each power consist of only the hard scales $k_{\pm}$,
respectively, Eq. (2.34) here introduces an additional dependence on
$\ln\left({\nu_{n({\bar{n}})}}/{q_{\mathrm{T}}}\right)$. In order to
facilitate numerical calculations of these coefficients, we rewrite the
infinite series in Eq. (2.34) in terms of the integrals over $k_{\pm}$, in
analogy to Eqs. (2.29) and (2.31). If $\omega<\sigma$ in Eq. (2.34), this
transformation is immediate,

$\begin{split}&\bar{\theta}(\omega-\sigma)\,\sum_{\eta=0}^{\infty}\,\frac{{F}^{(\eta+\alpha_{n})}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}}+\omega-\rho)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)\,}{(\omega-\rho)!\,\eta!}\,\frac{\nu_{n}^{\sigma+\eta-\omega}}{\sigma+\eta-\omega}\,\\\
&\qquad\qquad\qquad\Longrightarrow\,\frac{\bar{\theta}(\omega-\sigma)}{(\omega-\rho)!}\,\int^{\nu_{n}}_{0}\,\frac{{\rm
d}k_{+}}{k_{+}}\,k_{+}^{\sigma-\omega}\,{F}^{(\alpha_{n})}_{i/n,\beta_{n}}\left(k_{+}+m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}}+\omega-\rho)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)\,.\end{split}$
(2.35)

Here we have removed $\bar{\delta}^{\sigma+\eta}_{\omega}$ from the
expressions since it always equals unity while $\omega<\sigma$ and $\eta\geq
0$. The same reasoning can also be used to neglect the logarithmic
contributions in Eq. (2.34). However, it is not straightforward to apply Eq.
(2.35) in the remaining space with $\omega\geq\sigma$ as the integrand is
singular in the limit $k_{+}\to 0$. To this end, we make use of the higher-
ranked star distribution,

$\displaystyle\int_{0}^{\Lambda}\,{\rm
d}x\,\left[\frac{1}{x^{m}}\right]^{\nu}_{*}\,f(x)\equiv\int_{\nu}^{\Lambda}\,{\rm
d}x\,\frac{f(x)}{x^{m}}+\int_{0}^{\nu}\,{\rm
d}x\,\frac{1}{x^{m}}\,\left[f(x)-\sum_{n=0}^{m-1}\frac{x^{n}}{n!}f^{(n)}(0)\right]\,,$
(2.36)

to maintain the end point analyticity. Therein, $f(x)$ is assumed to be always
differentiable at $x=0$. With $m=1$, Eq. (2.36) reduces to the customary star
distribution proposed in [220, 221]. Equipped with eq. (2.36), we are now able
to convert the $n$-collinear-soft contribution in the case $\omega\geq\sigma$
and $\eta>\omega-\sigma$,

$\begin{split}&{\theta}(\omega-\sigma)\,\sum_{\eta=\omega-\sigma+1}^{\infty}\,\frac{{F}^{(\eta+\alpha_{n})}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}}+\omega-\rho)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)\,}{(\omega-\rho)!\,\eta!}\,\frac{\nu_{n}^{\sigma+\eta-\omega}}{\sigma+\eta-\omega}\,\\\
&\qquad\qquad\Longrightarrow\,\frac{{\theta}(\omega-\sigma)}{(\omega-\rho)!}\,\int^{\nu_{n}}_{0}\,{\rm
d}k_{+}\,\left[\frac{1}{k_{+}^{\omega-\sigma+1}}\right]^{\nu_{n}}_{*}\,{F}^{(\alpha_{n})}_{i/n,\beta_{n}}\left(k_{+}+m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}}+\omega-\rho)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)\,.\end{split}$
(2.37)

The remaining contributions in Eq. (2.34) include the logarithms
$\ln\left({\nu_{n}}/{q_{\mathrm{T}}}\right)$ as well as the first
$(\omega-\sigma)$ terms of the Taylor series. We will refrain from
transforming those contributions here since, at a given power precision, the
number of remaining terms are always finite.

Combining above results we can recast the $n$-collinear-soft contribution of
Eq. (2.34) as follows

$\begin{split}\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{cs}=&\phantom{{}+}\frac{\bar{\theta}(\omega-\sigma)}{(\omega-\rho)!}\,\int^{\nu_{n}}_{0}\,\frac{{\rm
d}k_{+}}{k_{+}}\,k_{+}^{\sigma-\omega}\,{F}^{(\alpha_{n})}_{i/n,\beta_{n}}\left(k_{+}+m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}}+\omega-\rho)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)\,\\\
&{}+\frac{{\theta}(\omega-\sigma)}{(\omega-\rho)!}\,\int^{\nu_{n}}_{0}\,{\rm
d}k_{+}\,\left[\frac{1}{k_{+}^{\omega-\sigma+1}}\right]^{\nu_{n}}_{*}\,{F}^{(\alpha_{n})}_{i/n,\beta_{n}}\left(k_{+}+m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}}+\omega-\rho)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)\,\\\
&{}+\frac{{\theta}(\omega-\sigma-1)}{(\omega-\rho)!}\,\sum_{\eta=0}^{\omega-\sigma-1}\,\frac{{F}^{(\alpha_{n}+\eta)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}}+\omega-\rho)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)\,}{\eta!}\,\frac{\nu_{n}^{\sigma+\eta-\omega}}{\sigma+\eta-\omega}\,\\\
&{}+\frac{{\theta}(\omega-\sigma)}{(\omega-\rho)!}\,\ln\left[\frac{\nu_{n}}{q_{\mathrm{T}}}\right]\,\frac{{F}^{(\alpha_{n}+\omega-\sigma)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}}+\omega-\rho)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)\,}{(\omega-\sigma)!}\,.\end{split}$
(2.38)

By analogy, the contribution of the ${\bar{n}}$-collinear-soft sector in Eq.
(2.34) can be transformed to finite integrals following an analogous path,
giving

$\begin{split}\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}s}=&\phantom{{}+}\frac{\bar{\theta}(\omega-\rho)}{(\omega-\sigma)!}\,\int^{\nu_{{\bar{n}}}}_{0}\,\frac{{\rm
d}k_{-}}{k_{-}}\,k_{-}^{\rho-\omega}\,{F}^{(\alpha_{n}+\omega-\sigma)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}\left(k_{-}+m_{H}e^{+Y_{H}}\right)\,\\\
&{}+\frac{{\theta}(\omega-\rho)}{(\omega-\sigma)!}\,\int^{\nu_{{\bar{n}}}}_{0}\,{\rm
d}k_{-}\,\left[\frac{1}{k_{-}^{\omega-\rho+1}}\right]^{\nu_{{\bar{n}}}}_{*}\,{F}^{(\alpha_{n}+\omega-\sigma)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}\left(k_{-}+m_{H}e^{+Y_{H}}\right)\,\\\
&{}+\frac{{\theta}(\omega-\rho-1)}{(\omega-\sigma)!}\,\sum_{\lambda=0}^{\omega-\rho-1}\,\frac{{F}^{(\alpha_{n}+\omega-\sigma)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}}+\lambda)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)\,}{\lambda!}\,\frac{\nu_{{\bar{n}}}^{\rho+\lambda-\omega}}{\rho+\lambda-\omega}\,\\\
&{}+\frac{{\theta}(\omega-\rho)}{(\omega-\sigma)!}\,\ln\left[\frac{\nu_{{\bar{n}}}}{q_{\mathrm{T}}}\right]\,\frac{{F}^{(\alpha_{n}+\omega-\sigma)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}}+\omega-\rho)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)\,}{(\omega-\rho)!}\,.\end{split}$
(2.39)

Using the above result we can now combine the $n$-collinear sector of Eq.
(2.28), the ${\bar{n}}$-collinear sector of Eq. (2.30), and the transitional
sector of Eqs. (2.38) and (2.39), arriving at the expansion of
$\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$ of Eq.
(2.17) in powers of $q_{\mathrm{T}}^{2}$ in the small $q_{\mathrm{T}}$ domain,

$\displaystyle\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}=\sum_{\omega}(q_{\mathrm{T}}^{2})^{\omega}\,\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\equiv\sum_{\omega}(q_{\mathrm{T}}^{2})^{\omega}\Bigg{\\{}\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{c}+\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}}+\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{cs}+\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}s}\Bigg{\\}}\,.$
(2.40)

At this point it merits reminding that the pair of auxiliary scales, $\nu_{n}$
and $\nu_{{\bar{n}}}$, that we have employed to define and separate the
dynamic regions in Eq. (2.26) appear explicitly in the expressions of all
three individual domains. For consistency’s sake, it is of the essence to
examine whether the combination of all domains also exhibits such a
dependence, or whether this dependence, in fact, vanishes after summing over
all the ingredients. To this end, we take the derivatives of the r.h.s. of Eq.
(2.40) with respect to $\nu_{n}$, and it evaluates to

$\begin{split}&\frac{\partial}{\partial\nu_{n}}\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{c}=-\,\frac{\partial}{\partial\nu_{n}}\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{cs}=-\frac{(\nu_{n})^{\sigma-\omega-1}}{(\omega-\rho)}{F}^{(\alpha_{n})}_{i/n,\beta_{n}}\left(\nu_{n}+m_{H}e^{-Y_{H}}\right){F}^{(\alpha_{\bar{n}}+\omega-\rho)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right),\\\
&\frac{\partial}{\partial\nu_{n}}\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}}=\frac{\partial}{\partial\nu_{n}}\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}s}=0\,,\end{split}$
(2.41)

where the results of Eqs. (2.29) and (2.38) are utilised to derive the r.h.s.
of Eq. (2.41). In addition, the terms arising from the derivative of the star
distribution, defined in Eq. (2.36), have been cancelled against the $\nu_{n}$
dependences arising from the logarithmic term and the truncated Taylor
polynomials in Eq. (2.38). A similar result can also be found for the
$\nu_{\bar{n}}$ dependence of Eq. (2.40), giving

$\begin{split}&\frac{\partial}{\partial\nu_{\bar{n}}}\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{c}=\frac{\partial}{\partial\nu_{\bar{n}}}\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{cs}=0\,,\\\
&\frac{\partial}{\partial\nu_{\bar{n}}}\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}}=-\frac{\partial}{\partial\nu_{\bar{n}}}\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}s}=-\frac{(\nu_{\bar{n}})^{\rho-\omega-1}}{(\omega-\sigma)!}{F}^{(\alpha_{n}+\omega-\sigma)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right){F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}\left(\nu_{\bar{n}}+m_{H}e^{+Y_{H}}\right).\end{split}$
(2.42)

Combining the results of Eqs. (2.41) and (2.42), respectively, we can conclude
that the coefficient
$\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}$
in Eq. (2.40) is indeed independent of the choice of the separators $\nu_{n}$
and $\nu_{\bar{n}}$ at each power. This agrees with the expectation from a
direct expansion of
$\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$,
currently beyond our capabilities as the PDFs are not rigorously known
analytically, that it should be a function of $q_{\mathrm{T}}$, $s$, and
$m_{H}$ only and independent of $\nu_{n}$ or $\nu_{\bar{n}}$.

####  2.4.2 Rapidity regularisation and zero-bin subtraction

In the previous section, with the help of the auxiliary scales $\nu_{n}$ and
$\nu_{\bar{n}}$, we have derived the power series for
$\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$ in Eq.
(2.40) in the vicinity $q_{\mathrm{T}}=0$ at central rapidities, $e^{\pm
Y_{H}}\sim\mathcal{O}(1)$. Even though this strategy has accomplished an
asymptotic expansion of
$\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$, and in
turn the $q_{\mathrm{T}}$ distribution, at NLO, it is not straightforward to
generalise this method to higher perturbative orders. The presence of
subdivisions of phase space through the auxiliary scales $\nu_{n}$ and
$\nu_{\bar{n}}$, generally inducing additional scales to phase space and loop
integrals, complicates the calculations substantially. Further, the absence of
additional constraints on the integration paths is also one of prerequisites
to establish the factorisation in a SCET-based analysis. To this end, we will
discuss the rearrangement of the ingredients defined in Eqs. (2.29) and (2.31)
as well as Eqs. (2.38) and (2.39) such that the dependence on the auxiliary
scales to subdivide the phase space is removed. In particular, recalling our
findings regarding the $\nu_{n({\bar{n}})}$-independence of
$\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$, but not
its components in Eqs. (2.41) and (2.42), such a rearrangement is warranted.
Although in the discussion below, we will still restrict ourselves to the NLO
constituent
$\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$
extracted from Eq. (2.3) as an illustrative example, it is expected that the
conclusion here can provide the theoretical baseline for an analysis in more
general situations.

We start by examining the collinear sectors. Here, eliminating the dependences
of the auxiliary scales $\nu_{n}$ and $\nu_{\bar{n}}$ amounts to reducing the
lower boundaries in Eq. (2.29) and Eq. (2.31) from their intrinsic domains
defined in Eq. (2.27) to the origin. To be precise, we split

$\displaystyle\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{c}\,=\,\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{c}\,-\,\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{c0}\,,\qquad\quad\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}}\,=\,\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}}\,-\,\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}0}\,,$
(2.43)

where
$\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}_{c}$
and
$\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}_{\bar{c}}$
collect the collinear contributions, integrating the light-cone momentum over
the entire range. Their expressions read

$\begin{split}\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{c}&\equiv\lim_{\tau\to
0}\int^{\tilde{k}^{\mathrm{max}}_{+}}_{0}\frac{{\rm
d}k_{+}}{k_{+}}\mathcal{R}(k_{-},k_{+},\tau)\frac{(k_{+})^{\sigma-\omega}}{(\omega-\rho)!}{F}^{(\alpha_{n})}_{i/n,\beta_{n}}\left(k_{+}+m_{H}e^{-Y_{H}}\right){F}^{(\alpha_{\bar{n}}+\omega-\rho)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right),\\\
\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}}&\equiv\lim_{\tau\to
0}\int^{\frac{q_{\mathrm{T}}^{2}}{\tilde{k}^{\mathrm{min}}_{+}}}_{0}\frac{{\rm
d}k_{-}}{k_{-}}\mathcal{R}(k_{-},k_{+},\tau)\frac{(k_{-})^{\rho-\omega}}{(\omega-\sigma)!}{F}^{(\alpha_{n}+\omega-\sigma)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right){F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}\left(k_{-}+m_{H}e^{+Y_{H}}\right).\end{split}$
(2.44)

Therein, in order to regularise the rapidity divergence in the limit
$k_{\pm}\to 0$, we multiply the original integrand by
$\mathcal{R}(k_{-},k_{+},\tau)$, which is a function of the light-cone momenta
$k_{\pm}$ and a real parameter $\tau$, satisfying,

$\displaystyle\lim_{\tau\to
0}\,\mathcal{R}(k_{-},k_{+},\tau)\,=\,\left\\{\begin{matrix}1&\text{$k_{\pm}\neq
0$}\,,\\\ 0&\text{$k_{+}=0$~{}or~{}$k_{-}=0$}\,.\end{matrix}\right.$ (2.45)

For the regions $k_{\pm}\neq 0$, the integrands in Eq. (2.44) are uniformly
convergent, which permits reducing the rapidity regulator
$\mathcal{R}(k_{-},k_{+},\tau)$ in the limit $\tau\to 0$ and restoring the
integrands to their original forms. Nevertheless, in the asymptotic area
$k_{+}=0$ or $k_{-}=0$, the role of the rapidity regulator becomes
indispensable, ensuring that the phase space integrals in Eq. (2.44) are
always well defined. The discussion in this section will emphasise rapidity
regulators that are able to preserve their expressions in both the
$n$-collinear and ${\bar{n}}$-collinear sectors for all powers, such as the
analytic regulator of [74, 76], the exponential regulator of [79, 80], and the
pure rapidity regulator of [113, 175]. We will dub them the conservative
rapidity regulator or conservative regularisation scheme (CRa) hereafter.
Other than that, there are also alternative proposals, including the
$\Delta$-regulator [162, 165, 169, 168, 176] and the $\eta$-regulator [77,
78]. At this point, either the rapidity regulators themselves or the
regulator-dressed propagators will participate in the power expansion such
that in practice the rapidity divergences can be regularised in varying
strategies in different ingredients from power to power. They will be called
dissipative rapidity regulator (DRa) in this work. We will postpone their
investigation for now but discuss the possible prescriptions in App. A.

In Eq. (2.43), we have also introduced
$\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}_{c0}$
and
$\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}_{\bar{c}0}$
in order to restore the equality after. Their expressions read,

$\begin{split}&\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{c0}=\lim_{\tau\to
0}\int^{\nu_{n}}_{0}\frac{{\rm
d}k_{+}}{k_{+}}\mathcal{R}(k_{-},k_{+},\tau)\frac{(k_{+})^{\sigma-\omega}}{(\omega-\rho)!}{F}^{(\alpha_{n})}_{i/n,\beta_{n}}\left(k_{+}+m_{H}e^{-Y_{H}}\right){F}^{(\alpha_{\bar{n}}+\omega-\rho)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right),\\\
&\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}0}=\lim_{\tau\to
0}\int^{\nu_{{\bar{n}}}}_{0}\frac{{\rm
d}k_{-}}{k_{-}}\mathcal{R}(k_{-},k_{+},\tau)\frac{(k_{-})^{\rho-\omega}}{(\omega-\sigma)!}{F}^{(\alpha_{n}+\omega-\sigma)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right){F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}\left(k_{-}+m_{H}e^{+Y_{H}}\right).\end{split}$
(2.46)

Since these functions are always associated with the $k_{\pm}$ integrals from
the origin point up to the cutoffs $\nu_{n({\bar{n}})}$, we will make use of
the subscripts $``c0"$ $(``\bar{c}0")$ to represent them hereafter. From Eq.
(2.46), noting that the integration variables $k_{\pm}$ are both well situated
within the convergence radii of ${F}_{i/n,\beta_{n}}$ and
${F}_{j/{\bar{n}},\beta_{\bar{n}}}$, according to Sec. 2.2, we can now perform
the power expansion in $k_{\pm}$. We then have,

$\begin{split}&\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{c0}\,=\,\frac{{F}^{(\alpha_{\bar{n}}+\omega-\rho)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)}{(\omega-\rho)!}\,\left\\{\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{c0r}\,+\,\theta\left(\omega-\sigma\right)\,\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{c0d}\,\right\\}\,,\\\
&\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}0}\,=\,\frac{{F}^{(\alpha_{n}+\omega-\sigma)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)}{(\omega-\sigma)!}\,\left\\{\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}0r}\,+\,\theta\left(\omega-\rho\right)\,\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}0d}\,\right\\}\,,\end{split}$
(2.47)

where

$\begin{split}\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{c0r}&=\,\sum^{\infty}_{\eta=\mathrm{max}\\{0,(\omega-\sigma+1)\\}}\frac{\nu^{\sigma-\omega+\eta}_{n}}{\sigma-\omega+\eta}\frac{{F}^{(\alpha_{n}+\eta)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)}{\eta!}\,,\\\
\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}0r}&=\,\sum^{\infty}_{\lambda=\mathrm{max}\\{0,(\omega-\rho+1)\\}}\frac{\nu^{\rho-\omega+\lambda}_{\bar{n}}}{\rho-\omega+\lambda}\frac{{F}^{(\alpha_{\bar{n}}+\lambda)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)}{\lambda!}\,,\\\
\end{split}$ (2.48)

are the regular components, where the regulator could be safely removed by the
means of Eq. (2.45), and

$\begin{split}\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{c0d}&=\,\sum_{\eta=0}^{\omega-\sigma}\,\frac{{F}^{(\alpha_{n}+\eta)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)}{\eta!}\,\left\\{\lim_{\tau\to
0}\int^{\nu_{n}}_{0}\frac{{\rm
d}k_{+}}{k_{+}}\mathcal{R}(k_{-},k_{+},\tau)(k_{+})^{\sigma-\omega+\eta}\right\\}\,,\\\
\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}0d}&=\,\sum_{\lambda=0}^{\omega-\rho}\,\frac{{F}^{(\alpha_{\bar{n}}+\lambda)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)}{\lambda!}\,\left\\{\lim_{\tau\to
0}\int^{\nu_{\bar{n}}}_{0}\frac{{\rm
d}k_{-}}{k_{-}}\mathcal{R}(k_{-},k_{+},\tau)(k_{-})^{\rho-\omega+\lambda}\right\\}\,.\end{split}$
(2.49)

contain the (regulated) rapidity divergences at $k_{\pm}\to 0$. Here, the
regulator $\mathcal{R}(k_{-},k_{+},\tau)$ is essential and the respective
integrals are therefore dependent on the choice of the regularisation scheme.
Substituting Eq. (2.48) into Eq. (2.47), we observe that the results for
$\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}_{c0}$
and
$\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}_{c0}$
still explicitly depend on the auxiliary scales $\nu_{n({\bar{n}})}$. However,
according to Eq. (2.41), all those dependences will drop out upon combining
them with the contributions from the transitional region of Eq. (2.34), and it
follows that

$\displaystyle\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{cs}+\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}s}-\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{c0}-\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}0}=\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{comb},$
(2.50)

where on the r.h.s.,
$\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}_{comb}$,
comprises a set of one-fold integrals of $k_{+}$ or $k_{-}$ over the range
$[0,+\infty]$. Its precise form will be specified below.

In the literature, there are two preferences for expressing this combined
result. The first one, used extensively when using the method of expansion by
regions [211, 212, 213, 214, 174] and SCET [222, 223, 224, 225, 226, 227, 228,
229, 230, 231], comprises the soft interaction
$k_{\pm}\sim\mathcal{O}(q_{\mathrm{T}})$ in full as well as the subtraction
terms to remove the overlap with the collinear regions. Hereafter, we will
refer to the results from this formulation as the full soft
$\langle\mathrm{FS}\rangle$ prescription. It yields that

$\displaystyle\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{comb}\,=\,\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}^{\langle\mathrm{FS}\rangle}_{s}\,-\,\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}^{\langle\mathrm{FS}\rangle}_{c\bar{c}s}\,.$
(2.51)

Here
$\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}^{\langle\mathrm{FS}\rangle}_{s}$
stands for the soft contribution, derived by expanding the integrand in Eq.
(2.17) in accordance with the soft scaling
$k_{\pm}\sim\mathcal{O}(q_{\mathrm{T}})$. At the $\omega$th-power, the result
reads,

$\displaystyle\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}^{\langle\mathrm{FS}\rangle}_{s}\,=$
$\displaystyle\,\theta(\omega-\rho-\sigma)\,\sum_{\eta=0}^{(\omega-\rho-\sigma)}\,{F}^{(\alpha_{n}+\omega-\rho-\sigma-\eta)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}}+\eta)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)\,\lim_{\tau\to
0}\,\int^{\infty}_{0}\,\frac{{\rm d}k_{+}}{k_{+}}\,$
$\displaystyle\times\mathcal{R}(k_{-},k_{+},\tau)\,\frac{(k_{-})^{\rho+\eta-\omega}}{\eta!}\,\frac{(k_{+})^{-\rho-\eta}}{(\omega-\rho-\sigma-\eta)!}\,,$
(2.52)

where we have extracted the factor of $(q_{\mathrm{T}}^{2})^{\omega}$ from the
integrals in Eq. (2.4.2), in line with the definition of the power series in
Eq. (2.40). Given the expression in Eq. (2.4.2), we can obtain an expression
for
$\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}^{\langle\mathrm{FS}\rangle}_{c\bar{c}s}$
by comparing Eq. (2.51) with Eq. (2.50). As this subtrahend at LP is
intimately associated with the soft limit of the collinear functions and thus
the zeroth-bin of the label momentum within the position-momentum space hybrid
representation,
$\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}^{\langle\mathrm{FS}\rangle}_{c\bar{c}s}$
is frequently dubbed the zero-bin contribution [159, 160, 161, 162] or the
soft subtraction term [163, 164, 160, 161, 165].

On the other hand, there are also alternative proposals in the literature
[106, 111, 166, 167, 114]. There, a decomposition of Eq. (2.50) without an
explicit soft sector is chosen. The virtue of such a scheme is that in absence
of the soft contribution the soft-collinear interaction vertices do not need
to be expanded using multiple scaling regimes, and the factorisation of the
distinct dynamic modes is thus considerably simplified. In this work, we will
dub this approach the one in the non-soft $\langle\mathrm{NS}\rangle$
prescription. It follows that,

$\displaystyle\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{comb}\,=\,-\,\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}^{\langle\mathrm{NS}\rangle}_{c\bar{c}}\,.$
(2.53)

Hence, the interior contribution
$\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}$
in Eq. (2.17) consists of only the collinear sectors in Eq. (2.44) as well as
$\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}^{\langle\mathrm{NS}\rangle}_{c\bar{c}}$.
Since
$\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}^{\langle\mathrm{NS}\rangle}_{c\bar{c}}$
here effectively removes the overlap between the $n$-collinear and
${\bar{n}}$-collinear contributions, it is also named the overlap subtraction
term in [166, 167, 114]. Notwithstanding, in view of the functional similarity
of
$\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}^{\langle\mathrm{NS}\rangle}_{c\bar{c}}$
and
$\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}^{\langle\mathrm{FS}\rangle}_{c\bar{c}s}$,
we will not further distinguish the terminologies in this paper, such as the
overlapping subtraction, the zero-bin contribution, or the soft remover.
Instead, they will be all regarded to be conceptually equivalent but evaluated
within different prescriptions.

Comparing Eq. (2.53) with Eq. (2.51), we note that the subtraction terms in
those two schemes possess the following relation,

$\displaystyle\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}^{\langle\mathrm{NS}\rangle}_{c\bar{c}}\,=\,\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}^{\langle\mathrm{FS}\rangle}_{c\bar{c}s}\,-\,\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}^{\langle\mathrm{FS}\rangle}_{s}\,.$
(2.54)

In the following, we will be focused on the derivation of
$\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}^{\langle\mathrm{NS}\rangle}_{c\bar{c}}$.
The expressions within the $\langle\mathrm{FS}\rangle$ prescription can be
obtained by means of the relationship above.

We begin with the transformation of
$\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}_{cs}$
and
$\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}_{\bar{c}s}$.
As illustrated in Eq. (2.34),
$\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}_{cs}$
comprises the logarithmic terms in the case of ${\omega}={(\sigma+\eta)}$ and
powers of the auxiliary scales $\nu_{n({\bar{n}})}$ otherwise. The latter
contribution can be further categorised based on the relation between $\omega$
and $(\sigma+\eta)$, more specifically,

$\begin{split}&\sum_{\eta=0}^{\infty}\frac{{F}^{(\eta+\alpha_{n})}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right){F}^{(\alpha_{\bar{n}}+\omega-\rho)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)\,}{(\omega-\rho)!\,\eta!}\frac{\nu_{n}^{\sigma+\eta-\omega}}{\sigma+\eta-\omega}\bar{\delta}^{\sigma+\eta}_{\omega}\\\
=\,&\frac{{F}^{(\alpha_{\bar{n}}+\omega-\rho)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)}{(\omega-\rho)!}\Bigg{\\{}\sum_{\eta=\mathrm{max}\\{0,(\omega-\sigma+1)\\}}^{\infty}\frac{{F}^{(\alpha_{n}+\eta)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)}{\eta!}\frac{\nu_{n}^{\sigma+\eta-\omega}}{\sigma+\eta-\omega}\\\
&\hskip
85.35826pt{}-\theta(\omega-\sigma-1)\sum_{\eta=0}^{\omega-\sigma-1}\frac{{F}^{(\alpha_{n}+\eta)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)}{\eta!\,}\lim_{\tau\to
0}\int^{\infty}_{\nu_{n}}\frac{{\rm
d}k_{+}}{k_{+}}\mathcal{R}(k_{-},k_{+},\tau)(k_{+})^{\sigma-\omega+\eta}\Bigg{\\}}\,.\end{split}$
(2.55)

Therein, within the curly brackets, the first term collects all the
contributions with $\omega<(\sigma+\eta)$, whilst the second term collects all
contributions with $\omega>(\sigma+\eta)$ and has been transformed into its
integral form. Aiming at a representation for
$\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}_{cs}$
that can manifestly cancel against the $\nu_{n}$ dependences in Eq. (2.48) and
Eq. (2.49), the rapidity regulator $\mathcal{R}(k_{-},k_{+},\tau)$ is
introduced here for bookkeeping purposes. Once the relationship
$\omega<(\sigma+\eta)$ is satisfied, $\mathcal{R}(k_{-},k_{+},\tau)$ always
stays inactive in the limit $\tau\to 0$.

An analogous rearrangement is also applicable to
$\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}_{\bar{c}s}$
in Eq. (2.34). It yields that

$\begin{split}&\sum_{\lambda=0}^{\infty}\frac{{F}^{(\alpha_{n}+\omega-\sigma)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right){F}^{(\alpha_{\bar{n}}+\lambda)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)}{(\omega-\sigma)!\,\lambda!}\frac{\nu_{{\bar{n}}}^{\rho+\lambda-\omega}}{\rho+\lambda-\omega}\bar{\delta}^{\rho+\lambda}_{\omega}\\\
=\,&\frac{{F}^{(\alpha_{n}+\omega-\sigma)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)}{(\omega-\sigma)!}\,\Bigg{\\{}\sum_{\lambda=\mathrm{max}\\{0,(\omega-\rho+1)\\}}^{\infty}\frac{{F}^{(\alpha_{\bar{n}}+\lambda)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)}{\lambda!}\frac{\nu_{{\bar{n}}}^{\rho+\lambda-\omega}}{\rho+\lambda-\omega}\\\
&\hskip
85.35826pt{}-\theta(\omega-\rho-1)\sum_{\lambda=0}^{\omega-\rho-1}\,\frac{{F}^{(\alpha_{\bar{n}}+\lambda)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)}{\lambda!}\lim_{\tau\to
0}\int^{\infty}_{\nu_{\bar{n}}}\frac{{\rm
d}k_{-}}{k_{-}}\mathcal{R}(k_{-},k_{+},\tau)(k_{-})^{\rho+\lambda-\omega}\Bigg{\\}}\,.\end{split}$
(2.56)

We can now recast the logarithmic contributions of
$\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}_{cs}$
and
$\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}_{\bar{c}s}$
in integral form, as presented in Eqs. (2.33) and (2.34). It is worth noting
at this point that the logarithmic terms from both sectors always and only
appear together once $\omega\geq\rho$ and $\omega\geq\sigma$ simultaneously.
Therefore, in the following, we cope with them simultaneously, i.e.,

$\begin{split}&\frac{{F}^{(\alpha_{n}+\omega-\sigma)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right){F}^{(\alpha_{\bar{n}}+\omega-\rho)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)\,}{(\omega-\rho)!\,(\omega-\sigma)!}\Bigg{\\{}\ln\left[\frac{\nu_{n}}{q_{\mathrm{T}}}\right]\,+\,\ln\left[\frac{\nu_{{\bar{n}}}}{q_{\mathrm{T}}}\right]\Bigg{\\}}\\\
=\,&\frac{{F}^{(\alpha_{n}+\omega-\sigma)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right){F}^{(\alpha_{\bar{n}}+\omega-\rho)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)\,}{(\omega-\rho)!\,(\omega-\sigma)!}\\\
&\times\;\Bigg{\\{}-\lim_{\tau\to 0}\int^{\infty}_{\nu_{n}}\frac{{\rm
d}k_{+}}{k_{+}}\,\mathcal{R}(k_{-},k_{+},\tau)+\lim_{\tau\to
0}\int^{\infty}_{0}\frac{{\rm
d}k_{+}}{k_{+}}\,\mathcal{R}(k_{-},k_{+},\tau)-\lim_{\tau\to
0}\int^{\infty}_{\nu_{\bar{n}}}\frac{{\rm
d}k_{-}}{k_{-}}\,\mathcal{R}(k_{-},k_{+},\tau)\Bigg{\\}}\,,\end{split}$ (2.57)

where in the last step we reconvert the logarithms into their integral form
and employ $\mathcal{R}(k_{-},k_{+},\tau)$ to regulate their possibly singular
behaviour in the limit $k_{\pm}\to 0$.

Comparing the above expressions with Eq. (2.48), we observe that terms
containing powers of the auxiliary scale $\nu_{n({\bar{n}})}$ take the same
form in Eqs. (2.55) and (2.56) and Eq. (2.48), and hence cancel upon
combination. Similarly, the dependences of the integral on the auxiliary
scales through their boundaries in Eqs. (2.55)-(2.57) can be assimilated in
their entirety by combining them with the corresponding integrals in Eq.
(2.49). Therefore, we arrive at,

$\displaystyle\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}^{\langle\mathrm{NS}\rangle}_{c\bar{c}}\,=\,\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}^{\langle\mathrm{NS}\rangle}_{c0}\,+\,\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}^{\langle\mathrm{NS}\rangle}_{\bar{c}0}\,,$
(2.58)

where

$\begin{split}\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}^{\langle\mathrm{NS}\rangle}_{c0}\,=\,&\theta\left(\omega-\sigma\right)\,\sum_{\eta=0}^{\omega-\sigma}\,\sum_{\lambda=0}^{\omega-\rho}\,\frac{{F}^{(\alpha_{n}+\eta)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}}+\lambda)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)}{\eta!\,\lambda!}\\\
&\hskip 85.35826pt\times\,\lim_{\tau\to 0}\int^{\infty}_{0}\frac{{\rm
d}k_{+}}{k_{+}}\left(1-\frac{\delta_{\omega}^{\sigma+\eta}}{2}\right)\,\mathcal{R}(k_{-},k_{+},\tau)\,(k_{+})^{\sigma-\omega+\eta}\,\delta^{\omega}_{\lambda+\rho}\,,\\\
\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}^{\langle\mathrm{NS}\rangle}_{\bar{c}0}\,=\,&\theta\left(\omega-\rho\right)\,\sum_{\eta=0}^{\omega-\sigma}\,\sum_{\lambda=0}^{\omega-\rho}\,\frac{{F}^{(\alpha_{n}+\eta)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right){F}^{(\alpha_{\bar{n}}+\lambda)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)}{\eta!\,\lambda!}\\\
&\hskip 85.35826pt\times\,\lim_{\tau\to 0}\int^{\infty}_{0}\frac{{\rm
d}k_{-}}{k_{-}}\left(1-\frac{\delta_{\omega}^{\rho+\lambda}}{2}\right)\,\mathcal{R}(k_{-},k_{+},\tau)\,(k_{-})^{\rho-\omega+\lambda}\,\delta^{\omega}_{\eta+\sigma}\,.\end{split}$
(2.59)

Here, we have divided
$\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}^{\langle\mathrm{NS}\rangle}_{c\bar{c}}$
into two pieces,
$\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}^{\langle\mathrm{NS}\rangle}_{c0}$
and
$\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}^{\langle\mathrm{NS}\rangle}_{\bar{c}0}$,
inheriting the structure from Eq. (2.49). At variance with Eq. (2.49),
however, additional factors $\delta_{\omega}^{\sigma+\eta}$ and
$\delta_{\omega}^{\rho+\lambda}$ are present in the integrands now, induced by
the unbounded integrals of Eq. (2.57). During the derivation, we have evenly
distributed these terms into the $n$-collinear and ${\bar{n}}$-collinear
sectors symmetrically, incurring a factor $\tfrac{1}{2}$ in Eq. (2.59).
Alternative assignments are in principle possible, which may impact the
expressions in the individual sectors but will leave the sum in Eq. (2.58)
invariant.

Equipped with the collinear functions of Eq. (2.44) and the subtrahends of Eq.
(2.59), we are now able to re-express the interior contribution
$\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$ as
follows,

$\displaystyle\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}=\,\sum^{\infty}_{\omega=\rho}\left(q_{\mathrm{T}}^{2}\right)^{\omega}\left[\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{c}-\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}^{\langle\mathrm{NS}\rangle}_{c0}\right]+\sum^{\infty}_{\omega=\sigma}\left(q_{\mathrm{T}}^{2}\right)^{\omega}\left[\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}}-\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}^{\langle\mathrm{NS}\rangle}_{\bar{c}0}\right].$
(2.60)

Herein, we only show the expressions in the NS scheme for brevity, which can
be straightforwardly converted into the FS ones through Eq. (2.54). Referring
back to Eq. (2.40), we have now found a formulation in which the individual
terms on the r.h.s. of Eq. (2.60) are free of any dependence on the auxiliary
scales $\nu_{n({\bar{n}})}$ along the integration path. However, in return,
they are now subject to the choice of rapidity regularisation scheme
implemented through $\mathcal{R}(k_{-},k_{+},\tau)$. As will be discussed in
Secs. 3.2 and 3.3, the regulator $\mathcal{R}(k_{-},k_{+},\tau)$ also depends
on scales $\tilde{\nu}_{n}$ and $\tilde{\nu}_{\bar{n}}$ in practice, which
bears resemblance to those in Eq. (2.26) and can effectively concentrate the
integrand in each sector onto its intrinsic domain, akin to the conventional
dimensional regularisation [232].

Finally, we can now analyse the scaling of $k_{\pm}$ in each sector on the
r.h.s. of Eq. (2.60). Since the collinear functions here are derived by
extrapolating the lower boundaries of Eq. (2.29) and Eq. (2.31), the
integration variables in Eq. (2.44) observe the same scaling rule as those in
momentum space. More explicitly, we have $k_{+}\sim\mathcal{O}(1)$ and
$k_{-}\sim\mathcal{O}(q_{\mathrm{T}}^{2}/m_{H})$ for the $n$-collinear element
and $k_{+}\sim\mathcal{O}(q_{\mathrm{T}}^{2}/m_{H})$ and
$k_{-}\sim\mathcal{O}(1)$ in the ${\bar{n}}$-collinear case. The subtraction
terms in Eq. (2.59), due to the facts that the integration variables $k_{\pm}$
therein have been both expanded in the arguments of ${F}_{i/n,\beta_{n}}$ and
${F}_{j/{\bar{n}},\beta_{\bar{n}}}$ and also that either ${F}_{i/n,\beta_{n}}$
or ${F}_{j/{\bar{n}},\beta_{\bar{n}}}$, respectively, is always synchronised
with that in the collinear function, on the other hand, allows to interpret
the variables $k_{\pm}$ in the dual scaling—one first of all expands the
integrand with the collinear scaling and fixes the $\omega$th-power
correction, which is followed by a posterior power expansion in accordance
with the soft scaling $k_{\pm}\sim\mathcal{O}(q_{\mathrm{T}})$ and the
retention of all contributions below or equal to the $(2\omega)$th-power.

(a) Subtrahends in the soft scaling

(b) Subtrahends in the dual scaling

Figure 1:  Organisation of zero-bin subtrahends in different scalings. Every
point in the map stands for a index-pair $(\lambda,\eta)$ in Eq. (2.59) which
characterises the subtraction procedure entailed by the expansion of
$\widetilde{\mathcal{I}}^{\lambda,\eta}_{[\kappa],\\{\alpha,\beta\\}}$ with
integer $\lambda,\eta\geq-1$ in the small $q_{\mathrm{T}}$ domain. As
examples, the red stars highlight the subleading (NLP) power corrections, the
blue triangles represent the sub-subleading (N2LP) ones, and the green dots
indicate the sub-sub-subleading (N3LP) ones.

It should be emphasised that at a lower power accuracy, the dual scaling can
be reduced to the soft-only one especially in the case $\rho=\sigma=\chi$ with
integer $\chi\geq-1$, which however does not hold in general. As an example,
Fig. 1 illustrates the case $\chi=-1$, which is the minimum of $\chi$ allowed
by the squared amplitudes in Eq. (2.3) and thus concerns the leading power
approximation in the small $q_{\mathrm{T}}$ domain. We observe that, on the
lowest level, both of schemes consist of the origin $\lambda=\eta=-1$ and thus
can be considered to be equivalent here. However, with the power accuracy
growing, the soft scaling forces the indices $\lambda$ and $\eta$ to align
along the diagonals, whereas the dual scaling organises the indices of the
zero-bin subtrahend along rectangular edges, clearly separating the two. Those
distinct patterns can subsequently evaluate to different zero-bin subtrahends
for general choices of rapidity regulators, which therefore defies our
interpretation of the integration variables of Eq. (2.59) in a purely soft
scaling prescription.

In order to systematically implement these scaling rules, we introduce the
power expansion operator $\widehat{\mathbf{T}}_{i}^{(\omega)}$, with the
subscript $i$ running over the sectors $\\{c,\bar{c},s\\}$, acting on the
object following it and projecting out the contribution at the $\omega$th-
power in line with the scalings indicated by $i$. In this way, the components
of Eq. (2.60) can be recast into,

$\begin{split}\left(q_{\mathrm{T}}^{2}\right)^{\omega}\,\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{c}\,&=\,\lim_{\tau\to
0}\,\int^{\tilde{k}^{\mathrm{max}}_{+}}_{0}\,\frac{{\rm
d}k_{+}}{k_{+}}\,\mathcal{R}(k_{-},k_{+},\tau)\,\widehat{\mathbf{T}}_{c}^{(\omega)}\,I^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}\,,\\\
\left(q_{\mathrm{T}}^{2}\right)^{\omega}\,\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}}\,&=\,\lim_{\tau\to
0}\,\int^{\frac{q_{\mathrm{T}}^{2}}{\tilde{k}^{\mathrm{min}}_{+}}}_{0}\,\frac{{\rm
d}k_{-}}{k_{-}}\,\mathcal{R}(k_{-},k_{+},\tau)\,\widehat{\mathbf{T}}_{\bar{c}}^{(\omega)}\,I^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}\,,\end{split}$
(2.61)

and

$\begin{split}\left(q_{\mathrm{T}}^{2}\right)^{\omega}\,\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{c0}\,&=\,\theta\left(\omega-\sigma\right)\lim_{\tau\to
0}\int^{\infty}_{0}\frac{{\rm
d}k_{+}}{k_{+}}\;\mathcal{R}(k_{-},k_{+},\tau)\sum_{\overline{\omega}=(\rho+\sigma)}^{2\omega}\left(1-\frac{\delta_{\overline{\omega}}^{2\omega}}{2}\right)\widehat{\mathbf{T}}_{s}^{(\overline{\omega})}\widehat{\mathbf{T}}_{c}^{(\omega)}I^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}\,,\\\
\left(q_{\mathrm{T}}^{2}\right)^{\omega}\,\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}0}\,&=\,\theta\left(\omega-\rho\right)\lim_{\tau\to
0}\int^{\infty}_{0}\frac{{\rm
d}k_{-}}{k_{-}}\;\mathcal{R}(k_{-},k_{+},\tau)\sum_{\overline{\omega}=(\rho+\sigma)}^{2\omega}\left(1-\frac{\delta_{\overline{\omega}}^{2\omega}}{2}\right)\widehat{\mathbf{T}}_{s}^{(\overline{\omega})}\widehat{\mathbf{T}}_{\bar{c}}^{(\omega)}I^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}\,,\end{split}$
(2.62)

where $I^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$ is the integrand of
$\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$ defined
in Eq. (2.17),

$\displaystyle
I^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}\equiv\,(k_{-})^{\rho}\,(k_{+})^{\sigma}\,{F}^{(\alpha_{n})}_{i/n,\beta_{n}}\left(k_{+}+m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}\left(k_{-}+m_{H}e^{+Y_{H}}\right)\,.$
(2.63)

We note that after the re-arrangements of Eqs. (2.46)-(2.50) and Eqs.
(2.55)-(2.57), the results in Eqs. (2.61)-(2.62) present a form that is
distinct from that observed Eqs. (2.29)-(2.31) as well as Eqs. (2.38)-(2.39)
obtained via the momentum cutoffs $\nu_{n}$ and $\nu_{\bar{n}}$. This is
relevant in particular for the “$sc$” and “$s\bar{c}$” sectors in Eqs.
(2.38)-(2.39), despite hybrid scaling having been used in the transitional
range. As illustrated in Eq. (2.32), the resulting coefficient
$\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}_{t}$
invokes only the integrands at the same power accuracy once the scaling for
the integration variable is given. Contrarily, the situation in Eq. (2.62) is
quite different, where, after applying the soft scaling, the results
$\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}_{c0}$
and
$\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\big{|}_{\bar{c}0}$
comprise not only corrections at the $\omega$th-power, but also those on a
lower level.

In order to interpret this structural difference, it merits recalling that in
deriving Eqs. (2.29) and (2.31) as well as Eqs. (2.38) and (2.39) we have
introduced a set of dynamic regions, see Eq. (2.26), such that within these
restricted phase space domains the assigned scalings are always effective
before and after the phase space integration. However, in calculating the
subtrahends in Eq. (2.62), in particular after the combination with Eqs.
(2.55)-(2.57), all boundaries that separate the different scalings of the
integration variables cancel out. As a result, the “$c0$” and “$\bar{c}0$”
expressions of Eq. (2.62) mainly comprise contributions from the rapidity
extremities $k_{\pm}\to\pm 0$, correlated to the same positions of the
collinear functions in Eq. (2.61), rather than emphasising any integration
segments in the physical domain. In light of this, the scaling laws
implemented in Eq. (2.62) should be conceived more of a prescription guiding
us to organise the zero-bin subtrahends in order to subdue any unphysical
asymptotic behaviour in the collinear sectors and averting the possibility of
these sectors generating any non-trivial remainders for the resulting power
series.

It is worth noting that the integrands in Eq. (2.62) appear to contradict the
homogeneity condition from the method of expansion by regions [211, 212, 213,
214, 174]. In this method, the asymptotic behaviour of the Feynman or phase
space integrals is associated with a set of regions along the integration
path. From this one is able to expand the integrands before completing the
integration and at a given power precision, the resulting power coefficients
only concern the expanded integrands of the same power accuracy, more
specifically,

$\displaystyle\widetilde{\mathcal{F}}^{(\omega)}\sim\sum_{i}\int{\rm
d}\Phi\,\widehat{\mathbf{T}}_{i}^{(\omega)}\,f\,,$ (2.64)

where $\widetilde{\mathcal{F}}^{(\omega)}$ denotes the $\omega$th term in the
power series, with $f$ standing for a given integrand and $\Phi$ encoding a
group of integration variables. In this paper, such a one-to-one
correspondence between the power accuracy of the integrands and that of
$\widetilde{\mathcal{F}}^{(\omega)}$ will be referred to as the homogeneity of
the asymptotic expansion hereafter.

On the contrary, a region analysis where the dimensional regulator uniquely
governs the divergences incurred by the reduced integrands, such homogeneity
is found to generally hold in a variety of kinematical limits, such as the
off-shell large-momentum expansion and the large-mass limits with either an
internal or external heavy parton [214]. The reason comes in part from the
fact that the overlapping sectors for those asymptotic regimes, according to
[174], admit the unbounded integrals over the multiple-expanded integrands,
which, using dimensional regularisation, are all scaleless and thus vanish. In
principle, if the homogeneity is also desired for the small $q_{\mathrm{T}}$
expansion, one can embed a special group of rapidity operators into Eq.
(2.62), such as the analytic [74, 76] and pure rapidity [113, 175] cases, from
which, mimicking the dimensional regulator, the integrals in the zero-bin
subtrahends of Eq. (2.62) become scaleless and are thus eliminated.
Subsequently, homogeneity is restored in the small $q_{\mathrm{T}}$ expansion
and the power coefficients are thereby recast as,

$\displaystyle\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\,=\,\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}^{\langle\mathrm{hom.rap.}\rangle}_{c}\,+\,\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}^{\langle\mathrm{hom.rap.}\rangle}_{\bar{c}}\,.$
(2.65)

This work will dub this set of rapidity regulators the homogeneous ones, as
labeled in the superscripts above. However, the above form is not generic. For
instance, once the exponential regulator [79, 80] is put in place, both the
“$c0$” and “$\bar{c}0$” sectors in Eq. (2.62) will make non-trivial and
indispensable contributions, from which all unphysical singularities generated
in the collinear sectors of Eq. (2.61) can be eliminated at each power,
according to Eq. (2.43). In this sense, the decomposition in Eqs. (2.65) can
only be appropriate for the asymptotic expansion with a homogeneous regulator,
while the generalised recipe, that will hold on a generic choice of rapidity
regulators, is presented in Eq. (2.60).

####  2.4.3 Comparison with existing results

We will now confront the zero-bin subtrahend derived in Eqs. (2.60)-(2.62)
with those proposed in the literature.

At LP, the zero-bin contribution at NLO involves only the master integral
$\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}$
with $\rho=\sigma=\omega=-1$ in the process $pp\to H+X$, and it has been
demonstrated in [159, 160, 161, 162, 165] that the zero-bin subtrahend can be
constructed by taking the soft limit of the integrand. It follows that,

$\displaystyle\frac{1}{q_{\mathrm{T}}^{2}}\widetilde{\mathcal{G}}^{-1,-1,(-1)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}^{\langle\mathrm{NS}\rangle}_{c\bar{c}}\xRightarrow{\text{~{}~{}soft~{}~{}}}\lim_{\tau\to
0}\,\int^{\infty}_{0}\,\frac{{\rm
d}k_{-}}{k_{-}}\,\mathcal{R}(k_{-},k_{+},\tau)\,\frac{{F}^{(\alpha_{n})}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right){F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)}{k_{-}\,k_{+}}\,.$
(2.66)

This result can also be reproduced from Eq. (2.62). At LP, the repeated
application of the expansion operators in Eq. (2.62) does not produce any
power of the integration variables $k_{\pm}$ and therefore yields the same
powers of $k_{\pm}$ as that in Eq. (2.66). Moreover, owing to the relation
$2\omega=\overline{\omega}=-2$, the $\delta$-functions in Eq. (2.62) are both
contributing a factor of $(1/2)$ from each sector, which add up to unity after
combining the results. This thus leads to the exact same expression of Eq.
(2.66).

An alternative scheme to calculate the zero-bin subtraction up to NLP was
proposed in [115, 114]. In this method, the asymptotic behaviour of the small
$q_{\mathrm{T}}$ regime is assumed to be entirely governed by the $n$\- and
${\bar{n}}$-collinear momenta. Subsequently, in order to remove the redundant
overlapping contributions, each collinear ingredient is expanded in accordance
with the scaling of the other one in the opposite direction and then their sum
is averaged, more specifically,

$\begin{split}\left(q_{\mathrm{T}}^{2}\right)^{\omega}\,\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}^{\langle\mathrm{NS}\rangle}_{c\bar{c}}\,&\xRightarrow{\text{~{}~{}opp.col.~{}~{}}}\,\frac{\theta\left(\omega-\rho\right)\theta\left(\omega-\sigma\right)}{2}\,\lim_{\tau\to
0}\,\int^{\infty}_{0}\,\frac{{\rm
d}k_{+}}{k_{+}}\,\mathcal{R}(k_{-},k_{+},\tau)\,\\\ &\hskip
56.9055pt\times\;\Bigg{\\{}\left(\sum_{\overline{\omega}=\sigma}^{\omega}\,\widehat{\mathbf{T}}_{\bar{c}}^{(\overline{\omega})}\right)\widehat{\mathbf{T}}_{c}^{(\omega)}+\left(\sum_{\overline{\omega}=\rho}^{\omega}\,\widehat{\mathbf{T}}_{c}^{(\overline{\omega})}\right)\widehat{\mathbf{T}}_{\bar{c}}^{(\omega)}\\\
&\hskip
79.66771pt{}+\widehat{\mathbf{T}}_{c}^{(\omega)}\,\left(\sum_{\overline{\omega}=\sigma}^{\omega-1}\widehat{\mathbf{T}}_{\bar{c}}^{(\overline{\omega})}\right)+\widehat{\mathbf{T}}_{\bar{c}}^{(\omega)}\left(\sum_{\overline{\omega}=\rho}^{\omega-1}\,\widehat{\mathbf{T}}_{c}^{(\overline{\omega})}\right)\Bigg{\\}}\;I^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}\,.\end{split}$
(2.67)

Using this result to calculate $I^{-1,-1}_{[\kappa],\\{\alpha,\beta\\}}$, we
observe that only the first two terms in the curly brackets can give non-
vanishing contributions at LP as the range of summation in the last two terms
is more restricted. In addition, due to the fact that the product of two
collinear expansion operators in the opposite collinear directions is
equivalent to one single soft expansion operator in this case, Eq. (2.67)
evaluates to the identical expression to that in Eq. (2.66).

At NLP, all terms of Eq. (2.67) contribute to the zero-bin subtraction. While
in the first two terms in the curly brackets the NLP collinear sectors have
been expanded in line with the counting rule in the opposite direction, the
last two terms recover the power corrections that did not contribute to the LP
result, but re-enter the zero-bin subtraction procedure here. Thus, retaining
all terms up to NLP, this yields,

$\begin{split}\widetilde{\mathcal{G}}^{-1,-1,(0)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}^{\langle\mathrm{NS}\rangle}_{c\bar{c}}\,\xRightarrow{\text{~{}~{}opp.col.~{}~{}}}\,&\theta\left(\omega-\rho\right)\,\theta\left(\omega-\sigma\right)\,\lim_{\tau\to
0}\,\int^{\infty}_{0}\,\frac{{\rm
d}k_{+}}{k_{+}}\,\mathcal{R}(k_{-},k_{+},\tau)\,\bigg{\\{}1\,+\,\frac{1}{k_{-}}\,+\frac{1}{k_{+}}\bigg{\\}}\,\\\
&\hskip
113.81102pt\times\;{F}^{(\alpha_{n})}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)\,.\end{split}$
(2.68)

This result agrees with the expectation from the dual scaling. Please note,
LP-like terms have now entered the subtraction term at NLP, constituting the
inhomogeneity we alluded to earlier. More explicitly, substituting the
integrand $I^{-1,-1}_{[\kappa],\\{\alpha,\beta\\}}$ into Eq. (2.62), we
obtain,

$\begin{split}\widetilde{\mathcal{G}}^{-1,-1,(0)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{c0}\,=\,&\theta\left(\omega-\sigma\right)\,\lim_{\tau\to
0}\,\int^{\infty}_{0}\,\frac{{\rm
d}k_{+}}{k_{+}}\,\mathcal{R}(k_{-},k_{+},\tau)\,\bigg{\\{}\frac{1}{2}\,+\,\frac{1}{k_{+}}\bigg{\\}}\,{F}^{(\alpha_{n})}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right),\\\
\widetilde{\mathcal{G}}^{-1,-1,(0)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}_{\bar{c}0}\,=\,&\theta\left(\omega-\rho\right)\,\lim_{\tau\to
0}\,\int^{\infty}_{0}\,\frac{{\rm
d}k_{-}}{k_{-}}\,\mathcal{R}(k_{-},k_{+},\tau)\,\bigg{\\{}\frac{1}{2}\,+\,\frac{1}{k_{-}}\bigg{\\}}\,{F}^{(\alpha_{n})}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}})}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right).\end{split}$
(2.69)

Recalling that additional constraints on the exponents $\rho$ and $\sigma$
have already been imposed in the definition of the power series in Eq. (2.60),
combining both contributions above precisely reproduces the results of Eq.
(2.68).

Furthermore, even though Eq. (2.68) is obtained by summarising the asymptotic
properties of the NLP ingredients in [115, 114], it is very interesting to
note that Eq. (2.68) is in fact still useful in organising the overlap removal
for N2LP and beyond. To see this, we note that the collinear expansion
operators $\widehat{\mathbf{T}}_{c}$ and $\widehat{\mathbf{T}}_{\bar{c}}$
commute. Hence, we can drop the factor of $\tfrac{1}{2}$ and recast the r.h.s.
of Eq. (2.68) below,

$\begin{split}\left(q_{\mathrm{T}}^{2}\right)^{\omega}\,\widetilde{\mathcal{G}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}\Bigg{|}^{\langle\mathrm{NS}\rangle}_{c\bar{c}}\,\xRightarrow{\text{~{}~{}opp.col.~{}~{}}}&\,\theta\left(\omega-\rho\right)\theta\left(\omega-\sigma\right)\,\lim_{\tau\to
0}\,\int^{\infty}_{0}\,\frac{{\rm
d}k_{+}}{k_{+}}\,\mathcal{R}(k_{-},k_{+},\tau)\,\\\
&\times\Bigg{\\{}\widehat{\mathbf{T}}_{\bar{c}}^{(\omega)}\widehat{\mathbf{T}}_{c}^{(\omega)}+\widehat{\mathbf{T}}_{c}^{(\omega)}\,\left(\sum_{\overline{\omega}=\sigma}^{\omega-1}\,\widehat{\mathbf{T}}_{\bar{c}}^{(\overline{\omega})}\,\right)+\widehat{\mathbf{T}}_{\bar{c}}^{(\omega)}\,\left(\sum_{\overline{\omega}=\rho}^{\omega-1}\,\widehat{\mathbf{T}}_{c}^{(\overline{\omega})}\,\right)\Bigg{\\}}I^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}\,,\end{split}$
(2.70)

Plugging the expressions of eq. (2.63) into eq. (2.70), we then obtain,

$\begin{split}\mathrm{r.h.s.~{}of~{}eq.~{}\eqref{eq:def:c0:cbar0:ns:oppcol2}}=&\;\theta\left(\omega-\rho\right)\theta\left(\omega-\sigma\right)\,\lim_{\tau\to
0}\,\int^{\infty}_{0}\,\frac{{\rm
d}k_{+}}{k_{+}}\,\mathcal{R}(k_{-},k_{+},\tau)\\\
&\times\Bigg{\\{}k_{-}^{\omega}\,k_{+}^{\omega}\,\frac{{F}^{(\alpha_{n}+\omega-\sigma)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)}{(\omega-\sigma)!}\frac{{F}^{(\alpha_{\bar{n}}+\omega-\rho)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)}{(\omega-\rho)!}\\\
&\hskip
17.07164pt{}+k_{-}^{\omega}\,\sum_{\eta=0}^{\omega-\sigma-1}\,k_{+}^{\sigma+\eta}\,\frac{{F}^{(\alpha_{n}+\eta)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}}+\omega-\rho)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)}{\eta!\,(\omega-\rho)!}\,\\\
&\hskip
17.07164pt{}+k_{+}^{\omega}\,\sum_{\eta=0}^{\omega-\sigma-1}\,k_{-}^{\rho+\lambda}\,\frac{{F}^{(\alpha_{n}+\omega-\sigma)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right){F}^{(\alpha_{\bar{n}}+\lambda)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)}{\lambda!(\omega-\sigma)!}\,\Bigg{\\}}\,.\end{split}$
(2.71)

Extracting a common factor of $(q_{\mathrm{T}}^{2})^{\omega}$ from the curly
bracket above, we observe its coincidence with the sum of Eq. (2.59) and thus
the equivalence with Eq. (2.69).

We will now compare our results with the results obtained through employing
the method of expansion by regions [211, 212, 213, 214, 174] and focus in
particular on the formalism proposed in [174]. Therein, the mathematical
foundation behind the expansion by regions has been discussed via examples for
the one-loop integrals in different kinematic limits. Of them, the analysis of
the Feynman integrals in the Sudakov limit [174] is intimately related to the
$q_{\mathrm{T}}$ spectra within the asymptotic regime. Therefore, after
appropriate adaptations, the techniques proposed in [174] can also be
exploited to compute the zero-bin subtraction here. In the following, we will
elucidate this application.

The starting point of the formalism in [174] is to work out a group of dynamic
regions that are able to encompass the entire phase space and also separately
accommodate the contributions from all involved scales. In our case, this goal
can be accomplished by the dynamic modes presented in Eqs. (2.26). The next
task is to extrapolate the bounded integrals contained in each region so as to
remove all auxiliary cutoff scales along the integration paths. This procedure
can be illustrated by using the expansion operators introduced in Eqs. (2.61)
and (2.62), i.e.

$\displaystyle\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}=\sum_{\omega}\left[\int^{\tilde{k}^{\mathrm{max}}_{+}}_{\nu_{n}}\widehat{\mathbf{T}}_{c}^{(\omega)}\,+\int^{\nu_{n}}_{\frac{q_{\mathrm{T}}^{2}}{\nu_{\bar{n}}}}\,\widehat{\mathbf{T}}_{s}^{(\omega)}\,+\int^{\frac{q_{\mathrm{T}}^{2}}{\nu_{\bar{n}}}}_{\tilde{k}^{\mathrm{min}}_{+}}\,\widehat{\mathbf{T}}_{\bar{c}}^{(\omega)}\right]\,\frac{{\rm
d}k_{+}}{k_{+}}\,I^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}\,,$ (2.72)

and then extrapolating the integral boundaries $\nu_{n}$ and $\nu_{\bar{n}}$,

$\begin{split}\sum_{\omega_{c}}\int^{\tilde{k}^{\mathrm{max}}_{+}}_{\nu_{n}}\widehat{\mathbf{T}}_{c}^{(\omega_{c})}&\rightarrow\lim_{\tau\to
0}\left[\sum_{\omega_{c}}\int^{\tilde{k}^{\mathrm{max}}_{+}}_{0}\mathcal{R}\,\widehat{\mathbf{T}}_{c}^{(\omega_{c})}-\sum_{\omega_{c},\omega_{s}}\int^{\nu_{n}}_{\frac{q_{\mathrm{T}}^{2}}{\nu_{\bar{n}}}}\mathcal{R}\,\widehat{\mathbf{T}}_{s}^{(\omega_{s})}\widehat{\mathbf{T}}_{c}^{(\omega_{c})}-\sum_{\omega_{c},\omega_{\bar{c}}}\int^{\frac{q_{\mathrm{T}}^{2}}{\nu_{\bar{n}}}}_{0}\mathcal{R}\,\widehat{\mathbf{T}}_{\bar{c}}^{(\omega_{\bar{c}})}\widehat{\mathbf{T}}_{c}^{(\omega_{c})}\right],\\\
\sum_{\omega_{s}}\int^{\nu_{n}}_{\frac{q_{\mathrm{T}}^{2}}{\nu_{\bar{n}}}}\widehat{\mathbf{T}}_{s}^{(\omega_{s})}&\rightarrow\lim_{\tau\to
0}\left[\sum_{\omega_{s}}\int^{\infty}_{0}\mathcal{R}\,\widehat{\mathbf{T}}_{s}^{(\omega_{s})}-\sum_{\omega_{c},\omega_{s}}\int^{\infty}_{\nu_{n}}\mathcal{R}\,\widehat{\mathbf{T}}_{c}^{(\omega_{c})}\widehat{\mathbf{T}}_{s}^{(\omega_{s})}-\sum_{\omega_{s},\omega_{\bar{c}}}\int^{\frac{q_{\mathrm{T}}^{2}}{\nu_{\bar{n}}}}_{0}\mathcal{R}\,\widehat{\mathbf{T}}_{\bar{c}}^{(\omega_{\bar{c}})}\widehat{\mathbf{T}}_{s}^{(\omega_{s})}\right],\\\
\sum_{\omega_{\bar{c}}}\int^{\frac{q_{\mathrm{T}}^{2}}{\nu_{\bar{n}}}}_{\tilde{k}^{\mathrm{min}}_{+}}\widehat{\mathbf{T}}_{\bar{c}}^{(\omega_{\bar{c}})}&\rightarrow\lim_{\tau\to
0}\left[\sum_{\omega_{\bar{c}}}\int^{\infty}_{\tilde{k}^{\mathrm{min}}_{+}}\mathcal{R}\,\widehat{\mathbf{T}}_{\bar{c}}^{(\omega_{\bar{c}})}-\sum_{\omega_{\bar{c}},\omega_{s}}\int^{\nu_{n}}_{\frac{q_{\mathrm{T}}^{2}}{\nu_{\bar{n}}}}\mathcal{R}\,\widehat{\mathbf{T}}_{s}^{(\omega_{s})}\widehat{\mathbf{T}}_{\bar{c}}^{(\omega_{\bar{c}})}-\sum_{\omega_{c},\omega_{\bar{c}}}\int^{\infty}_{\nu_{n}}\mathcal{R}\,\widehat{\mathbf{T}}_{c}^{(\omega_{c})}\widehat{\mathbf{T}}_{\bar{c}}^{(\omega_{\bar{c}})}\right].\end{split}$
(2.73)

Therein, in order to render the integrals still well-defined after
extrapolating the boundaries, the rapidity regulator $\mathcal{R}$ has been
put in place. Since the discussion here is focused on the conservative
rapidity regularisation prescription that preserves the expression of
$\mathcal{R}$ in all the involved sectors, in writing Eq. (2.73), we pull all
regulators $\mathcal{R}$ in front of the expansion operators. We note that the
r.h.s. of Eq. (2.73) still depends on the cutoff scales $\nu_{n}$ and
$\nu_{\bar{n}}$. This dependence can be eliminated through the following
identity,

$\begin{split}\sum_{\omega_{c},\omega_{\bar{c}},\omega_{s}}\int^{\infty}_{0}\mathcal{R}\,\widehat{\mathbf{T}}_{c}^{(\omega_{c})}\widehat{\mathbf{T}}_{s}^{(\omega_{s})}\widehat{\mathbf{T}}_{\bar{c}}^{(\omega_{\bar{c}})}\rightarrow&\sum_{\omega_{c},\omega_{s}}\int^{\frac{q_{\mathrm{T}}^{2}}{\nu_{\bar{n}}}}_{0}\mathcal{R}\,\widehat{\mathbf{T}}_{s}^{(\omega_{s})}\widehat{\mathbf{T}}_{c}^{(\omega_{c})}+\sum_{\omega_{c},\omega_{\bar{c}}}\int_{\frac{q_{\mathrm{T}}^{2}}{\nu_{\bar{n}}}}^{\nu_{n}}\mathcal{R}\,\widehat{\mathbf{T}}_{\bar{c}}^{(\omega_{\bar{c}})}\widehat{\mathbf{T}}_{c}^{(\omega_{c})}\\\
&{}+\sum_{\omega_{\bar{c}},\omega_{s}}\int^{\infty}_{\nu_{n}}\mathcal{R}\,\widehat{\mathbf{T}}_{s}^{(\omega_{s})}\widehat{\mathbf{T}}_{\bar{c}}^{(\omega_{\bar{c}})}\,.\end{split}$
(2.74)

In its derivation, we have used the commutativity of the power projection
operators $\widehat{\mathbf{T}}_{c}$, $\widehat{\mathbf{T}}_{s}$, and
$\widehat{\mathbf{T}}_{\bar{c}}$. Combining Eq. (2.74) with Eq. (2.73), we
have,

$\begin{split}\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}=&\;\lim_{\tau\to
0}\bigg{[}\sum_{\omega_{c}}\int^{\tilde{k}^{\mathrm{max}}_{+}}_{0}\mathcal{R}\,\widehat{\mathbf{T}}_{c}^{(\omega_{c})}\,+\sum_{\omega_{s}}\int^{\infty}_{0}\mathcal{R}\,\widehat{\mathbf{T}}_{s}^{(\omega_{s})}+\sum_{\omega_{\bar{c}}}\int^{\infty}_{\tilde{k}^{\mathrm{min}}_{+}}\mathcal{R}\,\widehat{\mathbf{T}}_{\bar{c}}^{(\omega_{\bar{c}})}\\\
&\hskip
28.45274pt{}-\sum_{\omega_{c},\omega_{s}}\int^{\infty}_{0}\mathcal{R}\,\widehat{\mathbf{T}}_{c}^{(\omega_{c})}\,\widehat{\mathbf{T}}_{s}^{(\omega_{s})}-\sum_{\omega_{c},\omega_{\bar{c}}}\int^{\infty}_{0}\mathcal{R}\,\widehat{\mathbf{T}}_{c}^{(\omega_{c})}\,\widehat{\mathbf{T}}_{\bar{c}}^{(\omega_{\bar{c}})}-\sum_{\omega_{s},\omega_{\bar{c}}}\int^{\infty}_{0}\mathcal{R}\,\widehat{\mathbf{T}}_{s}^{(\omega_{s})}\,\widehat{\mathbf{T}}_{\bar{c}}^{(\omega_{\bar{c}})}\\\
&\hskip
28.45274pt{}+\sum_{\omega_{c},\omega_{s},\omega_{\bar{c}}}\int^{\infty}_{0}\mathcal{R}\,\widehat{\mathbf{T}}_{c}^{(\omega_{c})}\,\widehat{\mathbf{T}}_{s}^{(\omega_{s})}\,\widehat{\mathbf{T}}_{\bar{c}}^{(\omega_{\bar{c}})}\,\bigg{]}\,\frac{{\rm
d}k_{+}}{k_{+}}\,I^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}\,.\end{split}$
(2.75)

Taking into account the facts that expanding in the collinear scalings
following the soft power-projection operator $\widehat{\mathbf{T}}_{s}$ can
only result in delta functions as well as that two consecutive expansions in
the collinear scalings are equivalent to one single expansion in the soft
scaling upon having summed all power contributions, we see that the doubly and
triply expanded terms in Eq. (2.75) are equal to an expansion in the soft
scaling. To this end, we obtain,

$\displaystyle\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}=$
$\displaystyle\;\sum_{\omega}\,\lim_{\tau\to
0}\,\bigg{[}\int^{\tilde{k}^{\mathrm{max}}_{+}}_{0}\mathcal{R}\,\widehat{\mathbf{T}}_{c}^{(\omega)}\,+\int^{\infty}_{\tilde{k}^{\mathrm{min}}_{+}}\mathcal{R}\,\widehat{\mathbf{T}}_{\bar{c}}^{(\omega)}-\int^{\infty}_{0}\mathcal{R}\,\widehat{\mathbf{T}}_{s}^{(\omega)}\bigg{]}\,\frac{{\rm
d}k_{+}}{k_{+}}\,I^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}\,.$ (2.76)

Comparing with Eqs. (2.61) and (2.62), we observe that the first two terms in
Eq. (2.76) and the collinear functions in Eq. (2.61) are equivalent. To
explore the relationship of the third term in Eq. (2.76) and the sum of Eq.
(2.62), we note that after bringing in the expression in Eq. (2.63) the soft
sector in Eq. (2.76) evaluates to an infinite sum of the unbounded integrals,
i.e.,

$\begin{split}\sum_{\eta,\lambda=0}^{\infty}\,\lim_{\tau\to
0}\,\int_{0}^{\infty}\,\frac{{\rm
d}k_{+}}{k_{+}}\,\mathcal{R}\left(k_{-},k_{+},\tau\right)\,\frac{k_{+}^{\sigma+\eta}\,k_{-}^{\rho+\lambda}}{\eta!\,\lambda!}\,{F}^{(\alpha_{n}+\eta)}_{i/n,\beta_{n}}\left(m_{H}e^{-Y_{H}}\right)\,{F}^{(\alpha_{\bar{n}}+\lambda)}_{j/{\bar{n}},\beta_{\bar{n}}}\left(m_{H}e^{+Y_{H}}\right)\,,\end{split}$
(2.77)

which is exactly equal to the result from Eq. (2.62) upon summing up all power
contributions. Therefore, we can conclude at least the correspondence of Eq.
(2.76) and Eqs. (2.61) and (2.62) on the cumulative level.

Notwithstanding, once a power-by-power comparison is being considered between
both results, the analysis becomes more subtle. In [174], the analytic
regulator is applied throughout the calculation of the Sudakov form factors.
Implementing this scheme into Eq. (2.75) (or Eq. (2.76)), the integrals that
are expanded in the soft scaling or comprise double or triple expansions all
vanish. Then, Eq. (2.75) and Eq. (2.76) both reduce into the same
decomposition we found in Eq. (2.65) induced by the homogeneous regulators.
However, if the inhomogeneous regularisation scheme is of one’s particular
interest, the soft function and the doubly and triply expanded constituents of
Eq. (2.75) can make non-trivial contributions at a given power accuracy, for
which an unambiguous counting rule is required in order to produce the correct
power coefficient
$\widetilde{\mathcal{I}}^{\rho,\sigma,(\omega)}_{[\kappa],\\{\alpha,\beta\\}}$.

###  2.5 Discussion and extrapolation

In the previous subsections, we have carried out the power expansion of the
double-differential observable ${{\rm d}\sigma_{H}}/{{\rm d}Y_{H}{\rm
d}q_{\mathrm{T}}^{2}}$ in the low $q_{\mathrm{T}}$ domain. The result reads,

$\displaystyle\frac{{\rm d}\sigma_{H}}{{\rm d}Y_{H}{\rm
d}q_{\mathrm{T}}^{2}}\,=\,\lambda_{t}^{2}\,\sum_{\omega,\overline{\omega}}\,\left({q_{\mathrm{T}}^{2}}\right)^{\omega+\overline{\omega}}\,\sum_{[\kappa]}\,\sum_{\\{\alpha,\beta\\},\rho,\sigma}\,\widetilde{\mathcal{H}}^{(\omega),\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}\left(m_{H},Y_{H},s\right)\,\bigg{\\{}\widetilde{\mathcal{I}}^{\rho,\sigma,(\overline{\omega})}_{[\kappa],\\{\alpha,\beta\\}}\,+\,\Delta\widetilde{\mathcal{I}}^{\rho,\sigma,(\overline{\omega})}_{[\kappa],\\{\alpha,\beta\\}}\bigg{\\}}\,.$
(2.78)

Therein,
$\widetilde{\mathcal{H}}^{(\omega),\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$
is a function of the hard scales $m_{H}$ and $s$. Its expression has been
presented in Eq. (2.12) on a power-by-power basis.
$\Delta\widetilde{\mathcal{I}}^{\rho,\sigma,(\overline{\omega})}_{[\kappa],\\{\alpha,\beta\\}}$
accounts for the power corrections induced by the boundary conditions of the
phase space integral. The results in the $n$-ultra-collinear limit can be
extracted from Eqs. (2.22)-(2.25), while those in the ${\bar{n}}$-direction
can be derived by exchanging the light-cone coordinates in Eqs. (2.24)-(2.25),
as appropriate. Finally, Eq. (2.78) contains
$\widetilde{\mathcal{I}}^{\rho,\sigma,(\overline{\omega})}_{[\kappa],\\{\alpha,\beta\\}}$
which comprises the contribution from the bulk of the phase space. In Sec.
2.4.1 and Sec. 2.4.2, we made use of two strategies to evaluate the expansion
coefficients
$\widetilde{\mathcal{I}}^{\rho,\sigma,(\overline{\omega})}_{[\kappa],\\{\alpha,\beta\\}}$.
In Sec. 2.4.1, owing to the fact that our integration region comprises
multiple scales, we re-categorised the integration path via a two of auxiliary
cutoff scales, $\nu_{n}$ and $\nu_{\bar{n}}$. A dedicated method is applied
for each subregion to perform the expansion within the low $q_{\mathrm{T}}$
domain. The quoted power series following this approach is illustrated in Eq.
(2.40). Then, in order to pave the way for future investigations of the power
expansion at N2LO and beyond, we recombined the individual terms of Eq.
(2.40), eliminating any dependence on the auxiliary scales. The rearranged
expansion coefficients are given in Eq. (2.60). Equipped therewith, we are
able to calculate the power series of ${{\rm d}\sigma_{H}}/{{\rm d}Y_{H}\,{\rm
d}q_{\mathrm{T}}^{2}}$ at an arbitrary power accuracy.

During our analysis, in order to facilitate the establishment of the relevant
scaling relations, we concentrated on the central rapidity region where
$e^{\pm Y_{H}}\sim\mathcal{O}(1)$. As it turns out, the expression in Eq.
(2.78) can also be extrapolated to the domain

$\displaystyle\mathrm{min}\\{\sqrt{s}-m_{H}e^{Y_{H}},\sqrt{s}-m_{H}e^{-Y_{H}}\\}>\mathrm{max}\\{m_{H}e^{Y_{H}},m_{H}e^{-Y_{H}}\\}\,.$
(2.79)

This extrapolation is straightforward to achieve in the kinematical reduction
in Sec. 2.2, since the expansion parameter therein is
$(q_{\mathrm{T}}^{2}/m_{H}^{2})$.

Regarding
$\Delta\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$,
it is worth reminding the reader that the power expansion in Sec. 2.3
primarily builds on the hierarchies,

$\displaystyle\sqrt{s}-m_{H}e^{\pm
Y_{H}}\gg\left|m_{H}-m_{\mathrm{T}}\right|e^{\pm Y_{H}}\,,\quad m_{H}e^{\pm
Y_{H}}\gg\frac{q_{\mathrm{T}}^{2}}{\sqrt{s}-m_{H}e^{\mp Y_{H}}}\,.$ (2.80)

Both relations are maintained by Eq. (2.79) in the low $q_{\mathrm{T}}$
domain, and therefore the power expansion of
$\Delta\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$
proceeds as before. A similar situation can be found in the interior
contribution
$\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$, in
which the categorisation and the subsequent expansion base on the following
relationship,

$\displaystyle\sqrt{s}-m_{H}e^{-Y_{H}}>m_{H}e^{-Y_{H}}\gtrsim\nu_{n}\gg\frac{q_{\mathrm{T}}^{2}}{\nu_{\bar{n}}}\gtrsim\frac{q_{\mathrm{T}}^{2}}{m_{H}e^{Y_{H}}}>\frac{q_{\mathrm{T}}^{2}}{\sqrt{s}-m_{H}e^{Y_{H}}}\,.$
(2.81)

These inequalities ensure the integration variable $k_{\pm}$ is always
comparable to $m_{H}e^{\mp Y_{H}}$ in magnitude in the collinear sector, while
that in the transitional domain is small enough to perform the power
expansion. Eq. (2.81) is manifestly satisfied by the domain of Eq. (2.79) for
the small $q_{\mathrm{T}}$ regime, from which there exists a sufficiently wide
window to introduce the auxiliary cutoff scales $\nu_{n}$ and $\nu_{\bar{n}}$.
Hence, the subsequent power series of
$\widetilde{\mathcal{I}}^{\rho,\sigma}_{[\kappa],\\{\alpha,\beta\\}}$ can be
derived as in Sec. 2.4.

Following the above consideration, we can implement the expansion of Eq.
(2.78) within the full range of Eq. (2.79). At the LHC with a colliding energy
$\sqrt{s}=13\,\text{TeV}$, this interval corresponds to the rapidity region
$|Y_{H}|\lesssim 4$ and, thus, enables a reliable power expansion in the bulk
of the accessible phase space. Moving from $|Y_{H}|\approx 0$ to
$|Y_{H}|\lesssim 4$, it is possible that novel hierarchies emerge from