# Projected transverse momentum resummation in top-antitop pair production at
LHC

Wan-Li Ju, Marek Schönherr

###### Abstract

The transverse momentum distribution of the ${t{\bar{t}}}$ system is of both
experimental and theoretical interests. In the presence of the azimuthally
asymmetric divergences, pursuing the resummation at the high logarithmic
precision is rather demanding in general. In this paper, we propose the
projected transverse momentum spectrum ${\rm d}\sigma_{{t{\bar{t}}}}/{\rm
d}q_{\tau}$, which is derived from the classical $\vec{q}_{\mathrm{T}}$
spectrum by integrating out the rejection component $q_{\tau_{\perp}}$ with
respect to a reference unit vector $\vec{\tau}$, to serve as an alternative
resolution to curb those asymmetric influences, in additional to the
azimuthally averaged case ${\rm d}\sigma_{{t{\bar{t}}}}/{\rm
d}|\vec{q}_{\mathrm{T}}|$. In the context of the effective field theories,
SCETII and HQET, we will demonstrate that in spite of the $q_{\tau_{\perp}}$
integrations, the leading asymptotic terms of ${\rm
d}\sigma_{{t{\bar{t}}}}/{\rm d}q_{\tau}$ still observe the factorization
pattern in terms of the hard, beam, and soft functions in the vicinity of
$q_{\tau}=0$ GeV. With the help of the renormalisation group equation
techniques, we evaluate the resummed spectra ${\rm
d}\sigma_{{t{\bar{t}}}}/{\rm d}q_{\mathrm{T,in}}$, ${\rm
d}\sigma_{{t{\bar{t}}}}/{\rm d}q_{\mathrm{T,out}}$, and ${\rm
d}\sigma_{{t{\bar{t}}}}/{\rm d}\Delta\phi_{{t{\bar{t}}}}$ at NLL+NLO,
N2LL+N2LO, and approximate N2LL′+N2LO accuracy. This is the first time the
azimuthal spectrum $\Delta\phi_{{t{\bar{t}}}}$ is appraised at or beyond the
N2LL level including a consistent treatment of both beam collinear and soft
radiation.

###### Contents

  1. 1 Introduction
  2. 2 Factorisation
    1. 2.1 Kinematics and the factorised cross section
    2. 2.2 Dynamic regions
    3. 2.3 The case of $N_{\mathrm{J}}=0$
    4. 2.4 The case of $N_{\mathrm{J}}\geq 1$
      1. 2.4.1 The $N_{\mathrm{J}}=1$ configuration
      2. 2.4.2 The $N_{\mathrm{J}}\geq 2$ configuration
      3. 2.4.3 Summary and discussion
    5. 2.5 The soft function with the exponential regulator
  3. 3 Resummation
    1. 3.1 Asymptotic behavior
    2. 3.2 Evolution equations
    3. 3.3 Observables
    4. 3.4 Matching to fixed-order QCD
  4. 4 Numerical Results
    1. 4.1 Parameters and uncertainty estimates
    2. 4.2 Validation
    3. 4.3 Resummation improved results
  5. 5 Conclusions
  6. A The double differential transverse momentum distribution ${\rm d}^{2}\sigma_{{t{\bar{t}}}}/{\rm d}^{2}\vec{q}_{\mathrm{T}}$
  7. B Evolution kernel for the non-diagonal anomalous dimension

##  1 Introduction

Hadroproduction of top-antitop pairs plays a pivotal part in the physics
programme of the LHC experiments due to its role in the precise extractions of
fundamental parameters of the Standard Model (SM). It has thus drawn plenty of
theoretical and experimental attention in the recent years. On the
experimental side, the inclusive top-pair production cross section has been
measured at the colliding energies $\sqrt{s}=5.02~{}$TeV [1, 2, 3] and
$\sqrt{s}=13~{}$TeV [4, 5, 6], respectively, whilst a large number of the
differential spectra have been published in the latest analyses [7, 8, 9, 10,
11, 12], including the transverse momentum of the ${t{\bar{t}}}$ system
$q_{\mathrm{T}}$, the invariant mass of the top quark pair $M_{t{\bar{t}}}$,
and the azimuthal opening angle of the top and antitop quarks
$\Delta\Phi_{{t{\bar{t}}}}$. Theoretical calculations of these spectra also
have since long attracted a lot of interest in the community. While NLO QCD
corrections to top-pair production were determined already some time ago [13,
14, 15, 16], recent advances have reached the NNLO accuracy [17, 18, 19, 20,
21, 22, 23]. Top-quark decay effects were considered in [24, 25, 26, 27, 28]
and the electroweak (EW) corrections in [29, 30, 31, 32, 33, 34, 35, 36, 21].
Along with the progress made in fixed-order calculations, in a bid to improve
the perturbative convergence and in turn the predictivity of the theoretical
results, resummed calculations have also been carried out within a variety of
frameworks and the kinematical limits. Examples include the mechanic threshold
[37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47], the top-quark pair production
threshold [48, 49, 50, 51, 52, 53, 54], the low transverse momentum domain
[55, 56, 57, 58, 59], and the narrow jettiness regime [60]. Very recently, the
combination of the fixed-order results with a parton shower has been discussed
in [61, 62].

This work will investigate the projection of the ${t{\bar{t}}}$ system’s
transverse momentum with respect to a reference unit vector $\vec{\tau}$ on
the azimuthal plane, more explicitly,
$q_{\tau}\equiv|\vec{q}_{\tau_{\|}}|\equiv|\vec{q}_{\mathrm{T}}\cdot\vec{\tau}|$.
In contrast to the traditional transverse momentum spectrum ${\rm
d}\sigma_{{t{\bar{t}}}}/{\rm d}q_{\mathrm{T}}$,
$q_{\mathrm{T}}=|\vec{q}_{\mathrm{T}}|$, where both components of
$\vec{q}_{\mathrm{T}}$ are measured and thereby constrained, the present
observable ${\rm d}\sigma_{{t{\bar{t}}}}/{\rm d}q_{\tau}$ only concerns the
projected piece $\vec{q}_{\tau_{\|}}$, leaving the perpendicular part
$\vec{q}_{\tau_{\perp}}$ unresolved and, hence, it should be integrated out.
As will be demonstrated in this paper, the act of integrating out the
perpendicular components will introduce new and distinguishing features to the
$q_{\tau}$ spectrum, particularly in regards to the treatment on the azimuthal
asymmetric contributions [63, 64, 57, 58].

Induced by the soft and collinear radiation, the fixed-order calculation of
the $q_{\tau}$ distribution exhibits substantial higher-order corrections in
the small $q_{\tau}$ region. This, thus, necessitates a resummation of the
dominant contributions in this regime to all orders to stabilise the
perturbative predictions. In order to accomplish this target, one of the
prerequisite conditions is to determine the dynamic regions driving the
asymptotic behaviour in the limit $q_{\tau}\to 0$. For the classic transverse
momentum resummation, this analysis was first presented for Drell-Yan
production in [65] by means of inspecting the power laws of a generic
configuration on the pinch singularity surface [66, 67, 68]. It was proven
that the leading singular behaviour in the small $q_{\mathrm{T}}$ domain was
well captured by the beam-collinear, soft, and hard regions. However, this
conclusion cannot be straightforwardly to be applied onto the $q_{\tau}$
resummation in top-pair production, as the deep recoil configuration
$|\vec{q}_{\tau_{\perp}}|\sim M_{{t{\bar{t}}}}\gg q_{\tau}$, which stems from
the integral over the perpendicular component, was kinematically excluded in
[65]. Therefore, for delivering an honest and self-consistent study on ${\rm
d}\sigma_{{t{\bar{t}}}}/{\rm d}q_{\tau}$, this work will reappraise the
scalings of the relevant configurations, comprising both the
$q_{\mathrm{T}}$-like configuration
$M_{t{\bar{t}}}\gg|\vec{q}_{\tau_{\perp}}|\sim q_{\tau}$ and the asymmetric
one $|\vec{q}_{\tau_{\perp}}|\sim M_{t{\bar{t}}}\gg q_{\tau}$.

To this end, we will exploit the strategy of expansion by regions [69, 70, 71,
72] to motivate the momentum modes governing the low $q_{\tau}$ regime, which
will cover the beam-collinear, soft, central-jet, and hard regions. Then, the
soft-collinear effective theory (SCET) [73, 74, 75, 76, 77, 78, 79, 80, 81,
82] and the heavy quark effective theory (HQET) [83, 84, 85, 86] are used to
embody those dynamic modes, thereby calculating the effective amplitudes and
the respective differential cross sections. After carrying out a multipole
expansion, the results constructed by those dynamic regions all reflect the
unambiguous scaling behaviors, which can be determined from the power
prescriptions of the relevant effective fields. From its results, we point out
that the leading asymptotic behavior of ${\rm d}\sigma_{{t{\bar{t}}}}/{\rm
d}q_{\tau}$ is still resultant of the symmetric configuration
$M_{t{\bar{t}}}\gg|\vec{q}_{\tau_{\perp}}|\sim q_{\tau}$, which is in practice
dictated by the beam-collinear, soft, and hard momenta, akin to the case of
${\rm d}\sigma_{{t{\bar{t}}}}/{\rm d}q_{\mathrm{T}}$, whereas the
contributions from the $|\vec{q}_{\tau_{\perp}}|\sim M_{t{\bar{t}}}\gg
q_{\tau}$ pattern are suppressed by at least one power of
$\lambda_{\tau}\equiv q_{\tau}/M_{t{\bar{t}}}$.

Upon the identification of the leading regions, we make use of the decoupling
properties of the soft modes [73, 49] to derive the factorisation formula for
${\rm d}\sigma_{{t{\bar{t}}}}/{\rm d}q_{\tau}$. Owing to the integration over
$\vec{q}_{\tau_{\perp}}$, the impact space integrals herein are all reduced
from 2D to 1D. Thus, the azimuthal asymmetric contributions, which in
principle contain divergent terms in the asymptotic regime in the general
${\rm d}\sigma_{{t{\bar{t}}}}/{\rm d}\vec{q}_{\mathrm{T}}$ cross section after
completing the inverse Fourier transformations, do not contribute any
divergences to the $q_{\tau}$ spectrum. This is the second
$\vec{q}_{\mathrm{T}}$-based observable free of asymmetric singularities in
addition to the azimuthally averaged spectra ${\rm
d}\sigma_{{t{\bar{t}}}}/{\rm d}q_{\mathrm{T}}$ [55, 56, 59].

To implement the resummation, we employ the renormalisation group equations
(RGE) and the rapidity renormalisation group equations (RaGE) to evolve the
intrinsic scales in the respective ingredients and in turn accomplish the
logarithmic exponentiations [87, 88, 89, 90]. Alternative approaches can also
be found in [91, 92, 93, 94, 95, 96, 97, 98, 99, 100]. For assessing the
resummation accuracy, we take the logarithmic counting rule
$L_{\mathrm{M}}\sim\alpha_{s}^{-1}\sim\lambda_{L}^{-1}$ throughout and
organise the perturbative corrections to the relevant sectors in line with the
following prescription,

$\begin{split}\frac{{\rm d}{\sigma}_{{t{\bar{t}}}}}{{\rm
d}q_{\tau}}\sim&\;\sigma^{\text{Born}}_{{t{\bar{t}}}}\,\exp\bigg{[}\;\underbrace{L_{\mathrm{M}}\,f_{0}(\alpha_{s}L_{\mathrm{M}})}_{(\text{LL})}+\underbrace{f_{1}(\alpha_{s}L_{\mathrm{M}})}_{(\text{NLL,NLL${}^{\prime}$})}+\underbrace{\alpha_{s}\,f_{2}(\alpha_{s}L_{\mathrm{M}})}_{(\text{N${}^{2}$LL,N${}^{2}$LL${}^{\prime}$})}+\underbrace{\alpha^{2}_{s}\,f_{3}(\alpha_{s}L_{\mathrm{M}})}_{(\text{N${}^{3}$LL,N${}^{3}$LL${}^{\prime}$})}+\ldots\;\bigg{]}\\\
&\quad\times\bigg{\\{}1(\text{LL,NLL})+\alpha_{s}(\text{NLL${}^{\prime}$,N${}^{2}$LL})+\alpha_{s}^{2}(\text{N${}^{2}$LL${}^{\prime}$,N${}^{3}$LL})+\alpha_{s}^{3}(\text{N${}^{3}$LL${}^{\prime}$,N${}^{4}$LL})+\ldots\bigg{\\}}\,.\end{split}$
(1.1)

Therein, the desired precisions of the anomalous dimensions are specified
between the square brackets in the exponent as for a given logarithmic
accuracy, while the according requirements on the fixed-order elements are
presented within the curly brackets. In this work, we will evaluate and
compare the resummed $q_{\tau}$ spectra on the next-to-leading-logarithmic
(NLL), N2LL, and approximate N2LL′ (aN2LL′) levels.

The paper is structured as follows. In Sec. 2, we will utilise the strategy of
expansion by dynamic regions [69, 70, 71, 72] and effective field theories,
i.e. SCETII [80, 81, 82] and HQET [83, 84, 85, 86], to derive the
factorisation formula governing the leading asymptotic behaviour of ${{\rm
d}{\sigma}_{{t{\bar{t}}}}}/{{\rm d}q_{\tau}}$ in the limit $q_{\tau}\to 0$.
Then, the (rapidity) renormalisation group equations will be solved in Sec. 3
for the respective sectors participating into the factorisation formula, from
which we exponentiate the characteristic logarithmic constituents in the
impact space and thereby accomplish the resummation of the singular terms in
the momentum space. Sec. 4 will be devoted to the numeric evaluations on the
spectra $q_{\tau}$. Therein, we will at first validate the approximations of
our factorisation formula up to N2LO, and then present the resummation
improved differential distributions for three particular observables, ${\rm
d}\sigma_{{t{\bar{t}}}}/q_{\mathrm{T},\mathrm{in}}$, ${\rm
d}\sigma_{{t{\bar{t}}}}/q_{\mathrm{T},\mathrm{out}}$, and ${\rm
d}\sigma_{{t{\bar{t}}}}/\Delta\phi_{{t{\bar{t}}}}$. Finally, we will offer
some concluding remarks in Sec. 5.

##  2 Factorisation

###  2.1 Kinematics and the factorised cross section

Figure 1:  The kinematics on the transverse plane in the laboratory reference
frame. $\vec{P}_{t({\bar{t}})}^{\perp}$ stands for the transverse momentum of
the (anti-)top quark. $\Delta\Phi_{{t{\bar{t}}}}$ is the azimuthal opening
angle between the top and anti-top quarks. $\vec{\tau}$ is a unit reference
vector in the transverse plane.

We start this section with the elaboration on the kinematics. As illustrated
in Fig. 1, the main concern of this work is on the interplay between a
reference unit vector $\vec{\tau}$ and the transverse momentum
$\vec{q}_{\mathrm{T}}$ of the ${t{\bar{t}}}$ system. By means of the reference
vector $\vec{\tau}$, $\vec{q}_{\mathrm{T}}$ can be decomposed into two parts,
the projection component $\vec{q}_{\tau_{\|}}$ and the rejection one
$\vec{q}_{\tau_{\perp}}$, i.e.,

$\displaystyle\vec{q}_{\mathrm{T}}=\vec{q}_{\tau_{\perp}}+\vec{q}_{\tau_{\|}}=q_{\tau_{\perp}}\vec{\tau}\times\vec{n}+q_{\tau_{\|}}\vec{\tau}\,.$
(2.1)

where $\vec{n}$ is another unit vector pointing to one of beam directions in
the laboratory reference frame. In the numeric implementation presented in
this paper, the magnitude of the projection $\vec{q}_{\tau_{\|}}$ is of
primary interest, which will hereafter be dubbed
$q_{\tau}\equiv|\vec{q}_{\tau_{\|}}|$.

The fixed-order calculation on the $q_{\tau}$ spectrum can be realised using
the QCD factorisation theorem of [101], that is,

$\begin{split}\frac{{\rm d}^{5}{\sigma_{{t{\bar{t}}}}}}{{\rm
d}M_{t{\bar{t}}}^{2}\,{\rm d}^{2}\vec{P}_{t}^{\perp}\,{\rm
d}Y_{t{\bar{t}}}\,{\rm
d}q_{\tau}}\,=\;&\sum_{\mathrm{sign}[P_{t}^{z}]}\frac{1}{16s(2\pi)^{6}}\int{\rm
d}^{2}\vec{q}_{\mathrm{T}}\,\Theta_{\mathrm{kin}}\,\delta\Big{[}q_{\tau}-|\vec{q}_{\mathrm{T}}\cdot\vec{\tau}|\Big{]}\,\frac{\Sigma_{{t{\bar{t}}}}}{M_{\mathrm{T}}^{t{\bar{t}}}\,|P_{t}^{z}|}\,,\end{split}$
(2.2)

where $M_{t{\bar{t}}}$ denotes the invariant mass of the ${t{\bar{t}}}$
system, and $s$ is the colliding energy. In this work will concentrate on
$\sqrt{s}=13$ TeV throughout. $Y_{t{\bar{t}}}$ and
$M_{\mathrm{T}}^{t{\bar{t}}}$ are for the pseudorapidity and the transverse
mass of the ${t{\bar{t}}}$ pair in the laboratory frame (LF), respectively.
$P_{t}^{z}$ represents the longitudinal components of the top quark momentum
measured from the $z$-direction rest frame ($z$RF) of the ${t{\bar{t}}}$ pair.
The $z$RF can be obtained through boosting the LF along one of the beam
directions until the longitudinal momentum of the ${t{\bar{t}}}$ pair has been
eliminated.

To perform the integral of $\vec{q}_{\mathrm{T}}$ in Eq. (2.2), it is of
essence to establish suitable kinematical boundaries to fulfill energy-
momentum conservation condition. To this end, we introduce the function
$\Theta_{\mathrm{kin}}$ to impose the following constraints,

$\begin{split}\Theta_{\mathrm{kin}}\,=\,&\Theta\Big{[}\sqrt{s}-M_{\mathrm{T}}^{t{\bar{t}}}-|\vec{q}_{\mathrm{T}}|\Big{]}\,\Theta\Big{[}M_{\mathrm{T}}^{t{\bar{t}}}-m_{\mathrm{T}}^{t}-m_{\mathrm{T}}^{{\bar{t}}}\Big{]}\,\Theta\left[\sinh^{-1}\left(\sqrt{\frac{(M_{t{\bar{t}}}^{2}+s)^{2}}{4\,s\,{M_{\mathrm{T}}^{t{\bar{t}}}}^{2}}-1}\right)-|Y_{t{\bar{t}}}|\right]\,,\end{split}$
(2.3)

where $\Theta[\dots]$ is the usual Heaviside function. Therein,
$m_{\mathrm{T}}^{t}$ and $m_{\mathrm{T}}^{{\bar{t}}}$ are the transverse
masses of the top and anti-top quarks in the LF, respectively. Finally, Eq.
(2.2) also entails ${\Sigma_{{t{\bar{t}}}}}$, encoding the contributions from
all partonic processes, defined as,

$\begin{split}{\Sigma_{{t{\bar{t}}}}}\,=&\;\sum_{i,j}\int^{1}_{0}\frac{{\rm
d}x_{n}}{x_{n}}\,\frac{{\rm
d}x_{{\bar{n}}}}{x_{{\bar{n}}}}\,f_{i/N}(x_{n})\,f_{j/\bar{N}}(x_{{\bar{n}}})\,\sum_{r}\,\int\,\prod_{m}^{r}{\rm
d}\Phi_{k_{m}}\;(2\pi)^{4}\delta^{4}\left(p_{i}+p_{j}-p_{t}-p_{\bar{t}}-\sum_{m}k_{m}\right)\\\
&\;\times\overline{\sum_{\mathrm{hel,col}}}\left|\mathcal{M}(i+j\to
t+\bar{t}+X)\right|^{2}\,.\end{split}$ (2.4)

Here $f_{i/N}(x)$ is the parton distribution functions (PDFs) for the parton
$i$ within the proton $N$ carrying the momentum fraction $x$. ${\rm
d}\Phi_{k_{m}}$ characterises the phase space integral of the $m$-th emitted
parton, that is,

$\begin{split}{\rm d}\Phi_{k_{m}}\,\equiv\,\frac{{\rm
d}^{4}k_{m}}{(2\pi)^{3}}\,\delta(k_{m}^{2})\,\Theta(k_{m}^{0})\,=\,\frac{1}{2}\;\frac{{\rm
d}y_{m}}{2\pi}\;\frac{{\rm
d}^{2}\vec{k}_{m}^{\perp}}{(2\pi)^{2}}\,,\end{split}$ (2.5)

where $y_{m}$ and $\vec{k}_{m}^{\perp}$ indicate the rapidity and transverse
momentum of the occurring emission, respectively. $|\mathcal{M}|^{2}$ is the
squared transition amplitude for the partonic processes of the indices
indicated. Substituting Eq. (2.4) into Eq. (2.2), it is ready to perform the
fixed-order calculations of the spectra of $q_{\tau}$. In the vicinity of
$q_{\tau}=0$ GeV, however, this perturbative expansion fails to converge, and
an asymptotic expansion of ${\rm d}\sigma_{{t{\bar{t}}}}/{\rm d}q_{\tau}$ can
be carried out in the small parameter $\lambda_{\tau}\equiv
q_{\tau}/M_{t{\bar{t}}}$,

$\displaystyle\frac{{\rm d}{\sigma_{{t{\bar{t}}}}}}{{\rm
d}q_{\tau}}\;\sim\;\sigma^{{t{\bar{t}}}}_{\mathrm{B}}\,\sum_{m,n}\,\left[\frac{\alpha_{s}(M_{t{\bar{t}}})}{4\pi}\right]^{m}\,\left[\underbrace{c^{(0)}_{m,n}\,\frac{\ln^{n}(\lambda_{\tau})}{\lambda_{\tau}}}_{\mathrm{LP}}+\,\underbrace{c^{(1)}_{m,n}\,\ln^{n}(\lambda_{\tau})}_{\mathrm{NLP}}+\,\underbrace{c^{(2)}_{m,n}\,\lambda_{\tau}\ln^{n}(\lambda_{\tau})}_{\mathrm{N}^{2}\mathrm{LP}}+\dots\right]\;,$
(2.6)

indicating the leading, next-to-leading and next-to-next-to-leading power
terms in $\lambda_{\tau}$, labelled LP, NLP, and N2LP, respectively. Therein,
$\sigma^{{t{\bar{t}}}}_{\mathrm{B}}$ is the Born level total cross of the
process $pp\to{t{\bar{t}}}+X$, $\alpha_{s}$ denotes the strong coupling
constant, and $c^{(k)}_{m,n}$ represents the coefficient for the asymptotic
constituent with the superscript $k$ specifying the occurring power. Thus,
conventionally, the leading power terms
$c^{(0)}_{m,n}\ln^{n}(\lambda_{\tau})/\lambda_{\tau}$ are associated with the
most singular behaviors in the low $q_{\tau}$ domain and also the main concern
of this work. It is important to note, however, that also the next-to-leading
power terms, $c^{(1)}_{m,n}\ln^{n}(\lambda_{\tau})$, are divergent as
$\lambda_{\tau}\to 0$.

###  2.2 Dynamic regions

Based on the strategy of expansion of dynamic regions [69, 70, 71, 72], the
asymptotic series of Eq. (2.6) can be interpreted with the aid of a set of
regions from the phase space and loop integrals. This work, in particular,
will choose the formalism of [72]. We base the definition of our regions on
the criterion of domain completeness, i.e. the existence of a set of non-
intersecting dynamic regions that cover the whole integration domain. This
criterion plays an essential role in consistently extrapolating the expanded
integrands from their own convergent domains to the entire integration ranges.
Other constraints are also imposed therein, including the regularisation of
the expanded integrands and the (at least partial) commutativity amongst the
asymptotic expansions. The former case can be fulfilled by introducing the
rapidity regulator [102, 98, 100, 87, 88, 89, 90] in implementing the SCETII
formalism [80, 81, 82]. However, for the latter criterion, we assume that all
the non-commutative dynamic regions, such as the collinear-plane modes [72],
will cancel out in the eventual $q_{\tau}$ spectra. It merits noting that,
this ansatz, together with the proposal of [72], has been scrutinised only
within the one-loop integrals in the various kinematical limits. We regard
their effectiveness on ${t{\bar{t}}}$ hadroproduction as the primary
hypothesis in this work. Recent developments on the criteria to implement the
region analysis can be found in [103, 104, 105].

As the first application of the domain completeness, we explore the relevant
modes for the $\vec{q}_{\tau_{\perp}}$ integral here. Following from Eq.
(2.2), two types of scales will be involved, one of $\mathcal{O}(q_{\tau})$
and another of $\mathcal{O}(M_{t{\bar{t}}})$. In order to disentangle the
influences of those two scales and to also fulfil the constraints of [72], we
identify the dynamic regions for the rejection component as follows

$\displaystyle\mathbf{Isotropic}\textrm{-}\mathbf{recoil}:$
$\displaystyle\;M_{t{\bar{t}}}\gg|\vec{q}_{\tau_{\perp}}|\sim q_{\tau}\,,$
(2.7) $\displaystyle\mathbf{Asymmetric}\textrm{-}\mathbf{recoil}:$
$\displaystyle\;|\vec{q}_{\tau_{\perp}}|\sim M_{t{\bar{t}}}\gg q_{\tau}\,.$
(2.8)

To precisely separate both regimes and still cover the complete integration
range, we introduce the auxiliary boundary $\Lambda_{\tau}$ satisfying
$M_{t{\bar{t}}}\gg\Lambda_{\tau}\gg q_{\tau}$, from which the two non-
intersecting intervals read $|\vec{q}_{\tau_{\perp}}|\lesssim\Lambda_{\tau}$
and $\Lambda_{\tau}\lesssim|\vec{q}_{\tau_{\perp}}|$. However, as demonstrated
in [72], those auxiliary boundary dependences will all drop out after
assembling all relevant domains. Their dependence is thus dropped in the
following. An analogous analysis is also applied to the phase space integral
for the real emissions in Eq. (2.5), which leads to,

$\displaystyle\mathbf{Soft}:\;$ $\displaystyle\;{y_{m}}\sim
0\,,\quad\quad\quad\quad~{}|\vec{k}_{m}^{\perp}|\sim q_{\tau}\,,$ (2.9)
$\displaystyle\mathbf{Beam}\textrm{-}\mathbf{collinear}\textrm{-}n:\,$
$\displaystyle\;{y_{m}}\sim+\ln\lambda_{\tau}\,,\quad\quad|\vec{k}_{m}^{\perp}|\sim
q_{\tau}\,,$ (2.10)
$\displaystyle\mathbf{Beam}\textrm{-}\mathbf{collinear}\textrm{-}{\bar{n}}:\;$
$\displaystyle\;{y_{m}}\sim-\ln\lambda_{\tau}\,,\quad\quad|\vec{k}_{m}^{\perp}|\sim
q_{\tau}\,,$ (2.11)
$\displaystyle\mathbf{Jet}\textrm{-}\mathbf{collinear}\textrm{-}{n_{\mathrm{J}}}:\;$
$\displaystyle\;{y_{m}}\sim
0\,,\quad\quad\quad\quad~{}|\vec{k}_{m}^{\perp}|\sim M_{t{\bar{t}}}\,.$ (2.12)

Therein, the decomposition of the $\vec{k}_{m}^{\perp}$ integration range
respects a similar pattern as the $\vec{q}_{\tau_{\perp}}$ case. The rapidity
integrals are broken down with respect to the reference points
$\pm\Lambda_{y_{m}}$ with the relationship
$0\ll\Lambda_{y_{m}}\ll\ln\lambda_{\tau}$, from which three non-overlapping
regions emerge, ${y_{m}}\lesssim-\Lambda_{y_{m}}$,
$-\Lambda_{y_{m}}\lesssim{y_{m}}\lesssim\Lambda_{y_{m}}$, and
$\Lambda_{y_{m}}\lesssim{y_{m}}$. As before, the dependence on the auxiliary
boundary $\Lambda_{y_{m}}$ cancels, and we omit it in the following. In
addition, please note that in deriving Eq. (2.12) the super-hard-collinear
domains ${y_{m}}\sim\pm\ln\lambda_{\tau}$ and $|\vec{k}_{m}^{\perp}|\sim
M_{t{\bar{t}}}$ are ignored as they explicitly contradict the energy-momentum
conservation condition in Eq. (2.4).

Apart from the above momentum modes, the amplitudes in Eq. (2.4) also contain
the loop integrals. A region analysis of this case can proceed in principle in
a similar way as for the real emission corrections. However, due to the
variability in the offshellness and the varying multiplicity of the external
particles, exhausting all relevant circumstances herein is usually more
challenging. To this end, this work will utilise momentum regions that can
lead to pinched singularities (PS) in a hadron collider process [67, 68, 66,
65, 106, 107], which in general consists of hard, collinear, soft, Coulomb,
and Glauber regions. During our calculation, to circumvent the complexity
induced by Coulomb singularities, we introduce the lower cutoff upon
$M_{t{\bar{t}}}$ in the phase space integral to stay clear of the
${t{\bar{t}}}$ production threshold. In order to cope with the remaining hard,
collinear, and soft components, we apply effective field theories, i.e. SCETII
[80, 81, 82] and HQET [83, 84, 85, 86], onto the massless and heavy partons,
respectively. Regarding Glauber gluon exchanges, even though their
cancellation for inclusive observables has been demonstrated in colour-singlet
production in hadronic collisions using a variety of approaches [108, 109,
110, 111], a systematic discussion of their effects on the production of
coloured systems, like hadronic ${t{\bar{t}}}$ pair production, is still
absent. While leftover soft contributions were observed in ${t{\bar{t}}}$
production in the context of light-cone ordered perturbation theory [112, 113,
114], it nevertheless deserves further confirmation from perturbative QCD, in
particular as to whether those soft remnants are eikonalisable or not [115,
116]. In this work, we will follow the approach of [57, 58], and assume the
irrelevance of the Glauber contributions. Recent developments in generalizing
SCET to encode the Glauber interactions can be found in [117, 118, 119, 115,
120].

We are now ready to summarise the dynamic modes presiding over the virtual and
real corrections to the $q_{\tau}$ spectra,

$\displaystyle\mathbf{Hard}:$ $\displaystyle\quad
k_{h}~{}\sim\;M_{t{\bar{t}}}\,\big{[}\mathcal{O}(1),\mathcal{O}(1),\mathcal{O}(1)\big{]}_{n}\,,$
(2.13) $\displaystyle\mathbf{Soft}:$ $\displaystyle\quad
k_{s}~{}\sim\;M_{t{\bar{t}}}\,\big{[}\mathcal{O}(\lambda_{\tau}),\mathcal{O}(\lambda_{\tau}),\mathcal{O}(\lambda_{\tau})\big{]}_{n}\,,$
(2.14) $\displaystyle\mathbf{Beam}\textrm{-}\mathbf{collinear}\textrm{-}{n}:$
$\displaystyle\quad
k_{c}~{}\sim\;M_{t{\bar{t}}}\,\big{[}\mathcal{O}(1),\mathcal{O}(\lambda_{\tau}^{2}),\mathcal{O}(\lambda_{\tau})\big{]}_{n}\,,$
(2.15)
$\displaystyle\mathbf{Beam}\textrm{-}\mathbf{collinear}\textrm{-}{{\bar{n}}}:$
$\displaystyle\quad
k_{\bar{c}}~{}\sim\;M_{t{\bar{t}}}\,\big{[}\mathcal{O}(\lambda_{\tau}^{2}),\mathcal{O}(1),\mathcal{O}(\lambda_{\tau})\big{]}_{n}\,,$
(2.16)
$\displaystyle\mathbf{Jet}\textrm{-}\mathbf{collinear}\textrm{-}{n_{\mathrm{J}}}:$
$\displaystyle\quad
k_{\mathrm{J}}~{}\sim\;M_{t{\bar{t}}}\,\big{[}\mathcal{O}(1),\mathcal{O}(\lambda_{\tau}^{2}),\mathcal{O}(\lambda_{\tau})\big{]}_{n_{\mathrm{J}}}\,.$
(2.17)

In writing those momentum modes, the light-cone coordinate system has been
applied, from which an arbitrary momentum $p^{\mu}$ are decomposed and
reforged as

$p^{\mu}=\frac{p\cdot{\bar{n}}_{k}}{2}n_{k}^{\mu}+\frac{p\cdot
n_{k}}{2}{\bar{n}}_{k}^{\mu}+{p}^{\perp,\mu}_{n_{k}}\equiv\frac{p_{{\bar{n}}_{k}}}{2}n_{k}^{\mu}+\frac{p_{n_{k}}}{2}{\bar{n}}_{k}^{\mu}+{p}^{\perp,\mu}_{n_{k}}\equiv\big{[}p_{{\bar{n}}_{k}},p_{n_{k}},{p}_{n_{k}}^{\perp}\big{]}_{n_{k}}.$
(2.18)

Therein, $n_{k}$ and ${\bar{n}}_{k}$ are two light-like vectors with the
normalisation conditions $n_{k}\cdot{\bar{n}}_{k}=2$ and
$n^{0}_{k}={\bar{n}}^{0}_{k}=1$. Throughout this paper, we will make use of
the symbols $n$ (${\bar{n}}$) and $n_{\mathrm{J}}$ to characterise the
positive (negative) beam and jet directions, respectively. In accordance with
Eq. (2.12), at least one of the transverse components of the vector
$n^{\mu}_{\mathrm{J}}$ should be non-vanishing. In establishing Eqs.
(2.14)-(2.17), the offshellnesses of the collinear and soft fluctuations are
assigned to be of $\mathcal{O}(\lambda^{2}_{\tau})$ so as to accommodate the
projected transverse momentum. The hard-collinear degree of freedom, e.g.
$k_{hc}\sim
M_{t{\bar{t}}}\,[\mathcal{O}(\lambda_{\tau}),\mathcal{O}(1),\mathcal{O}(\sqrt{\lambda_{\tau}})]_{n}$,
is neglected in the analysis as it always results in a scaleless integral for
the leading power accuracy investigated in this paper.

By comparison with the calculations on ${\rm d}\sigma_{{t{\bar{t}}}}/{\rm
d}q_{\mathrm{T}}$ of [55, 56, 57, 121, 58, 59], in addition to the common
pattern of the hard, beam-collinear, and soft regions, this work will also
take into account the central jet mode as shown in Eq. (2.17). For the low
$q_{\mathrm{T}}$ domain, the subleading nature of the jet region has been
demonstrated in [65, 122] by analyzing the powers of the relevant PS surfaces.
This conclusion is also expected to hold for our symmetric configuration in
eq. (2.7). Nevertheless, aside from the isotropic recoil, the region expansion
strategy in Ref. [72] also prompts the asymmetric posture of eq. (2.8) to
participate into the $q_{\tau}$ spectra. Owing to the kinematical differences
between the two, it is, however, not straightforward to apply the conclusions
of [65, 122] to the contributions induced by Eq. (2.8). In view of this, we
will revisit the power rules of all the relevant configurations conducted by
Eqs. (2.7)-(2.8) and Eqs. (2.13)-(2.17) in the following subsections.

For the sake of the later discussion, we categorise the $q_{\tau}$ spectrum
according to the number of the embedded jet modes $N_{\mathrm{J}}$, more
explicitly,

$\frac{{\rm d}^{5}{\sigma_{{t{\bar{t}}}}}}{{\rm d}M_{t{\bar{t}}}^{2}\,{\rm
d}^{2}\vec{P}_{t}^{\perp}\,{\rm d}Y_{t{\bar{t}}}\,{\rm
d}q_{\tau}}\longrightarrow\,\sum_{m=0}^{\infty}\,\frac{{\rm
d}^{5}{\sigma_{{t{\bar{t}}}}}}{{\rm d}M_{t{\bar{t}}}^{2}\,{\rm
d}^{2}\vec{P}_{t}^{\perp}\,{\rm d}Y_{t{\bar{t}}}\,{\rm
d}q_{\tau}}\Bigg{|}_{N_{\mathrm{J}}=m}\,.$ (2.19)

We will first of all concentrate on the $N_{\mathrm{J}}=0$ configuration in
Sec. 2.3, and then, in Sec. 2.4, we move on to the more general situation, the
$N_{\mathrm{J}}\geq 1$ contributions.

###  2.3 The case of $N_{\mathrm{J}}=0$

Figure 2:  Dynamic regions in the $N_{\mathrm{J}}=0$ configuration. Herein,
the black bubbles stand for the hard modes. The red bubble and springs encode
the soft fluctuations connecting the heavy partons represented by the black
double lines as well as the collinear modes depicted in the blue straight
lines.

In this part, we will discuss the dynamic regions contributing to the
$N_{\mathrm{J}}=0$ configuration, which comprises the hard, beam-collinear,
and soft regions, as exhibited in Eqs. (2.13)-(2.16). As opposed to the beam-
collinear and soft regions, which can be assigned to both the phase space and
loop integrals, the hard mode is only present in the virtual processes. In
this regard, the transverse momentum for the ${t{\bar{t}}}$ system observes
only the isotropic pattern in Eq. (2.7), which thus permits us to perform the
expansion on the kinematic variables in Eq. (2.2) as follows,

$\displaystyle P^{z}_{t}\Bigg{|}_{N_{\mathrm{J}}=0}$
$\displaystyle=\;\frac{x_{t}\beta_{t{\bar{t}}}}{\sqrt{1-\beta^{2}_{{t{\bar{t}}}}}}\,m_{t}\,+\mathcal{O}(\lambda_{\tau})\,,$
(2.20)

where $\beta_{t{\bar{t}}}\equiv\sqrt{1-4m_{t}^{2}/M_{t{\bar{t}}}^{2}}$ with
$m_{t}$ denoting the mass of the (anti-)top quark. $x_{t}$ is the cosine of
the scattering angle of the top quark in the ${t{\bar{t}}}$ rest frame.
Further,

$\displaystyle M_{\mathrm{T}}^{t{\bar{t}}}\Bigg{|}_{N_{\mathrm{J}}=0}$
$\displaystyle=\;M_{t{\bar{t}}}\,+\mathcal{O}(\lambda_{\tau})\,,$ (2.21)

and, thus,

$\displaystyle\Theta_{\mathrm{kin}}\Bigg{|}_{N_{\mathrm{J}}=0}$
$\displaystyle=\;\Theta\left(\sqrt{s}-M_{t{\bar{t}}}^{2}\right)\,\Theta\left(M_{t{\bar{t}}}-2\sqrt{m_{t}^{2}+(\vec{P}^{\perp}_{t})^{2}}\right)\,\Theta\left[\sinh^{-1}\left(\frac{s-M_{t{\bar{t}}}^{2}}{2M_{t{\bar{t}}}\sqrt{s}}\right)-|Y_{t{\bar{t}}}|\right]+\mathcal{O}(\lambda_{\tau})$
$\displaystyle\equiv\;\Theta_{\mathrm{kin}}^{(0)}+\mathcal{O}(\lambda_{\tau})\,.$
(2.22)

Please note, in the results above only the LP contributions are kept.
Similarly, in the following, we define the LP term in Eq. (2.3) as
$\Theta_{\mathrm{kin}}^{(0)}$ hereafter.

The remaining task is now to expand the partonic convolution function
$\Sigma_{{t{\bar{t}}}}$ in line with Eqs. (2.13)-(2.16). This can be achieved
by means of the effective field theories, i.e. SCETII and HQET. Therein, the
beam-collinear regions are embodied in terms of the gluon and quark fields,
$A_{n(\bar{n})}$ and $\xi_{n(\bar{n})}$ [79, 123, 81, 82] in SCETII. The soft
corrections on the heavy (anti-)quark are encoded by the field
$h_{v}(\chi_{v})$ from HQET, while those from the incoming gluons and massless
quarks are reflected by the fields $A_{s}$ and $q_{s}$ [81, 82] in SCETII. The
hard mode can be taken care of by the effective Hamiltonian
$\mathcal{H}_{{t{\bar{t}}}}^{\mathrm{eff}}$ [42], which consists of products
of Wilson coefficients and the corresponding field operators.

Up to LP, the interactions in the individual momentum regions are governed by
the effective Lagrangians [80, 81, 82, 86],

$\displaystyle\mathcal{L}_{n}=$ $\displaystyle\;\bar{\xi}_{n}\Big{(}in\cdot
D_{n}+i\not{D}_{n\perp}\frac{1}{i{\bar{n}}\cdot
D_{n}}i\not{D}_{n\perp}\Big{)}\frac{\not{{\bar{n}}}}{2}\xi_{n}-\frac{1}{2}\mathbf{Tr}\Big{[}F^{\mu\nu}_{n}F_{\mu\nu}^{n}\Big{]}\,,$
(2.23) $\displaystyle\mathcal{L}_{s}=$
$\displaystyle\;-\frac{1}{2}\mathbf{Tr}\Big{[}F^{\mu\nu}_{\text{s}}F_{\mu\nu}^{\text{s}}\Big{]}+\bar{q}_{\text{s}}i\not{D}_{\text{s}}q_{\text{s}}\,,$
(2.24) $\displaystyle\mathcal{L}_{v}=$
$\displaystyle\;\bar{h}_{v}\left(i\partial\cdot
v\right)h_{v}+\bar{\chi}_{v}\left(i\partial\cdot v\right)\chi_{v}\,,$ (2.25)

where $D^{\mu}_{n(s)}\equiv\partial^{\mu}-ig_{s}A^{\mu}_{n(s)}$ and
$F_{n(s)}^{\mu\nu}\equiv(i/g_{s})[D^{\mu}_{n(s)},D^{\nu}_{n(s)}]$ stand for
the covariant derivative and the field strength tensor for the collinear
(soft) fields $A^{\mu}_{n(s)}$, respectively. The Lagrangian
$\mathcal{L}_{{\bar{n}}}$ can be obtained by exchanging
$n\leftrightarrow{\bar{n}}$ in $\mathcal{L}_{n}$. In deriving Eqs.
(2.23)-(2.25), the decoupling transformations [73, 49] have already been
performed to strip the soft particle of the collinear and heavy partons. As a
result, the (anti-)top quark is free of any interactions at this point, while
the collinear partons only communicate with themselves. The LP contribution of
$\Sigma_{{t{\bar{t}}}}$ can be found through assembling the amplitudes induced
by $\mathcal{H}_{{t{\bar{t}}}}^{\mathrm{eff}}$ and the effective Lagrangians
in Eqs. (2.23)-(2.25), and we follow the scheme of [42, 56] to construct the
hard contributions. In light of the rapidity divergences arising from the soft
and collinear integrals, the exponential regulator proposed in [89, 90] is
applied throughout. Finally, the decoupling nature of the Lagrangians in Eqs.
(2.23)-(2.25) allows us to rewrite $\Sigma_{{t{\bar{t}}}}$ as,

$\displaystyle\Sigma_{{t{\bar{t}}}}\Bigg{|}_{N_{\mathrm{J}}=0}=$
$\displaystyle\;\frac{8\pi^{2}}{M_{t{\bar{t}}}^{2}}\,\sum_{\kappa}\,\int{{\rm
d}^{2}\vec{b}_{\mathrm{T}}}\,\exp\left(\mathrm{i}\,\vec{b}_{\mathrm{T}}\cdot\vec{q}_{\mathrm{T}}\right)\,\widetilde{\Sigma}_{{t{\bar{t}}}}^{[\kappa]}(\vec{b}_{\mathrm{T}},M_{t{\bar{t}}},\beta_{t{\bar{t}}},x_{t},Y_{t{\bar{t}}},\mu,\nu)\,\mathcal{W}_{t}\,\mathcal{W}_{{\bar{t}}}\,,$
(2.26)

where $\vec{b}_{\mathrm{T}}$ is the impact parameter introduced during the
Fourier transformation. $\Sigma_{{t{\bar{t}}}}$ includes the contributions
from all channels
$\kappa\in\\{g_{n}g_{{\bar{n}}},q^{i}_{n}\bar{q}^{j}_{{\bar{n}}},q^{i}_{{\bar{n}}}\bar{q}^{j}_{n}\\}$,
with $i,j\in\\{u,d,c,s,b\\}$ indicating the flavour of the quark field. $\mu$
and $\nu$ are the scales associated with the virtuality and rapidity
renormalisations, respectively. Finally, $\mathcal{W}_{t}$ and
$\mathcal{W}_{{\bar{t}}}$ are the heavy parton correlation functions,

$\displaystyle\mathcal{W}_{t}=$
$\displaystyle\frac{1}{4m_{t}N_{c}}\,\mathbf{Tr}\,\langle
0|{h}_{v_{t}}\left(0\right)|t\rangle\langle
t|\bar{h}_{v_{t}}\left(0\right)\frac{1+\not{v}_{t}}{2}|0\rangle=1\,,$ (2.27)
$\displaystyle\mathcal{W}_{{\bar{t}}}=$
$\displaystyle-\frac{1}{4m_{t}N_{c}}\,\mathbf{Tr}\,\langle
0|\bar{\chi}_{v_{\bar{t}}}\left(0\right)|{\bar{t}}\rangle\langle{\bar{t}}|\frac{1-\not{v}_{{\bar{t}}}}{2}{\chi}_{v_{{\bar{t}}}}\left(0\right)|0\rangle=1\,,$
(2.28)

where $N_{c}=3$ is the colour factor. $v_{t({\bar{t}})}$ is the velocity of
the (anti)top quark in the rest frame of the ${t{\bar{t}}}$ system.
Considering that the (anti-)top quark up to the LP accuracy amounts to a free
particle, the correlation functions $\mathcal{W}_{t({\bar{t}})}$ will never
receive any perturbative corrections. Therefore, in the second steps of Eqs.
(2.27)-(2.28), we evaluate them using the tree-level expressions. It should be
stressed that this result is subject to the absence of Coulomb interactions.
Otherwise, incorporating Coulomb potentials can manifestly break the
$\mathcal{W}_{t}$ and $\mathcal{W}_{{\bar{t}}}$ factorisation, thereby
inducing non-trivial perturbative corrections [48, 49, 50, 52, 53, 54]. To
this end, we will constrain the investigation in this work within the domain
$M_{t{\bar{t}}}\geq 400$ GeV, or, equivalently, $\beta_{t{\bar{t}}}\gtrsim
0.5$.111Including the Coulomb singularity in the low $q_{\tau}$ regime is
intricate since its super-renormalisable nature [124, 125] may induce
additional contributions from the power suppressed terms in SCETII. A
discussion of the similar scenario, i.e. the interplay between SCETI and
pNRQCD, can be found in [48, 126].

The partonic cross sections $\widetilde{\Sigma}_{{t{\bar{t}}}}^{[\kappa]}$ are
built by suitably combining beam-collinear, soft and hard functions,

$\displaystyle\widetilde{\Sigma}_{{t{\bar{t}}}}^{[q^{i}_{n}\bar{q}^{j}_{{\bar{n}}}]}(\vec{b}_{\mathrm{T}},M_{t{\bar{t}}},\beta_{t{\bar{t}}},x_{t},Y_{t{\bar{t}}},\mu,\nu)=$
$\displaystyle\left(\frac{1}{2N_{c}}\right)^{2}\,\mathcal{B}_{n}^{[q_{n}^{i}]}(\eta_{n},b_{\mathrm{T}},\mu,\nu)\,\mathcal{B}_{{\bar{n}}}^{[\bar{q}_{n}^{j}]}(\eta_{{\bar{n}}},b_{\mathrm{T}},\mu,\nu)\,$
$\displaystyle\times\,\sum_{\alpha,\beta}\left\\{\mathcal{H}_{\alpha\beta}^{[q^{i}_{n}\bar{q}^{j}_{{\bar{n}}}]}(M_{t{\bar{t}}},\beta_{t{\bar{t}}},x_{t},\mu)\,\mathcal{S}^{\alpha\beta}_{[q_{n}\bar{q}_{{\bar{n}}}]}(\vec{b}_{\mathrm{T}},\mu,\nu)\right\\}\,,$
(2.29)
$\displaystyle\widetilde{\Sigma}_{{t{\bar{t}}}}^{[g_{n}g_{{\bar{n}}}]}(\vec{b}_{\mathrm{T}},M_{t{\bar{t}}},\beta_{t{\bar{t}}},x_{t},Y_{t{\bar{t}}},\mu,\nu)=$
$\displaystyle\left(\frac{1}{N^{2}_{c}-1}\right)^{2}\,\sum_{\begin{subarray}{c}\alpha,\beta,{{h_{n}^{\phantom{\prime}}}},{{h_{n}^{\prime}}},{{h_{\bar{n}}^{\phantom{\prime}}}},{{h_{\bar{n}}^{\prime}}}\end{subarray}}\left\\{\mathcal{B}^{[g_{{\bar{n}}}]}_{{\bar{n}},{{h_{\bar{n}}^{\prime}}}{{h_{\bar{n}}^{\phantom{\prime}}}}}(\eta_{{\bar{n}}},\vec{b}_{\mathrm{T}},\mu,\nu)\,\mathcal{B}^{[g_{n}]}_{n,{{h_{n}^{\prime}}}{{h_{n}^{\phantom{\prime}}}}}(\eta_{n},\vec{b}_{\mathrm{T}},\mu,\nu)\,\right.$
$\displaystyle\left.\times\,\mathcal{H}_{\alpha\beta;{{{h_{\bar{n}}^{\prime}}}}{{h_{\bar{n}}^{\phantom{\prime}}}};{{h_{n}^{\prime}}}{{h_{n}^{\phantom{\prime}}}}}^{[g_{n}g_{{\bar{n}}}]}(M_{t{\bar{t}}},\beta_{t{\bar{t}}},x_{t},\mu)\,\mathcal{S}^{\alpha\beta}_{[g_{n}g_{{\bar{n}}}]}(\vec{b}_{\mathrm{T}},\mu,\nu)\right\\}\,,$
(2.30)

where the momentum fractions $\eta_{n}$ and $\eta_{{\bar{n}}}$ are defined as
$\eta_{n}=M_{t{\bar{t}}}\,e^{Y_{t{\bar{t}}}}/\sqrt{s}$ and
$\eta_{{\bar{n}}}=M_{t{\bar{t}}}\,e^{-Y_{t{\bar{t}}}}/\sqrt{s}$, respectively.
The expression for
$\widetilde{\Sigma}_{{t{\bar{t}}}}^{[q^{i}_{{\bar{n}}}\bar{q}^{j}_{{n}}]}$ can
be derived by exchanging the labels $n\leftrightarrow{\bar{n}}$.

$\mathcal{H}_{\alpha\beta}^{[q^{i}_{n}\bar{q}^{j}_{{\bar{n}}}]}$ and
$\mathcal{H}_{\alpha\beta;{{h_{\bar{n}}^{\prime}}}{{h_{\bar{n}}^{\phantom{\prime}}}};{{h_{n}^{\prime}}}{{h_{n}^{\phantom{\prime}}}}}^{[g_{n}g_{{\bar{n}}}]}$
are the hard functions for the quark and gluon initiated processes,
respectively. The indices $\alpha$ and $\beta$ label the colour states to
track the full colour correlation between the hard and the soft functions
detailed below. Similarly, the ${{h_{n}^{\phantom{\prime}}}}$,
${{h_{n}^{\prime}}}$, ${{h_{\bar{n}}^{\phantom{\prime}}}}$, and
${{h_{\bar{n}}^{\prime}}}$ denote the polarisation states of the incoming
gluons to capture all off-diagonal correlation effects of the beam-collinear
and hard functions. Please note, that the beam-hard function correlation for
massless quarks is devoid of off-diagonal contributions. The hard functions
now account for the LP contributions from the deep off-shell region in Eq.
(2.13). Their expressions read [42, 56],

$\displaystyle\mathcal{H}_{\alpha\beta}^{[q^{i}_{n}\bar{q}^{j}_{{\bar{n}}}]}(M_{t{\bar{t}}},\beta_{t{\bar{t}}},x_{t},\mu)=$
$\displaystyle\sum_{\alpha^{\prime},\beta^{\prime}}\,\Big{(}\mathcal{Z}^{h,{\alpha^{\prime}\alpha}}_{[{q_{n}\bar{q}_{{\bar{n}}}}]}\Big{)}^{*}\,\mathcal{Z}^{h,{\beta^{\prime}\beta}}_{[{q_{n}\bar{q}_{{\bar{n}}}}]}\,{\sum_{\\{a,a^{\prime},h\\}}}\,c_{a^{\prime}_{1}a^{\prime}_{2}a^{\prime}_{3}a^{\prime}_{4}}^{qq,(\alpha^{\prime})}\,\left[c_{a_{1}a_{2}a_{3}a_{4}}^{qq,(\beta^{\prime})}\right]^{*}\,$
$\displaystyle\quad\times\mathcal{M}^{*}(\bar{q}^{a^{\prime}_{1}h_{1}^{\phantom{\prime}}}_{{\bar{n}},j}{q}^{a^{\prime}_{2}h_{2}^{\phantom{\prime}}}_{n,i}\to
t_{a^{\prime}_{3}h_{3}^{\phantom{\prime}}}\bar{t}_{a^{\prime}_{4}h_{4}^{\phantom{\prime}}})\,\mathcal{M}(\bar{q}^{a_{1}h_{1}}_{{\bar{n}},j}{q}^{a_{2}h_{2}}_{n,i}\to
t_{a_{3}h_{3}}\bar{t}_{a_{4}h_{4}})\,,$ (2.31)
$\displaystyle\mathcal{H}_{\alpha\beta;h^{\prime}_{1}h_{1};h^{\prime}_{2}h_{2}}^{[g_{n}g_{{\bar{n}}}]}(M_{t{\bar{t}}},\beta_{t{\bar{t}}},x_{t},\mu)=$
$\displaystyle\sum_{\alpha^{\prime},\beta^{\prime}}\,\Big{(}\mathcal{Z}^{h,{\alpha^{\prime}\alpha}}_{[{g_{n}g_{{\bar{n}}}}]}\Big{)}^{*}\,\mathcal{Z}^{h,{\beta^{\prime}\beta}}_{[{g_{n}g_{{\bar{n}}}}]}\,{\sum_{\\{a,a^{\prime}\\},h_{3},h_{4}}}\,c_{a^{\prime}_{1}a^{\prime}_{2}a^{\prime}_{3}a^{\prime}_{4}}^{gg,(\alpha^{\prime})}\,\left[c_{a_{1}a_{2}a_{3}a_{4}}^{gg,(\beta^{\prime})}\right]^{*}\,$
$\displaystyle\quad\times\mathcal{M}^{*}(g^{a^{\prime}_{1}h^{\prime}_{1}}_{{\bar{n}}}g^{a^{\prime}_{2}h^{\prime}_{2}}_{n}\to
t_{a^{\prime}_{3}h_{3}^{\phantom{\prime}}}\bar{t}_{a^{\prime}_{4}h_{4}^{\phantom{\prime}}})\,\mathcal{M}(g^{a_{1}h_{1}}_{{\bar{n}}}g^{a_{2}h_{2}}_{n}\to
t_{a_{3}h_{3}}\bar{t}_{a_{4}h_{4}})$ (2.32)

where $\mathcal{M}$ denotes the amplitude for the partonic process
$\bar{q}q\to{t{\bar{t}}}$ or $gg\to{t{\bar{t}}}$. Therein, the variables
$a_{i}$$(i=1,\dots,4)$ are introduced for the colour states of the individual
external particles. In particular, every $a_{i}$ runs over the set
$\\{1,\dots,3\\}$ for quarks and anti-quarks, and $\\{1,\dots,8\\}$ for
gluons. Also, to facilitate our calculation, the orthonormal colour bases
$c^{qq}_{\\{a\\}}$ and $c^{gg}_{\\{a\\}}$ [127] are exploited in Eqs.
(2.31)-(2.32), more explicitly,

$\displaystyle c^{qq,(1)}_{a_{1}a_{2}a_{3}a_{4}}$
$\displaystyle=\frac{1}{3}\,\delta_{a_{1}a_{2}}\delta_{a_{3}a_{4}}\,,\;\;c^{qq,(2)}_{a_{1}a_{2}a_{3}a_{4}}=\frac{1}{\sqrt{2}}\sum_{c}T^{c}_{a_{1}a_{2}}T^{c}_{a_{3}a_{4}}\,,$
(2.33) $\displaystyle c^{gg,(1)}_{a_{1}a_{2}a_{3}a_{4}}$
$\displaystyle=\frac{1}{2\sqrt{6}}\,\delta_{a_{1}a_{2}}\delta_{a_{3}a_{4}}\,,\;\;c^{gg,(2)}_{a_{1}a_{2}a_{3}a_{4}}=\frac{i}{2\sqrt{3}}\,\sum_{c}f^{a_{1}ca_{2}}T^{c}_{a_{3}a_{4}}\,,\;\;c^{gg,(3)}_{a_{1}a_{2}a_{3}a_{4}}=\frac{1}{2}\sqrt{\frac{3}{5}}\sum_{c}d^{a_{1}ca_{2}}T^{c}_{a_{3}a_{4}}\,,$
(2.34)

where $T^{c}_{ab}$ stands for the generator in the fundamental representation
of the SU$(3)$ group. $f_{abc}$ and $d_{abc}$ mark the antisymmetric and
symmetric structure constants for the SU$(3)$ group, respectively. As alluded
to above, in calculating the hard functions, due to the absence of spin-
correlations for external quarks, we take the sum over all the helicity
configurations for the quark channel, while, to capture the non-diagonal gluon
polarisation effects, the explicit dependence of
$\mathcal{H}_{\alpha\beta;h^{\prime}_{1}h_{1};h^{\prime}_{2}h_{2}}^{[g_{n}g_{{\bar{n}}}]}$
on the gluon helicities is retained.

To evaluate the amplitudes in Eqs. (2.31)-(2.32), we make use of the on-shell
prescription to renormalise the top quark mass and the $n_{f}=5$ active flavor
scheme to handle the UV divergences from the strong coupling. The remaining
singularities are of infrared origin, and are subtracted by means of
$\mathcal{Z}^{h}_{[{g_{n}g_{{\bar{n}}}}]}$ and
$\mathcal{Z}^{h}_{[{q_{n}\bar{q}_{{\bar{n}}}}]}$ following the method in
[128]. Up to NLO, the amplitudes of all the helicity and colour configurations
can be extracted from R  ECOLA [129, 130]. The NNLO results are more involved,
and in consequence we will only address their logarithmic parts in this work.
For reference, grid-based NNLO results can be found in [131], and the progress
towards the full analytic evaluations is discussed in [132, 133, 134].

$\mathcal{B}_{n}^{[q_{n}^{i}]}$ and $\mathcal{B}^{[g_{n}]}_{n,h^{\prime}h}$
are the quark and gluon beam functions, respectively, collecting the
contributions from the region in Eq. (2.15). Their definitions in the
exponential regularisation scheme are [89, 135, 136],

$\displaystyle\mathcal{B}_{n}^{[q_{n}^{i}]}(\eta_{n},b_{\mathrm{T}},\mu,\nu)=$
$\displaystyle\;\mathcal{Z}^{c}_{[q_{n}]}\Big{(}\mathcal{Z}^{cs}_{[q_{n}]}\Big{)}^{-1}\,\lim_{\delta\to
0^{+}}\,\int\frac{{\rm
d}b_{n}}{4\pi}\,\exp\left(\frac{-\mathrm{i}\eta_{n}b_{n}P_{{\bar{n}}}}{2}\right)\,$
(2.35) $\displaystyle\times\mathbf{Tr}\,\langle
N(P)|\overline{\xi}^{i}_{n}W_{n}(-{ib_{0}\delta},-{ib_{0}\delta}+b_{n},\vec{b}_{\mathrm{T}})\frac{\not{{\bar{n}}}}{2}W^{\dagger}_{n}{\xi}^{i}_{n}(0)|N(P)\rangle\bigg{|}_{\delta=\frac{1}{\nu}}\,,$
$\displaystyle\mathcal{B}^{[g_{n}]}_{n,h^{\prime}h}(\eta_{n},\vec{b}_{\mathrm{T}},\mu,\nu)=$
$\displaystyle\;\mathcal{Z}^{c}_{[g_{n}]}\Big{(}{\mathcal{Z}^{cs}_{[g_{n}]}}\Big{)}^{-1}\,\lim_{\delta\to
0^{+}}\,\int\frac{{\rm
d}b_{n}}{4\pi}\,\exp\left(\frac{-\mathrm{i}\eta_{n}b_{n}P_{{\bar{n}}}}{2}\right)\,$
(2.36) $\displaystyle\times\mathbf{Tr}\,\langle
N(P)|\mathcal{A}^{\perp}_{n,\rho}(-{ib_{0}\delta},-{ib_{0}\delta}+b_{n},\vec{b}_{\mathrm{T}})\,\epsilon^{\rho}_{n,h^{\prime}}\,\left(-\eta_{n}P_{{\bar{n}}}\right)\,\epsilon^{\sigma*}_{n,h}\,\mathcal{A}^{\perp}_{n,\sigma}(0)|N(P)\rangle\bigg{|}_{\delta=\frac{1}{\nu}}\,,$

where $\mathcal{Z}^{c}_{[q_{n}]([g_{n}])}$ is the renormalisation constant for
the quark (gluon) beam function in the $\overline{\mathrm{MS}}$ scheme.
$\mathcal{Z}^{cs}_{[q_{n}]}$ and $\mathcal{Z}^{cs}_{[g_{n}]}$ represent the
ensuing zero-bin subtrahend to remove the soft-collinear overlapping terms.
$\delta$ is the exponential regulator suggested in [89], accompanied by the
constant $b_{0}\equiv 2\exp(-\gamma_{\mathrm{E}})$. $\mu$ and $\nu$ are the
scales associated with the virtuality and rapidity renormalisations,
respectively. Within the matrix elements, ${\xi}^{i}_{n}$ denotes the
collinear quark field of the flavour $i$ given in Eq. (2.23).
$\mathcal{A}_{n,\perp}^{\mu}\equiv\frac{1}{g_{s}}\,W_{n}^{\dagger}(iD_{n,\perp}^{\mu}W_{n})$
signifies the gauge invariant building block for the gluon field with $W_{n}$
denoting the collinear Wilson line [123]. Finally, $P$ is the momentum carried
by the initial proton with $P_{{\bar{n}}}=P\cdot{\bar{n}}$ being the largest
light-cone component. The anti-quark beam function
$\mathcal{B}_{n}^{[\bar{q}_{i}]}$ and those for the ${\bar{n}}$ direction can
be obtained by adjusting the labels and fields in Eqs. (2.35)-(2.36)
appropriately.

Comparing with the quark beam function, the gluon case possesses extra
indices, ${{h_{n}^{\phantom{\prime}}}},{{h_{n}^{\prime}}}\in\\{+,-\\}$, to
characterise its intrinsic polarisation effects [64]. In this work, the
following helicity basis is adopted,

$\displaystyle\epsilon^{\mu}_{n,\pm}\equiv\left\\{0,\frac{\mp
1}{\sqrt{2}},\frac{-\mathrm{i}}{\sqrt{2}},0\right\\}\,,\qquad\epsilon^{\mu}_{{\bar{n}},\pm}\equiv\left\\{0,\frac{\pm
1}{\sqrt{2}},\frac{-\mathrm{i}}{\sqrt{2}},0\right\\}\,.$ (2.37)

In principle, the representations of the helicity polarisation states
$\epsilon^{\mu}_{n,\pm}$ are not unique. Nevertheless, in order to avoid the
appearance of an unphysical phase factor in the cross section, the helicity
space utilised in the beam sector must synchronise with the one used for the
hard function in Eq. (2.30). In our case, since we use R  ECOLA to extract the
hard function, Eq. (2.37) is subject to the choice adopted in this program.
The quark beam function are known to N3LO accuracy [137, 138], while for the
gluon case only the helicity-conserving components
$\mathcal{B}_{++}^{[g_{n({\bar{n}})}]}$ and
$\mathcal{B}_{--}^{[g_{n({\bar{n}})}]}$ are known to this order [137, 138].
The helicity-flip components $\mathcal{B}_{+-}^{[g_{n({\bar{n}})}]}$ and
$\mathcal{B}_{-+}^{[g_{n({\bar{n}})}]}$ are only known to N2LO [136, 139,
140].

Finally, $\mathcal{S}^{\alpha\beta}_{[q_{n}\bar{q}_{{\bar{n}}}]}$ is the soft
function and covers the wide angle domain in Eq. (2.14). Its LP expression is,

$\displaystyle\mathcal{S}^{\alpha\beta}_{[q_{n}\bar{q}_{{\bar{n}}}]}(\vec{b}_{\mathrm{T}},\mu,\nu)$
$\displaystyle=$
$\displaystyle\;\sum_{\alpha^{\prime},\beta^{\prime}}\,\Big{(}\mathcal{Z}^{s,{\alpha\alpha^{\prime}}}_{[q_{n}\bar{q}_{{\bar{n}}}]}\Big{)}^{*}\,\mathcal{Z}^{s,\beta\beta^{\prime}}_{[q_{n}\bar{q}_{{\bar{n}}}]}\,\sum_{\\{a,a^{\prime},b\\}}\,\lim_{\delta\to
0^{+}}\,\left[c_{a^{\prime}_{1}a^{\prime}_{2}a^{\prime}_{3}a^{\prime}_{4}}^{qq,(\alpha^{\prime})}\right]^{*}\,c_{a_{1}a_{2}a_{3}a_{4}}^{qq,(\beta^{\prime})}\,$
(2.38) $\displaystyle\times\langle
0|\Big{[}\overline{Y}^{\dagger}_{n,b_{2}^{\phantom{\prime}}a^{\prime}_{2}}\overline{Y}_{{\bar{n}},a^{\prime}_{1}b_{1}^{\phantom{\prime}}}Y^{\dagger}_{v_{{\bar{t}}},b_{4}^{\phantom{\prime}}a^{\prime}_{4}}\,Y_{v_{t},a^{\prime}_{3}b_{3}^{\phantom{\prime}}}\Big{]}(-{ib_{0}\delta},-{ib_{0}\delta},\vec{b}_{\mathrm{T}})\,\Big{[}\overline{Y}^{\dagger}_{{\bar{n}},b_{1}a_{1}}\overline{Y}_{n,a_{2}b_{2}}Y^{\dagger}_{v_{t},b_{3}a_{3}}Y_{v_{{\bar{t}}},a_{4}b_{4}}\Big{]}(0)|0\rangle\,\bigg{|}_{\delta=\frac{1}{\nu}}\,,$
$\displaystyle\mathcal{S}^{\alpha\beta}_{[g_{n}g_{{\bar{n}}}]}(\vec{b}_{\mathrm{T}},\mu,\nu)$
$\displaystyle=$
$\displaystyle\;\sum_{\alpha^{\prime},\beta^{\prime}}\,\Big{(}\mathcal{Z}^{s,{\alpha\alpha^{\prime}}}_{[gg]}\Big{)}^{*}\,\mathcal{Z}^{s,{\beta\beta^{\prime}}}_{[gg]}\,\sum_{\\{a,a^{\prime},b\\}}\,\lim_{\delta\to
0^{+}}\,\,\left[c_{a^{\prime}_{1}a^{\prime}_{2}a^{\prime}_{3}a^{\prime}_{4}}^{gg,(\alpha^{\prime})}\right]^{*}\,c_{a_{1}a_{2}a_{3}a_{4}}^{gg,(\beta^{\prime})}$
(2.39) $\displaystyle\times\langle
0|\Big{[}\overline{\mathcal{Y}}^{\dagger}_{n,b_{2}^{\phantom{\prime}}a^{\prime}_{2}}\overline{\mathcal{Y}}_{{\bar{n}},a^{\prime}_{1}b_{1}^{\phantom{\prime}}}Y^{\dagger}_{v_{{\bar{t}}},b_{4}^{\phantom{\prime}}a^{\prime}_{4}}\,Y_{v_{t},a^{\prime}_{3}b_{3}^{\phantom{\prime}}}\Big{]}(-{ib_{0}\delta},-{ib_{0}\delta},\vec{b}_{\mathrm{T}})\,\Big{[}\overline{\mathcal{Y}}^{\dagger}_{{\bar{n}},b_{1}a_{1}}\overline{\mathcal{Y}}_{n,a_{2}b_{2}}Y^{\dagger}_{v_{t},b_{3}a_{3}}Y_{v_{{\bar{t}}},a_{4}b_{4}}\Big{]}(0)|0\rangle\,\bigg{|}_{\delta=\frac{1}{\nu}}\;,$

where $\mathcal{Z}^{s}_{[q_{n}\bar{q}_{{\bar{n}}}]}$ and
$\mathcal{Z}^{s}_{[gg]}$ are the renormalisation constants of the soft
function in the $\overline{\text{MS}}$ scheme. Again, $\delta$ denotes the
rapidity regulator [89, 90], and the $c_{\\{a\\}}^{qq(gg)}$ are the colour
coefficients defined in Eqs. (2.33)-(2.34). $\overline{Y}_{n({\bar{n}})}$ and
$\overline{\mathcal{Y}}_{n({\bar{n}})}$ describe the incoming soft Wilson
lines for the (anti-)quark and gluon fields, respectively, while the
$Y_{v_{t}(v_{{\bar{t}}})}$ are the outgoing soft Wilson lines of the
(anti-)top quark. Their specific expressions have been summarised in [127].
Even though the azimuthally averaged soft functions have been computed in the
recent years up to N2LO [55, 56, 141], the azimuthally resolved soft
functions, as displayed in Eqs. (2.38)-(2.39), are not yet available in the
context of effective field theories. The most relevant function at this point
is the soft correlation factor defined in [57]. However, the differing
strategies to renormalise the rapidity divergences between [57] and this
paper, make it challenging to thoroughly adapt their results to our work.
Hence, aiming at a self-consistent and independent study, in Sec. 2.5, we will
calculate the soft functions of Eqs. (2.38)-(2.39) at NLO accuracy,
disentangling the renormalisation scheme dependent logarithmic parts and the
constant contributions.

Now we can consider the differential cross section in Eq. (2.2) with the
reduced kinematic variables of Eqs. (2.3)-(2.20) as well as the expanded
partonic contributions in Eq. (2.26). We begin by disentangling the Fourier
integrals in Eq. (2.26) with the help of the reference vector $\vec{\tau}$,
using the different scaling behaviours of the components of the vector
$\vec{b}_{\mathrm{T}}$ (or correspondingly, $\vec{q}_{\mathrm{T}}$),

$\displaystyle\Sigma_{{t{\bar{t}}}}\Bigg{|}_{N_{\mathrm{J}}=0}\propto\,\sum_{\kappa}\int^{\infty}_{-\infty}{\rm
d}b_{\tau_{\|}}\exp\left(\mathrm{i}\,{q}_{\tau_{\|}}b_{\tau_{\|}}\right)\,\int^{\infty}_{-\infty}{\rm
d}b_{\tau_{\perp}}\exp\left(\mathrm{i}\,{q}_{\tau_{\perp}}b_{\tau_{\perp}}\right)\,\widetilde{\Sigma}_{{t{\bar{t}}}}^{[\kappa]}(\vec{b}_{\mathrm{T}},M_{t{\bar{t}}},\beta_{t{\bar{t}}},x_{t},Y_{t{\bar{t}}},\mu,\nu)+\dots\,,$
(2.40)

where in analogy to the case in Fig. 1, we apply the relationships,

$\displaystyle\vec{b}_{\mathrm{T}}=\vec{b}_{\tau_{\perp}}+\vec{b}_{\tau_{\|}}=b_{\tau_{\perp}}\vec{\tau}\times\vec{n}+b_{\tau_{\|}}\vec{\tau}\,,\quad\vec{q}_{\mathrm{T}}=\vec{q}_{\tau_{\perp}}+\vec{q}_{\tau_{\|}}=q_{\tau_{\perp}}\vec{\tau}\times\vec{n}+q_{\tau_{\|}}\vec{\tau}\,.$
(2.41)

Substituting the separated expression in Eq. (2.40) into Eq. (2.2) yields,

$\displaystyle\frac{{\rm d}^{5}{\sigma_{{t{\bar{t}}}}}}{{\rm
d}M_{t{\bar{t}}}^{2}\,{\rm d}^{2}\vec{P}_{t}^{\perp}\,{\rm
d}Y_{t{\bar{t}}}\,{\rm d}q_{\tau}}\Bigg{|}_{N_{\mathrm{J}}=0}\propto$
$\displaystyle\,\sum_{\kappa}\int^{\infty}_{-\infty}{\rm
d}b_{\tau_{\|}}\,\int^{\infty}_{-\infty}\,{\rm
d}q_{\tau_{\|}}\,\delta(q_{\tau}-|q_{\tau_{\|}}|)\,\exp\left(\mathrm{i}\,{q}_{\tau_{\|}}b_{\tau_{\|}}\right)\,$
(2.42) $\displaystyle\times\int^{\infty}_{-\infty}{\rm
d}b_{\tau_{\perp}}\,\widetilde{\Sigma}_{{t{\bar{t}}}}^{[\kappa]}(\vec{b}_{\mathrm{T}},M_{t{\bar{t}}},\beta_{t{\bar{t}}},x_{t},Y_{t{\bar{t}}},\mu,\nu)\,\int^{\infty}_{-\infty}\,{\rm
d}q_{\tau_{\perp}}\,\exp\left(\mathrm{i}\,{q}_{\tau_{\perp}}b_{\tau_{\perp}}\right)+\dots\,,$

where the variables that are independent of $\vec{b}_{\mathrm{T}}$ or
$\vec{q}_{\mathrm{T}}$ are omitted for simplicity. In comparison to Eq. (2.2),
one of the main differences in Eq. (2.42) resides in the absence of the
boundaries on the $q_{\tau_{\perp}}$ integral, which is a consequence of the
$\vec{q}_{\mathrm{T}}$ independence in the expanded function
$\Theta_{\mathrm{kin}}^{(0)}$ of Eq. (2.3). In this way, the integral over
$q_{\tau_{\perp}}$ in the second line of Eq. (2.42) can be completed before
the inverse Fourier transformation, thereby leading to the Dirac delta
function $\delta(b_{\tau_{\perp}})$. 222 It is worth emphasizing that this
operation is prohibited in the original expression of Eq. (2.2) due to the
explicit dependences on $\vec{q}_{\mathrm{T}}$ in the boundary function of Eq.
(2.3).  The integration over $b_{\tau_{\perp}}$ is then straightforward,
reducing the $\vec{b}_{\mathrm{T}}$ dependence of
$\widetilde{\Sigma}_{{t{\bar{t}}}}$ to a dependence on $\vec{b}_{\tau_{\|}}$
only. This leaves the integrals over $b_{\tau_{\|}}$ and $q_{\tau_{\|}}$. The
evaluation on the latter is immediate via the function
$\delta(q_{\tau}-|q_{\tau_{\|}}|)$ in the first line of Eq. (2.42),
eliminating the imaginary part of
$\exp(\mathrm{i}\,{q}_{\tau_{\|}}b_{\tau_{\|}})$. The integration over
$b_{\tau_{\|}}$, finally, is less straightforward to perform, and we turn to
the numeric solutions in Sec. 4.

To summarise, specifying the kinematic factors in Eq. (2.42), we are able to
establish the LP contribution from the $N_{\mathrm{J}}=0$ configuration,

$\displaystyle\frac{{\rm d}^{5}{\sigma_{{t{\bar{t}}}}}}{{\rm
d}M_{t{\bar{t}}}^{2}\,{\rm d}^{2}\vec{P}_{t}^{\perp}\,{\rm
d}Y_{t{\bar{t}}}\,{\rm
d}q_{\tau}}\Bigg{|}_{N_{\mathrm{J}}=0}=\sum_{\mathrm{sign}[x_{t}]}\frac{\Theta_{\mathrm{kin}}^{(0)}}{16\pi^{3}\,\beta_{t{\bar{t}}}\,|x_{t}|\,M^{4}_{{t{\bar{t}}}}\,s}\,\sum_{\kappa}\int^{\infty}_{-\infty}$
$\displaystyle{\rm
d}{b}_{\tau_{\|}}\,\cos\left({b}_{\tau_{\|}}{q}_{\tau}\right)\,$ (2.43)
$\displaystyle\times\widetilde{\Sigma}_{{t{\bar{t}}}}^{[\kappa]}(b_{\tau_{\|}}\vec{\tau},M_{t{\bar{t}}},\beta_{t{\bar{t}}},x_{t},Y_{t{\bar{t}}},\mu,\nu)\,,$

where the index $\kappa$ runs over
$\\{g_{n}g_{{\bar{n}}},q^{i}_{n}\bar{q}^{j}_{{\bar{n}}},q^{i}_{{\bar{n}}}\bar{q}^{j}_{n}\\}$
as before. $\Theta_{\mathrm{kin}}^{(0)}$ imposes the kinematic constraints as
defined in Eq. (2.20), and the $\widetilde{\Sigma}_{{t{\bar{t}}}}$ are the
contributions of the individual partonic processes, as presented in Eqs.
(2.29)-(2.30).

From Eq. (2.43), we are ready to assess the power accuracy of the
$N_{\mathrm{J}}=0$ configuration, which in fact can be fulfilled by counting
the power of each of its ingredients. First of all, since the kinematic
constraint $M_{t{\bar{t}}}>400$ GeV (or equivalently,
$\beta_{t{\bar{t}}}\gtrsim 0.5$) has been imposed, the prefactor in front of
the integral does not induce any power-like behaviour for the bulk of phase
space, and is thus of $\mathcal{O}{(1)}$. Next, in order to determine the
powers of the impact parameter ${b}_{\tau_{\|}}$ and the Fourier basis
$\cos\left({b}_{\tau_{\|}}{q}_{\tau}\right)$, note that the ${b}_{\tau_{\|}}$
integral serves in part as the momentum conservation condition on the beam and
soft radiations. Considering that the transverse momenta for the real
emissions are all of $\mathcal{O}(\lambda_{\tau})$ in the $N_{\mathrm{J}}=0$
configuration, we arrive at [78],

$\displaystyle{b}_{\tau_{\|}}\sim\mathcal{O}(\lambda^{-1}_{\tau})\,,\quad\quad\cos\left({b}_{\tau_{\|}}{q}_{\tau}\right)\sim\mathcal{O}{(1)}\,.$
(2.44)

The remaining task is to ascertain the power of
$\widetilde{\Sigma}_{{t{\bar{t}}}}$. As defined in Eqs. (2.29)-(2.30),
$\widetilde{\Sigma}_{{t{\bar{t}}}}$ entails the hard, soft, and beam functions
from the various partonic transitions. The hard sector consists of nothing but
the Wilson coefficients multiplied by the colour (helicity) bases, which
invokes no power-like behaviour and is, hence, of $\mathcal{O}{(1)}$. The soft
functions are defined in Eqs. (2.38)-(2.39) as the products of the soft Wilson
lines sandwiched by the vacuum states. From the scaling prescriptions in
SCETII [80, 81, 82], the soft Wilson lines are also of $\mathcal{O}{(1)}$,
regardless of the quark or gluon channels, which in turn yields,

$\displaystyle\mathcal{S}^{\alpha\beta}_{[q_{n}\bar{q}_{{\bar{n}}}]}\sim\mathcal{S}^{\alpha\beta}_{[g_{n}{g}_{{\bar{n}}}]}\sim\mathcal{O}{(1)}\,.$
(2.45)

The last piece to examine is the beam function, see Eqs. (2.35)-(2.36). It
contains the integral of the collinear building blocks,
$(W^{\dagger}_{n}\xi_{n})$ and $\mathcal{A}_{n}^{\perp}$ sandwiched between
the proton states, with respect to $b_{n}$. Akin to the case of
${b}_{\tau_{\|}}$, the power of $b_{n}$ is now related to $P_{{\bar{n}}}$,
giving $b_{n}\sim P_{{\bar{n}}}\sim\mathcal{O}{(1)}$. For the integrand, the
scaling rules of those operators and the corresponding external states are
$(W^{\dagger}_{n}\xi_{n})\sim\mathcal{A}_{n}^{\perp}\sim\mathcal{O}(\lambda_{\tau})$
and $|q_{n}\rangle\sim|g_{n}(\pm)\rangle\sim\mathcal{O}(\lambda^{-1}_{\tau})$
[80, 81, 82], respectively. We can thus conclude,

$\displaystyle\mathcal{B}_{n}^{[q_{i}]}\sim\mathcal{B}_{n,h^{\prime}h}^{[g_{n}]}\sim\mathcal{O}{(1)}\,.$
(2.46)

Combining the above scaling relationships, we observe that the only ingredient
of Eq. (2.43) that can bring about a power-like behaviour is the differential
${\rm d}{b}_{\tau_{\|}}$, from which the power of the $N_{\mathrm{J}}=0$
configuration is lowered by $\lambda_{\tau}$. In accordance, the $q_{\tau}$
spectrum behaves as,

$\displaystyle\frac{{\rm d}^{5}{\sigma_{{t{\bar{t}}}}}}{{\rm
d}M_{t{\bar{t}}}^{2}\,{\rm d}^{2}\vec{P}_{t}^{\perp}\,{\rm
d}Y_{t{\bar{t}}}\,{\rm
d}q_{\tau}}\Bigg{|}_{N_{\mathrm{J}}=0}\sim\mathcal{O}(\lambda^{-1}_{\tau})\,.$
(2.47)

Confronting Eq. (2.47) with the series in Eq. (2.6), it illustrates that Eq.
(2.43) can deliver at least in part the most singular behaviors of Eq. (2.6).
In order to assess the existence of other contributions to the leading
asymptotic terms in Eq. (2.6), we devote the following section to investigate
all possible $N_{\mathrm{J}}\geq 1$ configurations.

###  2.4 The case of $N_{\mathrm{J}}\geq 1$

Figure 3:  Dynamic regions in the $N_{\mathrm{J}}\geq 1$ configuration.
Herein, the black bubbles stand for the hard modes. The red bubble and curly
lines encode the soft corrections connecting the heavy partons represented by
the black double lines as well as both the beam- and jet-collinear modes
depicted in the blue straight lines.

This part will discuss the contributions induced by the hard, beam-collinear,
jet-collinear, and soft modes in the layouts with at least one hard jet. A
representative diagram of the associated dynamic regions is displayed in Fig.
3.

####  2.4.1 The $N_{\mathrm{J}}=1$ configuration

We start our analysis with a single insertion of a jet region. So far, as the
jet momentum is the sole source of energetic transverse recoil for the
${t{\bar{t}}}$ system, the components $q_{\tau_{\|}}$ and $q_{\tau_{\perp}}$
must admit the asymmetric configuration in Eq. (2.8) as a result of momentum
conservation. With this in mind, we can now expand the kinematic variables and
the boundary conditions of Eq. (2.2),

$\displaystyle M_{\mathrm{T}}^{t{\bar{t}}}\Bigg{|}_{N_{\mathrm{J}}=1}=$
$\displaystyle\;\sqrt{M_{t{\bar{t}}}^{2}+q^{2}_{\tau_{\perp}}}+\mathcal{O}(\lambda_{\tau})\,\equiv\,M_{\mathrm{T},\tau}^{t{\bar{t}}}\,+\mathcal{O}(\lambda_{\tau})\,,$
(2.48) $\displaystyle\Theta_{\mathrm{kin}}\Bigg{|}_{N_{\mathrm{J}}=1}=$
$\displaystyle\;\Theta\Big{[}\sqrt{s}-M_{\mathrm{T},\tau}^{t{\bar{t}}}-|q_{\tau_{\perp}}|\Big{]}\Theta\Big{[}M_{\mathrm{T},\tau}^{t{\bar{t}}}-m_{\mathrm{T},\tau}^{t}-m_{\mathrm{T},\tau}^{{\bar{t}}}\Big{]}$
$\displaystyle\hskip
56.9055pt\times\Theta\left[\sinh^{-1}\left(\sqrt{\frac{(M_{t{\bar{t}}}^{2}+s)^{2}}{4\,s\,{M_{\mathrm{T},\tau}^{t{\bar{t}}}}^{\\!\\!\\!2}}-1}\right)-|Y_{t{\bar{t}}}|\right]+\mathcal{O}(\lambda_{\tau})$
$\displaystyle\equiv$
$\displaystyle\;\Theta^{(1)}_{\mathrm{kin}}+\mathcal{O}(\lambda_{\tau})\,,$
(2.49) $\displaystyle P^{z}_{t}\Bigg{|}_{N_{\mathrm{J}}=1}=$
$\displaystyle\;\pm\,\frac{M_{\mathrm{T},\tau}^{t{\bar{t}}}}{2}\sqrt{1+\frac{m_{\mathrm{T},\tau}^{t}+m_{\mathrm{T},\tau}^{{\bar{t}}}}{M_{\mathrm{T},\tau}^{t{\bar{t}}}}}\,\sqrt{1+\frac{m_{\mathrm{T},\tau}^{t}-m_{\mathrm{T},\tau}^{{\bar{t}}}}{M_{\mathrm{T},\tau}^{t{\bar{t}}}}}\,\sqrt{1-\frac{m_{\mathrm{T},\tau}^{t}+m_{\mathrm{T},\tau}^{{\bar{t}}}}{M_{\mathrm{T},\tau}^{t{\bar{t}}}}}\,\sqrt{1-\frac{m_{\mathrm{T},\tau}^{t}-m_{\mathrm{T},\tau}^{{\bar{t}}}}{M_{\mathrm{T},\tau}^{t{\bar{t}}}}}\,$
$\displaystyle\;{}+\mathcal{O}(\lambda_{\tau})\,$ $\displaystyle\equiv$
$\displaystyle\;\widehat{P}^{z}_{t}\,+\mathcal{O}(\lambda_{\tau})\,,$ (2.50)

where the approximate transverse masses of the top and antitop quarks are
defined as
$m_{\mathrm{T},\tau}^{t}\equiv\sqrt{m^{2}_{t}+(\vec{q}_{\tau_{\perp}})^{2}}$
and
$m_{\mathrm{T},\tau}^{{\bar{t}}}\equiv\sqrt{m^{2}_{t}+(\vec{P}_{t}^{\perp}-\vec{q}_{\tau_{\perp}})^{2}}$,
respectively. The approximate ${t{\bar{t}}}$ invariant mass
$M_{\mathrm{T},\tau}^{t{\bar{t}}}$ is defined in Eq. (2.48). While Eqs.
(2.48)-(2.4.1) follow immediately from the definition of the transverse mass
${M}_{\mathrm{T}}$ and the boundary condition of Eq. (2.3), deriving Eq.
(2.4.1) necessitates solving the energy conservation equation in the zRF,

$\displaystyle{M}_{\mathrm{T}}=\sqrt{{m}_{\mathrm{T}}^{2}+(P^{z}_{t})^{2}}+\sqrt{({m}^{\prime}_{\mathrm{T}})^{2}+(P^{z}_{t})^{2}}\,.$
(2.51)

The solution is then expanded using the power counting of Eq. (2.8), keeping
the lowest power contributions. Using these results, we can now evaluate the
partonic contributions for $N_{\mathrm{J}}=1$. We follow the same steps as in
the derivation of Eq. (2.26) with the addition of embedding the soft Wilson
lines and the jet functions as appropriate here. The soft-collinear
decomposition as illustrated in Eqs. (2.23)-(2.25) is independent of the
specific configurations and, thus, still holds at present. After combining all
contributions and omitting the unrelated higher power correction terms, we
arrive at the partonic function $\Sigma_{{t{\bar{t}}}}$ of the
$N_{\mathrm{J}}=1$ case,

$\displaystyle\Sigma_{{t{\bar{t}}}}\Bigg{|}_{N_{\mathrm{J}}=1}=$
$\displaystyle\;\frac{1}{2\pi
s{\eta}^{\prime}_{n}{\eta}^{\prime}_{{\bar{n}}}}\sum_{\kappa,\lambda}\int{\rm
d}^{2}\vec{k}_{\mathrm{T}}\,{\rm d}y_{k}\,{\rm
d}^{2}\vec{b}_{\mathrm{T}}\,\exp\left(\mathrm{i}\,\vec{b}_{\mathrm{T}}\cdot\vec{q}_{\mathrm{T}}+\mathrm{i}\,\vec{b}_{\mathrm{T}}\cdot\vec{k}_{\mathrm{T}}\right)\widetilde{\Sigma}^{[\kappa]}_{{t{\bar{t}}},[\lambda]}(\vec{b}_{\mathrm{T}},\vec{P}_{t},\vec{k}_{\mathrm{T}},Y_{t{\bar{t}}},y_{k},\mu,\nu)\;,$
(2.52)

where now the indeces
$\kappa\in\\{g_{n}g_{{\bar{n}}},q^{i}_{n}\bar{q}^{j}_{{\bar{n}}},\bar{q}^{j}_{n}q^{i}_{{\bar{n}}},q^{i}_{n}g_{{\bar{n}}},g_{n}q^{i}_{{\bar{n}}},\bar{q}^{i}_{n}g_{{\bar{n}}},g_{n}\bar{q}^{i}_{{\bar{n}}}\\}$
and
$\lambda\in\\{g_{n_{\mathrm{J}}},{q}^{j}_{n_{\mathrm{J}}},\bar{q}^{j}_{n_{\mathrm{J}}}\\}$
label the initial and final state light partons, respectively.
$\vec{P}_{t}\equiv(\vec{P}_{t}^{\perp},\widehat{P}_{t}^{z})$ denotes the
spatial momentum of the top quark as measured in the zRF.
$\vec{k}_{\mathrm{T}}$ and $y_{k}$ stand for the transverse momentum and
rapidity of the jet in the LF, respectively. The variables
${\eta}^{\prime}_{n}$ and ${\eta}^{\prime}_{{\bar{n}}}$ are the momentum
fractions of the active partons in the beam functions with respect to the
colliding protons, namely,

$\displaystyle{\eta}^{\prime}_{n}$
$\displaystyle=(M_{\mathrm{T},\tau}^{t{\bar{t}}}\,e^{Y_{t{\bar{t}}}}+k_{\mathrm{T}}\,e^{y_{k}})/\sqrt{s}\,,$
(2.53) $\displaystyle{\eta}^{\prime}_{{\bar{n}}}$
$\displaystyle=(M_{\mathrm{T},\tau}^{t{\bar{t}}}\,e^{-Y_{t{\bar{t}}}}+k_{\mathrm{T}}\,e^{-y_{k}})/\sqrt{s}\,.$
(2.54)

In the following discussion, we will particularly focus on the
$q^{i}_{n}\bar{q}^{j}_{{\bar{n}}}\to{t{\bar{t}}}\,g_{n_{\mathrm{J}}}$ process.
All other partonic processes can be accessed through exchanging the labels
$n\leftrightarrow{\bar{n}}$ or the active partons herein, as appropriate. The
expression of
$\widetilde{\Sigma}^{[q^{i}_{n}\bar{q}^{j}_{{\bar{n}}}]}_{{t{\bar{t}}},[g_{n_{\mathrm{J}}}]}$
reads,

$\displaystyle\widetilde{\Sigma}^{[q^{i}_{n}\bar{q}^{j}_{{\bar{n}}}]}_{{t{\bar{t}}},[g_{n_{\mathrm{J}}}]}(\vec{b}_{\mathrm{T}},\vec{P}_{t},\vec{k}_{\mathrm{T}},Y_{t{\bar{t}}},y_{k},\mu,\nu)=$
$\displaystyle\;\left(\frac{1}{2N_{c}}\right)^{2}\,\mathcal{B}_{n}^{[q_{i}]}({\eta}^{\prime}_{n},b_{\mathrm{T}},\mu,\nu)\,\mathcal{B}_{{\bar{n}}}^{[\bar{q}_{j}]}({\eta}^{\prime}_{{\bar{n}}},b_{\mathrm{T}},\mu,\nu)\,\mathcal{J}^{[g]}_{{n_{\mathrm{J}}}}(\vec{k}_{\mathrm{T}},y_{k})\,$
$\displaystyle\;\;\times\sum_{\alpha,\beta}\left\\{\mathcal{H}_{\alpha\beta,[g_{n_{\mathrm{J}}}]}^{[q^{i}_{n}\bar{q}^{j}_{{\bar{n}}}]}(\vec{P}_{t},\vec{k}_{\mathrm{T}},Y_{t{\bar{t}}},y_{k},\mu)\,\mathcal{S}^{\alpha\beta,[g_{n_{\mathrm{J}}}]}_{[q_{n}\bar{q}_{{\bar{n}}}]}(\vec{b}_{\mathrm{T}},\mu,\nu)\,\right\\}\,,$
(2.55)

where the beam functions $\mathcal{B}_{n}^{[q_{i}]}$ and
$\mathcal{B}_{{\bar{n}}}^{[\bar{q}_{j}]}$ are defined in the same way as those
in Eqs. (2.35)-(2.36). The soft sector
$\mathcal{S}^{\alpha\beta,[g_{n_{\mathrm{J}}}]}_{[q_{n}\bar{q}_{{\bar{n}}}]}$
takes the similar appearance to that in eq. (2.38), except for the necessary
adaptation in the colour bases and the inclusion of the jet Wilson lines. The
jet function $\mathcal{J}^{[g]}_{{n_{\mathrm{J}}}}$ is the novel ingredient in
the $N_{\mathrm{J}}=1$ configuration,

$\displaystyle\mathcal{J}^{[g]}_{{n_{\mathrm{J}}}}(\vec{k}_{\mathrm{T}},y_{k})=-\frac{1}{4\pi(N_{c}^{2}-1)}\,\int{\rm
d}m^{2}_{k}\,{\rm d}^{4}x\,\exp\left(ik\cdot x\right)\,\mathbf{Tr}\,\langle
0|\mathcal{A}^{\bot}_{n_{\mathrm{J}}}(x)\mathcal{A}_{n_{\mathrm{J}}}^{\bot}(0)|0\rangle\,=\,1\,,$
(2.56)

where $m^{2}_{k}=k^{2}$ measures the offshellness of the jet radiations.
$\mathcal{A}_{n_{\mathrm{J}}}^{\bot}$ stands for the gauge-invariant collinear
fields in SCETII, see Eq. (2.36). In calculating Eq. (2.56), after using the
dimensional regulator to regulate the UV divergences, completing the
coordinate space integral in Eq. (2.56) always results in contributions of the
form $(m_{k}^{2}/\mu^{2})^{\epsilon}$ [142, 143, 144, 145, 146, 147]. The
following integration over the complete $m_{k}$ range, however, turns out to
be scaleless and, thus, the unmeasured jet functions involved in this paper
will never receive any perturbative corrections in $\alpha_{s}$. Therefore, we
equate $\mathcal{J}^{[g]}_{{n_{\mathrm{J}}}}$ in Eq. (2.56) with its tree-
level result. Further, it is worth pointing out that, since we are not
observing the jet itself, but only its recoil on the ${t{\bar{t}}}$ system,
the relative transverse momenta amongst the collinear emissions as well as
their helicity dependence have been integrated out.

The quark jet function shows the same behaviour [146, 147],

$\displaystyle\mathcal{J}^{[q]}_{{n_{\mathrm{J}}}}(\vec{k}_{\mathrm{T}},y_{k})=\,\frac{1}{8\pi
N_{c}{\bar{n}}_{\mathrm{J}}\cdot p}\,\int{\rm d}m^{2}_{k}\,{\rm
d}^{4}x\,\exp\left(ik\cdot x\right)\,\mathbf{Tr}\,\langle
0|\xi_{n_{\mathrm{J}}}(x)\bar{\xi}_{n_{\mathrm{J}}}(0)\not{{\bar{n}}}_{\mathrm{J}}|0\rangle\,=\,1\,.$
(2.57)

In addition to the above functions capturing the low-offshellness effects,
$\widetilde{\Sigma}^{[q^{i}_{n}\bar{q}^{j}_{{\bar{n}}}]}_{{t{\bar{t}}},[g_{n_{\mathrm{J}}}]}$
also needs the hard function
$\mathcal{H}_{\alpha\beta,[g_{n_{\mathrm{J}}}]}^{[q^{i}_{n}\bar{q}^{j}_{{\bar{n}}}]}$,
consisting of the UV-renormalised and IRC-subtracted partonic amplitude
$\mathcal{M}(q\bar{q}\to{t{\bar{t}}}\,g)$ constructed in analogy to Eqs.
(2.31)-(2.32). The presence of the jet mode, however, complicates the IRC-
subtraction procedure as singularities arising from the region
$|\vec{k}_{\mathrm{T}}|\to 0$ will need to be treated. This necessitates the
inclusion of lower-multiplicity partonic processes, e.g. by using

$\displaystyle\mathcal{Z}_{3}\otimes\mathcal{M}(q\bar{q}\to{t{\bar{t}}}\,g)+\mathcal{Z}_{2}\otimes\mathcal{M}(q\bar{q}\to{t{\bar{t}}})$
(2.58)

as the IRC-finite quantity [148, 149]. Its details, however, are unimportant
in the following as the only quantity of interest in this paper is the scaling
behaviour of the hard function
$\mathcal{H}_{\alpha\beta,[g_{n_{\mathrm{J}}}]}^{[q^{i}_{n}\bar{q}^{j}_{{\bar{n}}}]}$
itself.

Substituting Eq. (2.52) into Eq. (2.2), we obtain the master formula for the
$N_{\mathrm{J}}=1$ configuration,

$\displaystyle\frac{{\rm d}^{5}{\sigma_{{t{\bar{t}}}}}}{{\rm
d}M_{t{\bar{t}}}^{2}\,{\rm d}^{2}\vec{P}_{t}^{\perp}\,{\rm
d}Y_{t{\bar{t}}}\,{\rm d}q_{\tau}}\Bigg{|}_{N_{\mathrm{J}}=1}=$
$\displaystyle\;\sum_{\mathrm{sign}[\widehat{P}^{z}_{t}]}\sum_{\kappa,\lambda}\int{\rm
d}q_{\tau_{\perp}}\,{\rm
d}y_{k}\;\frac{\Theta^{(1)}_{\mathrm{kin}}}{8(2\pi)^{5}}\,\frac{\widetilde{\Sigma}_{{t{\bar{t}}},[\lambda]}^{[\kappa]}(\vec{0},\vec{P}_{t},-\vec{q}_{\tau_{\perp}},Y_{t{\bar{t}}},y_{k},\mu,\nu)}{M_{\mathrm{T},\tau}^{t{\bar{t}}}|\widehat{P}^{z}_{t}|{\eta}^{\prime}_{n}{\eta}^{\prime}_{{\bar{n}}}s^{2}}\,.$
(2.59)

In deriving this formula, the argument $\vec{b}_{\mathrm{T}}$ of the function
$\widetilde{\Sigma}_{{t{\bar{t}}},[\lambda]}^{[\kappa]}$ has been integrated
out following the multipole expansion. To see this, we apply the decomposition
in Eq. (2.41) again onto the impact parameter $\vec{b}_{\mathrm{T}}$ and the
jet transverse momentum $\vec{k}_{\mathrm{T}}$ in eq. (2.52),

$\displaystyle\vec{b}_{\mathrm{T}}=\vec{b}_{\tau_{\perp}}+\vec{b}_{\tau_{\|}}=b_{\tau_{\perp}}\vec{\tau}\times\vec{n}+b_{\tau_{\|}}\vec{\tau}\;,\qquad\vec{k}_{\mathrm{T}}=\vec{k}_{\tau_{\perp}}+\vec{k}_{\tau_{\|}}=k_{\tau_{\perp}}\vec{\tau}\times\vec{n}+k_{\tau_{\|}}\vec{\tau}\;.$
(2.60)

The parallel component ${k}_{\tau_{\|}}$ drops out of the lowest power hard
sector $\mathcal{H}_{\alpha\beta,[\lambda]}^{[\kappa]}$ during the asymptotic
expansion, due to the scaling hierarchy

$\displaystyle{k}_{\tau_{\|}}\sim\mathcal{O}(\lambda_{\tau})\ll{k}_{\tau_{\perp}}\sim\mathcal{O}(1)\,.$
(2.61)

Therefore, the ${k}_{\tau_{\|}}$ integral therein can be calculated
immediately, resulting in a $\delta(b_{\tau_{\|}})$ distribution in the
integral. The following ${b}_{\tau_{\|}}$ integral in
$\widetilde{\Sigma}_{{t{\bar{t}}},[\lambda]}^{[\kappa]}$ can then also be
carried out without further complications. On the other hand, in the
perpendicular direction, the expansion in $\lambda_{\tau}$ is subject to the
relationship

$\displaystyle q_{\tau_{\perp}}\sim k_{\tau_{\perp}}\sim\mathcal{O}(1)\gg
k^{\perp}_{s}\sim k^{\perp}_{c}\sim
k^{\perp}_{\bar{c}}\sim\mathcal{O}(\lambda_{\tau})\,,$ (2.62)

where $k^{\perp}_{s}$ and $k^{\perp}_{c(\bar{c})}$ represent the transverse
momenta of the soft and beam-collinear regions, as given in Eqs.
(2.14)-(2.16). Thus, the multipole expansion at the hard vertices can
eliminate the argument $\vec{b}_{\tau_{\perp}}$ of the beam-collinear and soft
functions up to $\mathcal{O}(\lambda_{\tau})$ [78], such that the inverse
Fourier transformation in Eq. (2.52) can be completed prior to the integral
over $k_{\tau_{\perp}}$, giving rise to
$\delta(k_{\tau_{\perp}}+q_{\tau_{\perp}})$ and, thus, sets
$\vec{k}_{\mathrm{T}}=-\vec{q}_{\tau_{\perp}}$ upon integration.

Using Eq. (2.59), we are now ready to determine the power accuracy of the
$N_{\mathrm{J}}=1$ configuration. Following the assessment of Eq. (2.43), for
the bulk of the phase space, the kinematic factors in Eq. (2.59), such as
$M_{\mathrm{T},\tau}^{t{\bar{t}}}$ and $\widehat{P}^{z}_{t}$, invoke no power-
like behaviour, and thus all belong to $\mathcal{O}(1)$. The hard, beam-
collinear, and soft functions are of $\mathcal{O}(1)$ by construction, and so
are the jet functions of Eqs. (2.56)-(2.57). Then, the remaining factors that
matter to the power counting are the differentials ${\rm d}q_{\tau_{\perp}}$
and ${\rm d}y_{k}$, which are of $\mathcal{O}(1)$ as well due to the scaling
laws in Eq. (2.8) and the definition of Eq. (2.17). Hence, all the ingredients
from Eq. (2.59) are characterised by the $\mathcal{O}(1)$ behaviour, which
then allows us to establish,

$\displaystyle\frac{{\rm d}{\sigma_{{t{\bar{t}}}}}}{{\rm
d}M_{t{\bar{t}}}^{2}\,{\rm d}^{2}\vec{P}_{t}^{\perp}\,{\rm
d}Y_{t{\bar{t}}}\,{\rm
d}q_{\tau}}\Bigg{|}_{N_{\mathrm{J}}=1}\sim\;\mathcal{O}(1)\,.$ (2.63)

Comparing Eq. (2.63) with the asymptotic series in Eq. (2.6) and the
${N_{\mathrm{J}}=0}$ scaling rule in Eq. (2.47), it is noted that the
$N_{\mathrm{J}}=1$ configuration here gives rise to the regular behaviors of
Eq. (2.6), which pertains to the NLP corrections and is one power higher than
the ${N_{\mathrm{J}}=0}$ influences.

####  2.4.2 The $N_{\mathrm{J}}\geq 2$ configuration

The case with (at least) two hard jet insertions differs from the pervious
case where the scaling laws for $q_{\tau_{\|}}$ and $q_{\tau_{\perp}}$ are
uniquely determined. The variety of the jet transverse momenta in the
$N_{\mathrm{J}}\geq 2$ configuration can accommodate both the isotropic and
asymmetric recoil configurations in Eqs. (2.7)-(2.8). In the following
paragraphs, we will discuss them individually.

##### Isotropic recoil.

We start with the isotropic recoil configuration. In this case, the transverse
components $q_{\tau_{\|}}$ and $q_{\tau_{\perp}}$ respect the power
prescription in Eq. (2.7), from which the reduced kinematics variables have
been presented in Eqs. (2.3)-(2.20). As demonstrated in Eq. (2.47), none of
them bears any power-like behaviour. Then, the problem reduces to the scaling
behaviour of the partonic function $\Sigma_{{t{\bar{t}}}}$ in Eq. (2.4) in the
$N_{\mathrm{J}}\geq 2$ configuration. The calculation of
$\Sigma_{{t{\bar{t}}}}$ now follows similarly to Eq. (2.52) in the
$N_{\mathrm{J}}=1$ case, aside from duplicating the jet function and
generalizing the hard and soft sectors in accordance, giving

$\displaystyle\Sigma_{{t{\bar{t}}}}\Bigg{|}_{N_{\mathrm{J}}\geq
2}\,\sim\,\int\,\bigg{[}\prod^{N_{\mathrm{J}}}_{i=1}{\rm d}y_{i}\,{\rm
d}^{2}\vec{k}^{\perp}_{i}\bigg{]}\,{\rm
d}^{2}\vec{b}_{\mathrm{T}}\,\exp\bigg{[}\mathrm{i}\;\vec{b}_{\mathrm{T}}\cdot\Big{(}\vec{q}_{\mathrm{T}}+\sum^{N_{\mathrm{J}}}_{j=1}\vec{k}^{\perp}_{j}\Big{)}\bigg{]}\,\mathcal{H}\,\otimes\mathcal{S}\,\otimes\mathcal{B}_{n}\,\otimes\mathcal{B}_{{\bar{n}}}\,\otimes\prod^{N_{\mathrm{J}}}_{l=1}J_{n_{l}}\,,$
(2.64)

where $y_{i}$ and $\vec{k}^{\perp}_{i}$ denote the pseudo-rapidity and
transverse momentum of the $i$-th jet. For simplicity, the indices specifying
the partonic channels have been omitted here. Akin to the $N_{\mathrm{J}}=1$
case, the hard, soft, beam, and jet functions in Eq. (2.64) are all of
$\mathcal{O}(1)$. To appraise the power accuracy of the pseudo-rapidity
$y_{i}$, it merits noting that due to the scaling rules in Eq. (2.12) and Eq.
(2.17), the differential ${\rm d}y_{i}$ always evaluates to $\mathcal{O}(1)$
and thus does not influence the scaling behaviour of $\Sigma_{{t{\bar{t}}}}$.

Now the remaining task is to determine the scaling behaviour of
$\vec{k}^{\perp}_{i}$ and $\vec{b}_{\mathrm{T}}$. It is noted that the scaling
behaviour of $\vec{b}_{\mathrm{T}}$ depends on the regions of
$\vec{k}^{\perp}_{i}$. To exhaust all the possibilities, we regroup the jets
here according to the scaling behaviour of the transverse components,

$\displaystyle\mathcal{P}_{k}\,\equiv$
$\displaystyle\;\\{\vec{k}_{i}^{\perp}\in\mathcal{R}_{k};\;\mathcal{O}(\lambda_{\tau})\sim
q_{\tau_{\perp}}\sim
q_{\tau_{\|}}\sim\vec{k}_{i,\tau_{\perp}}\ll\vec{k}_{i,\tau_{\|}}\sim\mathcal{O}(\lambda^{0}_{\tau})\\}\,,$
(2.65) $\displaystyle\mathcal{O}_{k}\,\equiv$
$\displaystyle\;\\{\vec{k}_{i}^{\perp}\in\mathcal{R}_{k};\;\mathcal{O}(\lambda_{\tau})\sim
q_{\tau_{\perp}}\sim
q_{\tau_{\|}}\sim\vec{k}_{i,\tau_{\|}}\ll\vec{k}_{i,\tau_{\perp}}\sim\mathcal{O}(\lambda^{0}_{\tau})\\}\,,$
(2.66) $\displaystyle\mathcal{I}_{k}\,\equiv$
$\displaystyle\;\\{\vec{k}_{i}^{\perp}\in\mathcal{R}_{k};\;\mathcal{O}(\lambda_{\tau})\sim
q_{\tau_{\perp}}\sim
q_{\tau_{\|}}\ll\vec{k}_{i,\tau_{\perp}}\sim\vec{k}_{i,\tau_{\|}}\sim\mathcal{O}(\lambda^{0}_{\tau})\\}\,,$
(2.67)

where the full set $\mathcal{R}_{k}$ collects all transverse momenta for the
$N_{\mathrm{J}}$ jets, more explicitly,

$\displaystyle\mathbf{card}[\mathcal{R}_{k}]=N_{\mathrm{J}}\,.$ (2.68)

Here the operator $\mathbf{card}$ evaluates the cardinality of a set.
$\mathcal{P}_{k}$, $\mathcal{O}_{k}$, and $\mathcal{I}_{k}$ are three non-
intersecting subsets of $\mathcal{R}_{k}$ consisting of different types of the
jet directions. For instance, $\mathcal{P}_{k}$ contains $\vec{k}_{i}^{\perp}$
parallel or antiparallel to the reference vector $\vec{\tau}$, while
$\mathcal{O}_{k}$ encompasses the orthogonal ones. $\mathcal{I}_{k}$ comprises
all the other configurations.

We are now ready to investigate the scalings for $\vec{b}_{\mathrm{T}}$. As
required by the label momentum conservation [78], the scaling power of
$b_{\tau_{\|}}$ ($b_{\tau_{\perp}}$) is subject to the strongest momenta in
the $\vec{\tau}$ ($\vec{n}\times\vec{\tau}$) direction. To this end, if there
are label momenta dictating both sides, namely,
$\mathbf{card}[\mathcal{P}_{k}]\,\mathbf{card}[\mathcal{O}_{k}]+\mathbf{card}[\mathcal{I}_{k}]\geq
1$, we have ${\rm d}\vec{b}_{\mathrm{T}}\sim\mathcal{O}(1)$. Otherwise, either
the $\vec{\tau}$ or the $\vec{n}\times\vec{\tau}$ direction will be governed
by $\mathcal{O}(\lambda_{\tau})$ fluctuations, which gives rise to ${\rm
d}\vec{b}_{\mathrm{T}}\sim\mathcal{O}(\lambda^{-1}_{\tau})$. Summarising these
relationships, we are capable of establishing the power counting for
$\Sigma_{{t{\bar{t}}}}$,

$\displaystyle\Sigma_{{t{\bar{t}}}}\Bigg{|}_{N_{\mathrm{J}}\geq
2,\mathbf{Iso}}\,\sim\,\left\\{\begin{array}[]{lr}\mathcal{O}{\left(\lambda_{\tau}^{\mathbf{card}[\mathcal{P}_{k}]+\mathbf{card}[\mathcal{O}_{k}]}\right)},&\mathbf{card}[\mathcal{P}_{k}]\,\mathbf{card}[\mathcal{O}_{k}]+\mathbf{card}[\mathcal{I}_{k}]\geq
1\,,\\\
\mathcal{O}{(\lambda_{\tau}^{N_{\mathrm{J}}-1})},&\mathbf{card}[\mathcal{P}_{k}]\,\mathbf{card}[\mathcal{O}_{k}]+\mathbf{card}[\mathcal{I}_{k}]=0\,,\end{array}\right.$
(2.71)

where the extra powers of $\lambda_{\tau}$ come from the differentials ${\rm
d}k_{i,\tau_{\perp}}\sim\mathcal{O}(\lambda_{\tau})$ of $\mathcal{P}_{k}$ and
${\rm d}k_{i,\tau_{\|}}\sim\mathcal{O}(\lambda_{\tau})$ of $\mathcal{O}_{k}$.
In light of the non-negative nature of the cardinality, the lowest power Eq.
(2.71) can reach is $\mathcal{O}(1)$, where the sets $\mathcal{P}_{k}$ and
$\mathcal{O}_{k}$ are both empty and thus $\mathcal{I}_{k}=\mathcal{R}_{k}$.
It should be emphasised that this finding is not dependent on the number of
the embedded jet modes, or, more specifically, the result of
$\mathbf{card}[\mathcal{I}_{k}]$. This differs from the naïve expectation from
the scaling behaviour of the effective Hamiltonian
$\mathcal{H}^{\mathrm{eff}}_{{t{\bar{t}}}}$, where it appears that the power
accuracy of $\mathcal{H}^{\mathrm{eff}}_{{t{\bar{t}}}}$ grows along with the
increase in the number of jets. The reason is that every jet in our
calculation is unmeasured and thus participates in the $q_{\tau}$ spectrum
such that, when calculating the contributions in Eqs. (2.56)-(2.57), the
4-dimensional coordinate space integral ${\rm
d}^{4}x\sim\mathcal{O}(\lambda^{-4}_{\tau})$ balances the power suppression
from the collinear field operators and the differential ${\rm d}m_{k}^{2}$.
Hence, the insertion of the jet modes invokes no power-like behaviour unless
kinematical constraints are imposed on the jet directions, such as those in
$\mathcal{P}_{k}$ or $\mathcal{O}_{k}$.

Substituting the result of Eq. (2.71) into Eq. (2.2), we arrive at

$\displaystyle\mathrm{Min}\left[\frac{{\rm d}{\sigma_{{t{\bar{t}}}}}}{{\rm
d}M_{t{\bar{t}}}^{2}\,{\rm d}^{2}\vec{P}_{t}^{\perp}\,{\rm
d}Y_{t{\bar{t}}}\,{\rm d}q_{\tau}}\Bigg{|}_{N_{\mathrm{J}}\geq
2,\mathbf{Iso}}\right]\,\sim\,\mathcal{O}(\lambda_{\tau})\,,$ (2.72)

where the integral over $q_{\tau_{\perp}}$ has increased the power of Eq.
(2.71) by one order of $\lambda_{\tau}$. In previous investigations on the
$q_{\mathrm{T}}$ spectrum, the scaling of the isotropic pattern was also
addressed in [65] and the outcome in Eq. (2.72) is in agreement with their
findings. Comparing Eq. (2.72) with the $N_{\mathrm{J}}=0$ and
$N_{\mathrm{J}}=1$ configurations, it is found that the result here is one
order higher than $N_{\mathrm{J}}=1$ case in Eq. (2.63), and two orders higher
with respect to the $N_{\mathrm{J}}=0$ one from Eq. (2.47).

##### Asymmetric recoil.

Turning now to the asymmetric recoil configuration, the transverse components
$q_{\tau_{\|}}$ and $q_{\tau_{\perp}}$ observe the scaling rules in Eq. (2.8).
The accordingly expanded boundary conditions can be found in Eqs.
(2.4.1)-(2.4.1). As demonstrated in Eq. (2.63), those kinematic factors are
all of $\mathcal{O}(1)$. Thus, the investigation of the scaling behaviour
driven by the asymmetric recoil configuration relies again on the analysis of
the convolution function in Eq. (2.64). As before, to denominate all the
possible configurations of the jet transverse momenta, we introduce the sets,

$\displaystyle\mathcal{P}_{k}\,\equiv$
$\displaystyle\;\\{\vec{k}_{i}^{\perp}\in\mathcal{R}_{k};\;\mathcal{O}(\lambda_{\tau})\sim
q_{\tau_{\|}}\sim\vec{k}_{i,\tau_{\perp}}\ll
q_{\tau_{\perp}}\sim\vec{k}_{i,\tau_{\|}}\sim\mathcal{O}(\lambda^{0}_{\tau})\\}\,,$
(2.73) $\displaystyle\mathcal{O}_{k}\,\equiv$
$\displaystyle\;\\{\vec{k}_{i}^{\perp}\in\mathcal{R}_{k};\;\mathcal{O}(\lambda_{\tau})\sim
q_{\tau_{\|}}\sim\vec{k}_{i,\tau_{\|}}\ll
q_{\tau_{\perp}}\sim\vec{k}_{i,\tau_{\perp}}\sim\mathcal{O}(\lambda^{0}_{\tau})\\}\,,$
(2.74) $\displaystyle\mathcal{I}_{k}\,\equiv$
$\displaystyle\;\\{\vec{k}_{i}^{\perp}\in\mathcal{R}_{k};\;\mathcal{O}(\lambda_{\tau})\sim
q_{\tau_{\|}}\ll
q_{\tau_{\perp}}\sim\vec{k}_{i,\tau_{\perp}}\sim\vec{k}_{i,\tau_{\|}}\sim\mathcal{O}(\lambda^{0}_{\tau})\\}\,.$
(2.75)

Since there are at least one pair of label momenta presiding over the
direction $\vec{n}\times\vec{\tau}$ now, the differential ${\rm
d}b_{\tau_{\perp}}$ is always of $\mathcal{O}(\lambda^{0}_{\tau})$. As for the
$\vec{\tau}$ direction, if there also appear any label momenta, namely
$\mathbf{card}[\mathcal{O}_{k}]+\mathbf{card}[\mathcal{I}_{k}]\geq 1$, the
power scaling of $b_{\tau_{\perp}}$ will give the same result as the
perpendicular piece, i.e. ${\rm
d}b_{\tau_{\perp}}\sim\mathcal{O}(\lambda^{0}_{\tau})$. Otherwise, this
direction will still be occupied by the soft and beam fluctuations, which
leads to ${\rm d}b_{\tau_{\perp}}\sim\mathcal{O}(\lambda^{-1}_{\tau})$.
Summarizing those observations, it follows,

$\displaystyle\Sigma_{{t{\bar{t}}}}\Bigg{|}_{N_{\mathrm{J}}\geq
2,\mathbf{Asy}}\,\sim\,\left\\{\begin{array}[]{lr}\mathcal{O}{\left(\lambda_{\tau}^{\mathbf{card}[\mathcal{P}_{k}]+\mathbf{card}[\mathcal{O}_{k}]}\right)},&\mathbf{card}[\mathcal{O}_{k}]+\mathbf{card}[\mathcal{I}_{k}]\geq
1\,,\\\
\mathcal{O}{(\lambda_{\tau}^{N_{\mathrm{J}}-1})},&\mathbf{card}[\mathcal{O}_{k}]+\mathbf{card}[\mathcal{I}_{k}]=0\,.\end{array}\right.$
(2.78)

Combining this result with Eq. (2.2) and exploiting the non-negative nature of
the cardinality, the minimal power of the cross section in the asymmetric
recoil configuration can be derived,

$\displaystyle\mathrm{Min}\left[\frac{{\rm d}{\sigma_{{t{\bar{t}}}}}}{{\rm
d}M_{t{\bar{t}}}^{2}\,{\rm d}^{2}\vec{P}_{t}^{\perp}\,{\rm
d}Y_{t{\bar{t}}}\,{\rm d}q_{\tau}}\Bigg{|}_{N_{\mathrm{J}}\geq
2,\mathbf{Asy}}\right]\,\sim\,\mathcal{O}(1)\,,$ (2.79)

where the scaling rule $q_{\tau_{\perp}}\sim\mathcal{O}(\lambda^{0}_{\tau})$
has been utilised in line with Eqs. (2.8). Comparing the outcome with Eq.
(2.63) of the $N_{\mathrm{J}}=1$ case, it is interesting to note that the
lowest power behaviour for the asymmetric configuration is not impacted by the
increase in the number of jets $N_{\mathrm{J}}$, both are of $\mathcal{O}(1)$.
However, in light of the asymptotic series in Eq. (2.6) and the power accuracy
of the $N_{\mathrm{J}}=0$ configuration in Eq. (2.47), it is noted that Eq.
(2.79) is only able to account for the regular terms in part, which belongs to
at least NLP and thus will not be the main concern in the latter numeric
evaluations.

####  2.4.3 Summary and discussion

Now that we have examined all possible configurations that are relevant to the
regime $q_{\tau}\to 0$ and determined the power precision of each contribution
with the help of the effective field theories SCETII and HQET. We find that
the leading-power contribution of $\mathcal{O}(\lambda^{-1}_{\tau})$ is given
by the $N_{\mathrm{J}}=0$ configuration in Eq. (2.47). It is followed at next-
to-leading power at $\mathcal{O}(\lambda^{0}_{\tau})$ by the asymmetric recoil
configurations regardless of the total number of the embedded jet modes, see
Eq. (2.63) for $N_{\mathrm{J}}=1$ and Eq. (2.79) for $N_{\mathrm{J}}\geq 2$,
respectively. The highest order contribution at $\mathcal{O}(\lambda_{\tau})$
is produced by the isotropic recoil configuration in Eq. (2.72) with the
$N_{\mathrm{J}}\geq 2$ jets having been incorporated.

Comparing those power rules with the asymptotic series in Eq. (2.6), we
observe that the leading singular terms are solely governed by the
$N_{\mathrm{J}}=0$ configuration. As will be illustrated below, this finding
imposes non-trivial constraints on the ingredients of Eqs. (2.29)-(2.30).
Firstly, due to the IRC-safe nature of the observable ${\rm
d}\sigma_{{t{\bar{t}}}}/{\rm d}q_{\tau}$, the asymptotic terms in Eq. (2.6)
are finite at each power in $\lambda_{\tau}$. In terms of the EFT ingredients,
this means the IRC subtraction factor of the hard sector in Eqs. (2.31)-(2.32)
and the renormalisation constants of the beam and soft functions in Eqs.
(2.35)-(2.36) and Eqs. (2.38)-(2.39) must cancel after the combination,

$\displaystyle\sum_{\alpha,\beta}\,\Big{(}\mathcal{Z}^{h,{\alpha^{\prime}\alpha}}_{[\kappa]}\mathcal{Z}^{s,{\alpha\alpha^{\prime\prime}}}_{[\kappa]}\Big{)}^{*}\,\mathcal{Z}^{h,{\beta^{\prime}\beta}}_{[\kappa]}\,\mathcal{Z}^{s,\beta\beta^{\prime\prime}}_{[\kappa]}\,\mathcal{Z}^{c}_{[\kappa_{n}]}\,\mathcal{Z}^{c}_{[\kappa_{{\bar{n}}}]}\,=\delta_{\alpha^{\prime}\alpha^{\prime\prime}}\delta_{\beta^{\prime}\beta^{\prime\prime}}\,,$
(2.80)

where the superscript $\kappa$ here runs again over the partonic channels
$\\{g_{n}g_{{\bar{n}}},q^{i}_{n}\bar{q}^{j}_{{\bar{n}}},q^{i}_{{\bar{n}}}\bar{q}^{j}_{n}\\}$
as in Eq. (2.26). $\kappa_{n}$ and $\kappa_{{\bar{n}}}$ mark the incoming
active partons along the $n$ and ${\bar{n}}$ directions, respectively. We have
omitted the contributions from $\mathcal{Z}^{cs}_{[\kappa_{n({\bar{n}})}]}$ of
Eq. (2.35) as it was solely introduced to remove the soft-collinear overlap
[135, 136] and thus is not involved with the renormalisation of the beam
function. The relationship in Eq. (2.80) will be examined in Sec. 2.5 through
an explicit NLO calculation.

From Eq. (2.80), we can infer the scale evolution relationship amongst the
elements of Eqs. (2.29)-(2.30). It is worth reminding that the scale
dependences in the hard, beam, and soft functions are all brought about
through the IRC subtraction or the UV/rapidity renormalisation. To this end,
the cancellation of the those subtraction factors and the renormalisation
constants renders the partonic convolution functions
$\widetilde{\Sigma}_{{t{\bar{t}}}}^{[\kappa]}$ in Eqs. (2.31)-(2.32)
independent of the scale $\mu$ or $\nu$, more specifically,

$\displaystyle\frac{{\rm d}}{{\rm
d}\ln\mu}\widetilde{\Sigma}_{{t{\bar{t}}}}^{[\kappa]}(\vec{b}_{\mathrm{T}},M_{t{\bar{t}}},\beta_{t{\bar{t}}},x_{t},Y_{t{\bar{t}}},\mu,\nu)=\frac{{\rm
d}}{{\rm
d}\ln\nu}\widetilde{\Sigma}_{{t{\bar{t}}}}^{[\kappa]}(\vec{b}_{\mathrm{T}},M_{t{\bar{t}}},\beta_{t{\bar{t}}},x_{t},Y_{t{\bar{t}}},\mu,\nu)=0\,.$
(2.81)

This result correlates the RGEs and RaGEs of the relevant ingredients therein,
and also permits us to leave out the scales $\mu$ and $\nu$ from the arguments
of $\widetilde{\Sigma}_{{t{\bar{t}}}}^{[\kappa]}$ hereafter. In Sec. 3, we
will utilise Eq. (2.81) to derive the evolution equations for the soft
interactions.

Apart from the leading singularities in Eq. (2.6), our analyses in Sec. 2.4.1
and Sec. 2.4.2 demonstrate that the following subleading terms entail the
participations of the jet modes. In the previous investigations, a similar
point was first addressed in [65]. There, the power scaling of the jet
contributions is determined by assessing the relevant pinch singularity
surfaces in the low $q_{\mathrm{T}}$ domain, from which the asymptotic
behaviors of the $q_{\mathrm{T}}$ spectra are related to the power laws of the
hard vertices333The focus of [65] is only on the Drell-Yan process. However,
in absence of the super renormalisable vertices, such as the Coulomb
exchanges, the findings of [65] are, however, generalisable to top quark pair
production at LHC. ,

$\displaystyle\frac{{\rm d}\sigma_{{t{\bar{t}}}}}{{\rm
d}^{2}\vec{q}_{\mathrm{T}}}\propto\delta^{2}(\vec{q}_{\mathrm{T}}-\vec{K}_{1}^{\perp}-\vec{K}_{2}^{\perp}),$
(2.82)

where $\vec{K}_{1,2}^{\perp}$ are the transverse momenta of two active partons
connecting the hard vertex producing the top quark pair. If
$\vec{K}_{1,2}^{\perp}$ and $\vec{q}_{\mathrm{T}}$, including all their
components, are of $\mathcal{O}(\lambda_{\tau})$, it yields ${\rm
d}\sigma_{{t{\bar{t}}}}/{\rm
d}^{2}\vec{q}_{\mathrm{T}}\sim\mathcal{O}(\lambda^{-2}_{\tau})$, which is in
agreement with our result in Eq. (2.47) from the $N_{\mathrm{J}}=0$
configuration after integrating out $q_{\tau_{\perp}}$. When the jet modes are
taken into account in the isotropic recoil configuration in Eq. (2.7), it
follows that $\vec{q}_{\mathrm{T}}\sim\mathcal{O}(\lambda_{\tau})$ and
$\vec{K}_{1,2}^{\perp}\sim\mathcal{O}(1)$, which gives rise to the scaling
behaviour ${\rm d}\sigma_{{t{\bar{t}}}}/{\rm
d}^{2}\vec{q}_{\mathrm{T}}\sim\mathcal{O}(1)$ from Eq. (2.82), also coinciding
with the expression in Eq. (2.72) once the integral over $q_{\tau_{\perp}}$ is
performed. This congruence of findings is a consequence of the equivalence of
Eq. (2.82) and our approach in regards to the isotropic configuration in Eq.
(2.7). During our derivations, the power behaviour of the $q_{\tau}$ spectrum
is extracted in part from the scaling laws of the impact parameter
$\vec{b}_{\mathrm{T}}$, which is in practice correlated to the delta function
in Eq. (2.82) by means of the inverse Fourier transformation, as illustrated
in (2.44).

Nevertheless, once the asymmetric recoil configuration of Eq. (2.8) is
encountered, it is not straightforward to apply Eq. (2.82). The power rules
for the transverse components are
$q_{\tau_{\|}}\sim\mathcal{O}(\lambda_{\tau})$ and
$q_{\tau_{\perp}}\sim\mathcal{O}(1)$ here, from which the right-hand side of
Eq. (2.82) is of $\mathcal{O}(\lambda^{-1}_{\tau})$ in the $N_{\mathrm{J}}=1$
situation and $\mathcal{O}(1)$ in the $N_{\mathrm{J}}\geq 2$ case. Even though
this accidentally agrees with our EFT-based derivation in the
$N_{\mathrm{J}}\geq 2$ case, see Eq. (2.79), its prediction in the
$N_{\mathrm{J}}=1$ case is one power lower than our EFT-based outcome in Eq.
(2.63). This mismatch originates in part in the fact that in [65] all the jet
transverse momenta possess homogenous components, as is the case our
$\mathcal{I}_{k}$ classification in Eq. (2.67) and Eq. (2.75). This
arrangement works well in the isotropic configuration of Eq. (2.7). However,
if the jet orientations are of particular concern, such as in the
$N_{\mathrm{J}}=1$ case, the extra power suppression will come into play by
means of the integrals over the transverse components, e.g. the projected
element $k_{\tau_{\|}}$ detailed in Eq. (2.61).

###  2.5 The soft function with the exponential regulator

(a) Light-light correlation
H

(b) Light-heavy correlation
H

(c) heavy-heavy correlation
H

Figure 4:  Representative Feynman diagrams for the NLO soft function. The
double line represents the gauge links attached to the top or antitop quark.
The single line stands for that from the incoming massless parton, such as
gluon, quark, and antiquark.

In the previous parts, we have derived the factorisation formula for the
leading singular behaviour of the $q_{\tau}$ spectrum in Eqs. (2.29)-(2.30),
which entails the soft functions $\mathcal{S}^{\alpha\beta}_{[\dots]}$ to
accommodate the wide angle correlations amongst the active partons. The field-
operator definitions of the soft elements have been presented in Eqs.
(2.38)-(2.39) in terms of the soft Wilson lines, see [127], sandwiched by the
orthonormal colour basis in Eqs. (2.33)-(2.34). In this part, we will
calculate the soft function to NLO accuracy. Their three typical contributions
are depicted in Fig. 4.

Without any loss of generality, the fixed-order results can be parameterised
as

$\displaystyle\mathcal{S}^{\alpha\beta}_{[\kappa]}(\vec{b}_{\mathrm{T}},\mu,\nu)\equiv\sum_{m=0}^{\infty}\,\left[\frac{\alpha_{s}(\mu)}{4\pi}\right]^{m}\,\mathcal{S}^{\alpha\beta,(m)}_{[\kappa]}(\vec{b}_{\mathrm{T}},\mu,\nu)\,,$
(2.83)

where
$\kappa\in\\{g_{n}g_{{\bar{n}}},q^{i}_{n}\bar{q}^{j}_{{\bar{n}}},q^{i}_{{\bar{n}}}\bar{q}^{j}_{n}\\}$
as that in Eq. (2.26). The coefficients for the first two orders in
$\alpha_{s}$ read

$\displaystyle\mathcal{S}^{\alpha\beta,(0)}_{[\kappa]}(\vec{b}_{\mathrm{T}},\mu,\nu)\,=$
$\displaystyle\;\delta_{\alpha\beta}\,,$ (2.84)
$\displaystyle\mathcal{S}^{\alpha\beta,(1)}_{[\kappa]}(\vec{b}_{\mathrm{T}},\mu,\nu)\,=$
$\displaystyle\;\sum_{a,b}\,\langle
c^{\alpha}_{\kappa}|\,\mathbf{T}_{a}\\!\cdot\mathbf{T}_{b}\,|c^{\beta}_{\kappa}\rangle\;\mathcal{I}_{ab}(\vec{b}_{\mathrm{T}},\mu,\nu)\,.$
(2.85)

In absence of any soft interactions, the LO result
$\mathcal{S}^{\alpha\beta,(0)}_{[\kappa]}$ is equal to the identity matrix as
a consequence of the orthonormality of our colour basis in Eqs. (2.33)-(2.34).
For the NLO expression, we employ the colour algebra formalism suggested in
[150] to illustrate the outcome. Therein, the flavour subscripts
$a,b\in\\{q,\bar{q},g,t,\bar{t}\\}$ denote the active partons participating in
the hard kernel. $\mathbf{T}_{a}$ signifies the colour charge operator for the
parton $a$, and $\mathbf{T}_{b}$ likewise. $|c_{\kappa}\rangle$ is the vector
representation of $c^{\kappa}_{\\{a\\}}$ in Eqs. (2.33)-(2.34) in colour
space. The coefficient function $\mathcal{I}_{ab}$ contains the contributions
from the squared soft amplitudes induced by the Wilson lines in Eqs.
(2.38)-(2.39). Its expression in the exponential regularisation scheme [89,
90] is

$\displaystyle\mathcal{I}_{ab}(\vec{b}_{\mathrm{T}},\mu,\nu)\equiv\,-\frac{\left(\mu^{2}e^{\gamma_{\mathrm{E}}}\pi\right)^{\epsilon}}{\pi}\,\lim_{\delta\to
0}\int{\rm d}k_{n}\,{\rm d}k_{{\bar{n}}}\,{\rm
d}^{2-2\epsilon}\;\vec{k}_{\mathrm{T}}\,\delta^{+}(k^{2})\,\frac{v_{a}\cdot
v_{b}}{(k\cdot v_{a})(k\cdot
v_{b})}\,e^{\mathrm{i}\vec{k}_{\mathrm{T}}\cdot\vec{b}_{\mathrm{T}}-b_{0}E_{k}\delta}\Bigg{|}_{\delta\to\frac{1}{\nu}}\,+\mathcal{Z}_{\mathrm{c}.\mathrm{t}.}^{ab}\,,$
(2.86)

where $b_{0}=2\exp(-\gamma_{\mathrm{E}})$ with $\gamma_{\mathrm{E}}$ being the
Euler constant. $\epsilon$ and $\delta$ are the regulators for the virtuality
and rapidity divergences, respectively. $v_{a}$ denotes the four-velocity of
parton $a$, for instance $n^{\mu}({\bar{n}}^{\mu})$ for $a=q_{n({\bar{n}})}$
and $v^{\mu}_{t(\bar{t})}=P^{\mu}_{t(\bar{t})}/m_{t}$ for the
$a=t({\bar{t}})$. $\mathcal{Z}_{\mathrm{c}.\mathrm{t}.}$ stems from the
perturbative expansion of the renormalisation constant
$\mathcal{Z}_{[\kappa]}^{s}$ and its complex conjugate in Eqs. (2.38-2.39),
defined in the $\overline{\text{MS}}$ scheme throughout this paper.

The result of Eq. (2.86) depends on the partons $a$ and $b$. We will first
examine the case of $a$ and $b$ being light flavours,
$a,b=l\in\\{q,\bar{q},g\\}$, depicted in Fig. 4(a). Starting with the case
$a=b=l$, the light-like nature of the vectors $n$ and ${\bar{n}}$ trivialises
the calculation and we have,

$\displaystyle\mathcal{I}_{ll}(\vec{b}_{\mathrm{T}},\mu,\nu)=0\,,\qquad\mathcal{Z}_{\mathrm{c}.\mathrm{t}.}^{ll}(\vec{b}_{\mathrm{T}},\mu,\nu,\epsilon)=0\,.$
(2.87)

However, if the participants consist of different light particles, $a=l\neq
b=l^{\prime}$, $\mathcal{I}_{ll^{\prime}}$ does not vanish in general. To
calculate its contribution, we follow [89] and first integrate out the
longitudinal components $k_{n}$ and $k_{{\bar{n}}}$ in line with the
exponential regularisation. Then, we expand the result in $\delta$ and
truncate the series to $\mathcal{O}(\delta^{0})$. Finally, we carry out the
integral over $\vec{k}_{\mathrm{T}}$ with the aid of the dimensional
regulator, giving

$\displaystyle\mathcal{I}_{ll^{\prime}}(\vec{b}_{\mathrm{T}},\mu,\nu)=L^{2}_{\mathrm{T}}-2\,L_{\mathrm{T}}\,L_{\nu}+\frac{\pi^{2}}{6}\,,\qquad\mathcal{Z}_{\mathrm{c}.\mathrm{t}.}^{ll^{\prime}}(\vec{b}_{\mathrm{T}},\mu,\nu,\epsilon)=\frac{2}{\epsilon^{2}}+\frac{2\,}{\epsilon}L_{\nu}\,,$
(2.88)

where $L_{\nu}=\ln[\mu^{2}/\nu^{2}]$ and
$L_{\mathrm{T}}=\ln[b_{\mathrm{T}}^{2}\mu^{2}/b_{0}^{2}]$. We have compared
Eq. (2.88) to the expressions in [89] and find the agreement after
synchronising the overall colour factors.

In addition to those correlations between the incoming particles, Eq. (2.86)
also involves the contributions from the heavy top quark,
$h\in\\{t,\bar{t}\\}$, see Fig. 4(b). The presence of the massive partons
complicates the calculation substantially. Since the denominators at this
moment are not homogeneous in $k_{n}$ and $k_{{\bar{n}}}$, performing the
integral over those longitudinal components is not straightforward, especially
when involving the exponent $E_{k}\delta$. To this end, we resort to the
Mellin-Barnes (MB) representation [151, 152] to recast the inhomogeneous
propagators in a first step. For $\mathcal{I}_{lh}$, we thus apply the
following substitution [71, 153],

$\displaystyle\frac{1}{(k\cdot v_{h})^{\lambda}}\,\xRightarrow[\text{
}]{\textbf{M.B. }}\,-\frac{1}{\Gamma(\lambda)}\frac{1}{4\pi^{2}}$
$\displaystyle\int^{c_{i}+\mathrm{i}\infty}_{c_{i}-\mathrm{i}\infty}{\rm
d}z_{1}\,{\rm
d}z_{2}\;\Gamma(-z_{1})\,\Gamma(-z_{2})\,\Gamma(\lambda+z_{1}+z_{2})$ (2.89)
$\displaystyle\qquad\times\,\left[\frac{(k\cdot
n)(v_{h}\cdot{\bar{n}})}{2}\right]^{z_{1}}\,\left[\frac{(k\cdot{{\bar{n}}})(v_{h}\cdot{n})}{2}\right]^{z_{2}}\,\left(-\vec{k}_{\mathrm{T}}\cdot\vec{v}^{\perp}_{h}\right)^{-\lambda-
z_{1}-z_{2}}\,.$

where the contours (or the values $c_{i}$) are chosen such that the poles from
$\Gamma(\lambda+z_{1}+z_{2})$ are to the left of the path, while those of
$\Gamma(-z_{1})\Gamma(-z_{2})$ are to the right. Integrating the MB-
transformed propagators over $k_{n{({\bar{n}})}}$ and $\vec{k}_{\mathrm{T}}$
now follows the similar pattern to that in deriving Eq. (2.88). Then, to
perform the $\delta$-expansion, we first make use of the package MB [154] to
determine the contours for $z_{1}$ and $z_{2}$, and then feed the outputs to
MBasymptotics [155] for deriving the asymptotic series in $\delta$. The
remaining integrals are those over $z_{1}$ and $z_{2}$, for which we utilise
MBsums [156] to implement Cauchy’s residue theorem and sum up the ensuing
residues from the software Mathematica. The final expression is,

$\displaystyle\mathcal{I}_{lh}(\vec{b}_{\mathrm{T}},\mu,\nu)\,=$
$\displaystyle\;\frac{L^{2}_{\mathrm{T}}}{2}-L_{\mathrm{T}}\,L_{\nu}+2\,L_{\mathrm{T}}\,\ln\left(v_{l}\cdot
v_{h}\right)$
$\displaystyle\;{}+2\,\mathrm{arcsinh}^{2}\Big{[}|\vec{v}_{h}^{\perp}|\sin(\varphi_{hb})\Big{]}-(2\pi\mathrm{i})\,\mathrm{arcsinh}\Big{[}|\vec{v}_{h}^{\perp}|\cos(\varphi_{hb})\Big{]}\,$
$\displaystyle\;{}+\frac{\pi^{2}}{12}+\mathcal{D}_{lh}(\varphi_{hb})+\mathcal{D}_{lh}(-\varphi_{hb})\,,$
(2.90)
$\displaystyle\mathcal{Z}_{\mathrm{c}.\mathrm{t}.}^{lh}(\vec{b}_{\mathrm{T}},\mu,\nu,\epsilon)\,=$
$\displaystyle\;\frac{1}{\epsilon^{2}}+\frac{1}{\epsilon}\Big{[}L_{\nu}-2\ln(v_{l}\cdot
v_{h})\Big{]}\,,$ (2.91)

where $\varphi_{hb}$ denotes the azimuthal opening angle between the vectors
$\vec{v}_{h}^{\perp}$ and $\vec{b}_{\mathrm{T}}$. The function
$\mathcal{D}_{lh}$ is defined as

$\displaystyle\mathcal{D}_{lh}(\theta)=\frac{4|\vec{v}_{h}^{\perp}|}{\pi}\,\int_{0}^{\frac{\pi}{2}}{\rm
d}{\phi}\,\ln\left[\cot(\phi)\right]\,\Bigg{\\{}\frac{\cos(\theta+\phi)}{\sqrt{1+|\vec{v}_{h}^{\perp}|^{2}\sin^{2}(\theta+\phi)}}\arcsin\bigg{[}\frac{|\vec{v}_{h}^{\perp}|\cos(\theta+\phi)}{\sqrt{1+|\vec{v}_{h}^{\perp}|^{2}}}\bigg{]}\Bigg{\\}}\,.$
(2.92)

An analogous strategy is also applicable for the evaluations of the $a=b=h$
and $a=h\neq b=h^{\prime}$ cases, $\mathcal{I}_{hh}$ and
$\mathcal{I}_{hh^{\prime}}$, of Fig. 4(c). Since rapidity divergences do not
emerge from those time-like gauge links, we can set the regulator $\delta=0$
from the beginning and solve the $k^{\mu}$-integrals in the conventional
dimensional regularisation. The results exhibit explicit $z_{1}$ and $z_{2}$
dependences, which are treated with the above MB-Tools [154, 156] to complete
the inverse MB transformations. Eventually, they yield,

$\displaystyle\mathcal{I}_{hh}(\vec{b}_{\mathrm{T}},\mu,\nu)\,=$
$\displaystyle\;2\,L_{\mathrm{T}}\,+\frac{{r}^{h}_{b}}{\sqrt{1+({r}^{h}_{b})^{2}}}\Bigg{[}4\,\mathrm{arcsinh}({r}^{h}_{b})-(2\pi\mathrm{i})\Bigg{]}\,,$
(2.93)
$\displaystyle\mathcal{Z}^{hh}_{\mathrm{c.t.}}(\vec{b}_{\mathrm{T}},\mu,\nu,\epsilon)\,=$
$\displaystyle\;-\frac{2}{\epsilon}\,,$ (2.94)
$\displaystyle\mathcal{I}_{hh^{\prime}}(\vec{b}_{\mathrm{T}},\mu,\nu)\,=$
$\displaystyle\;L_{\mathrm{T}}\,\frac{4\,v_{h}\cdot
v_{h^{\prime}}}{\sqrt{(v_{h}\cdot
v_{h^{\prime}})^{2}-1}}\,\mathrm{arctanh}\bigg{[}\sqrt{\frac{v_{h}\cdot
v_{h^{\prime}}-1}{v_{h}\cdot v_{h^{\prime}}+1}}\bigg{]}$
$\displaystyle\;+\int_{0}^{1}{\rm
d}\zeta\,\frac{{r}^{hh^{\prime}}_{b}(\zeta)}{\sqrt{1+\left({r}^{hh^{\prime}}_{b}(\zeta)\right)^{2}}}\,\left(\frac{v_{h}\cdot
v_{h^{\prime}}}{{v}^{2}_{hh^{\prime}}(\zeta)}\right)\,\Bigg{[}4\,\mathrm{arcsinh}\left({r}^{hh^{\prime}}_{b}(\zeta)\right)-(2\pi\mathrm{i})\Bigg{]}\,,$
(2.95)
$\displaystyle\mathcal{Z}^{hh^{\prime}}_{\mathrm{c.t.}}(\vec{b}_{\mathrm{T}},\mu,\nu,\epsilon)\,=$
$\displaystyle\;-\frac{1}{\epsilon}\frac{4\,v_{h}\cdot
v_{h^{\prime}}}{\sqrt{(v_{h}\cdot
v_{h^{\prime}})^{2}-1}}\,\mathrm{arctanh}\bigg{[}\sqrt{\frac{v_{h}\cdot
v_{h^{\prime}}-1}{v_{h}\cdot v_{h^{\prime}}+1}}\bigg{]}\,,$ (2.96)

where

$\displaystyle{v}^{\mu}_{hh^{\prime}}(\zeta)\,=\;\zeta\,v_{h}^{\mu}+(1-\zeta)\,v_{h^{\prime}}^{\mu}\,,\qquad{r}^{h}_{b}\,=\;\widehat{b}_{\mathrm{T}}\cdot\vec{v}^{\perp}_{h}\,,\qquad{r}^{hh^{\prime}}_{b}(\zeta)\,=\;\widehat{b}_{\mathrm{T}}\cdot\frac{\vec{v}_{hh^{\prime}}^{\perp}(\zeta)}{\sqrt{v^{2}_{hh^{\prime}}(\zeta)}}\,.$
(2.97)

During the calculations, we have exploited the on-shell conditions
$v_{t}^{2}=v_{\bar{t}}^{2}=1$ and introduced the definition
$\widehat{b}_{\mathrm{T}}\equiv\vec{b}_{\mathrm{T}}/b_{\mathrm{T}}$.

With all relevant coefficient functions $\mathcal{I}_{ab}$ at hand, together
with the colour factors of Eq. (2.85), we can now establish the NLO soft
function together with the renormalisation constants
$\mathcal{Z}^{s}_{[\kappa]}$. For the latter case, we have checked that the
results in Eqs. (2.88)-(2.96) indeed satisfy the identity in Eq. (2.80), where
the hard and beam renormalisation constants are extracted from [128] and [135,
136], respectively. Regarding the renormalised finite parts
$\mathcal{I}_{ab}$, we compare our expressions with those from the CSS
framework [57] at the scale $\mu=\nu=b_{0}^{2}/b_{\mathrm{T}}^{2}$ and find
full agreement. Furthermore, since during the derivation we did not utilise
the relationship $\vec{v}_{t}^{\perp}+\vec{v}_{\bar{t}}^{\perp}=\vec{0}$ for
simplification, the results in Eqs. (2.88-2.5) are also comparable with the
soft function in the process $pp\to{t{\bar{t}}}H$ [157]. However, at this
moment, although the real parts of $\mathcal{I}_{ab}$ still coincide with
those in [157], the signs in front of the terms $2\pi\mathrm{i}$ in Eq. (2.5),
Eq. (2.93), and Eq. (2.5) are inverted. This purely imaginary term does not
enter the present calculation on the ${t{\bar{t}}}$ production or the
comparison with [57], as these imaginary contributions cancel through the
momentum conservation $\vec{v}_{t}^{\perp}+\vec{v}_{\bar{t}}^{\perp}=\vec{0}$.
Nonetheless, this difference can influence the transverse momentum spectrum in
${t{\bar{t}}}H$ production, or similar processes, especially within the domain
$M_{{t{\bar{t}}}H}\gtrsim M_{t{\bar{t}}}\sim 2m_{t}$. We leave it to a
forthcoming publication to elaborate on those processes and, in particular,
deliver a numeric comparison of our results, to those from [157], and a fixed-
order QCD calculation from S  HERPA [158, 159, 160].

##  3 Resummation

In the previous section we have analysed the dynamic regions contributing to
the leading singularities in Eq. (2.6) and established the corresponding
factorisation formula in impact parameter space. This section will now be
devoted to the resummation of these asymptotic behaviours.

###  3.1 Asymptotic behavior

We start with identifying the singular terms in Eq. (2.43). Therefore, without
loss of generality, we parametrise the perturbative expansion of
$\widetilde{\Sigma}_{{t{\bar{t}}}}$ as follows,

$\displaystyle\widetilde{\Sigma}_{{t{\bar{t}}}}\,\sim\,\sum_{m,n}\,\alpha^{m}_{s}(M_{t{\bar{t}}})\;L^{n}_{\mathrm{M}}\,\Big{\\{}s_{m,n}(\beta_{t{\bar{t}}},x_{t},Y_{t{\bar{t}}})+a_{m,n}(\mathrm{sign}[b_{\tau_{\|}}],\beta_{t{\bar{t}}},x_{t},Y_{t{\bar{t}}})\Big{\\}}\;,$
(3.1)

where
$L_{\mathrm{M}}\equiv\log\left[{{b}^{2}_{\tau_{\|}}M_{t{\bar{t}}}^{2}}/{b_{0}^{2}}\right]$
collects all the relevant dimensionful quantities. Order by order, there are
two dimensionless coefficients, $s_{m,n}$ and $a_{m,n}$. While both depend on
$\beta_{t{\bar{t}}}$, $x_{t}$, and $Y_{t{\bar{t}}}$, $s_{m,n}$ is independent
of the magnitude and orientation of the impact parameter. Thus, we will refer
them to the azimuthal _symmetric_ terms (AST) hereafter. On the other hand, in
presence of the helicity-flipping beam radiation in Eq. (2.36) and the wide
angle soft correlations from Eqs. (2.38)-(2.39),
$\widetilde{\Sigma}_{{t{\bar{t}}}}$ also contains contributions sensitive to
the orientation of the impact parameter, i.e. the azimuthal asymmetric term
(AAT) $a_{m,n}$ [64, 58]. Since the rejection component $b_{\tau_{\perp}}$ has
been integrated out in Eq. (2.42), the orientation dependence is reduced to a
dependence on the sign of $b_{\tau_{\|}}$ in Eq. (3.1). Please note, in
deriving eq. (3.1), we have set the scales $\mu=\nu=M_{t{\bar{t}}}$ in Eqs.
(2.29)-(2.30) and then extracted the $L_{\mathrm{M}}$ terms through
dimensional analysis. Other scales choices $\\{\mu,\nu\\}$ lead to the same
result after regrouping the dimensional quantities appropriately, on account
of the (rapidity) renormalisation group discussed in Sec. 2.4.3.

From Eq. (3.1), we can evaluate the $q_{\tau}$ spectrum by completing the
inverse Fourier transformation in Eq. (2.43). Exploiting the fact that the
function $\cos(b_{\tau_{\|}}q_{\tau})$ is even under the integral over
$b_{\tau_{\|}}$, we have

$\displaystyle\frac{{\rm d}{\sigma_{{t{\bar{t}}}}}}{{\rm
d}M_{t{\bar{t}}}^{2}\,{\rm d}^{2}\vec{P}_{t}^{\perp}\,{\rm
d}Y_{t{\bar{t}}}\,{\rm d}q_{\tau}}\,\sim$
$\displaystyle\;\sum\limits_{m,n}\int^{\infty}_{-\infty}{\rm
d}b_{\tau_{\|}}\,\cos(b_{\tau_{\|}}q_{\tau})\,L^{n}_{\mathrm{M}}\,\Big{\\{}s_{m,n}(\beta_{t{\bar{t}}},x_{t},Y_{t{\bar{t}}})+a_{m,n}(\mathrm{sign}[b_{\tau_{\|}}],\beta_{t{\bar{t}}},x_{t},Y_{t{\bar{t}}})\Big{\\}}$
(3.2) $\displaystyle\,=$
$\displaystyle\;\sum\limits_{m,n}\Big{\\{}2s_{m,n}(\beta_{t{\bar{t}}},x_{t},Y_{t{\bar{t}}})+a_{m,n}^{+}(\beta_{t{\bar{t}}},x_{t},Y_{t{\bar{t}}})+a_{m,n}^{-}(\beta_{t{\bar{t}}},x_{t},Y_{t{\bar{t}}})\Big{\\}}\;\mathcal{F}^{(n)}_{\tau}(q_{\tau},M_{t{\bar{t}}})\,,$

where $a_{m,n}^{\pm}$ are the decomposition of the AAT components according to
$\mathrm{sign}[b_{\tau_{\|}}]$ and the function $\mathcal{F}^{(n)}_{\tau}$ is
defined as,

$\displaystyle\mathcal{F}^{(n)}_{\tau}(q_{\tau},M_{t{\bar{t}}})=\int^{\infty}_{0}{\rm
d}b_{\tau_{\|}}\,\cos(b_{\tau_{\|}}q_{\tau})\,L^{n}_{\mathrm{M}}\,.$ (3.3)

To appraise $\mathcal{F}^{(n)}_{\tau}(q_{\tau},M_{t{\bar{t}}})$, we follow the
strategy of [93] and introduce the generating function,

$\displaystyle\mathbf{F}_{\tau}(\eta,q_{\tau},M_{t{\bar{t}}})\,=$
$\displaystyle\;\sum_{n=0}^{\infty}\,\frac{\eta^{n}}{n!}\,\mathcal{F}^{(n)}_{\tau}(q_{\tau},M_{t{\bar{t}}})\,=\;\int^{\infty}_{0}{\rm
d}{b}_{\tau_{\|}}\,\cos\left({b}_{\tau_{\|}}{q}_{\tau}\right)\,\left(\frac{{b}^{2}_{\tau_{\|}}M_{t{\bar{t}}}^{2}}{b_{0}^{2}}\right)^{\eta}$
$\displaystyle=$
$\displaystyle\;-\frac{4^{-\eta}e^{2\gamma_{\mathrm{E}}\eta}\,\sin(\pi\eta)\,\Gamma[2\eta+1]}{{q}_{\tau}}\,\left(\frac{M_{t{\bar{t}}}^{2}}{q^{2}_{\tau}}\right)^{\eta}\,.$
(3.4)

Hence, $\mathcal{F}^{(n)}_{\tau}(q_{\tau},M_{t{\bar{t}}})$ corresponds to the
$n$-th derivative of $\mathbf{F}_{\tau}(\eta,q_{\tau},M_{t{\bar{t}}})$ with
respect to $\eta$ at the point $\eta=0$. The results of the first few ranks
read

$\begin{split}\mathcal{F}^{(0)}_{\tau}(q_{\tau},M_{t{\bar{t}}})\,=&\;\;0\,,\\\
\mathcal{F}^{(1)}_{\tau}(q_{\tau},M_{t{\bar{t}}})\,=&\;-\frac{\pi}{q_{\tau}}\,,\\\
\mathcal{F}^{(2)}_{\tau}(q_{\tau},M_{t{\bar{t}}})\,=&\;-\frac{2\pi}{q_{\tau}}\ln\left[\frac{M_{t{\bar{t}}}^{2}}{4q_{\tau}^{2}}\right]\,,\\\
\mathcal{F}^{(3)}_{\tau}(q_{\tau},M_{t{\bar{t}}})\,=&\;-\frac{3\pi}{q_{\tau}}\ln^{2}\left[\frac{M_{t{\bar{t}}}^{2}}{4q_{\tau}^{2}}\right]\,-\frac{\pi^{3}}{q_{\tau}},\\\
\mathcal{F}^{(4)}_{\tau}(q_{\tau},M_{t{\bar{t}}})\,=&\;-\frac{4\pi}{q_{\tau}}\ln^{3}\left[\frac{M_{t{\bar{t}}}^{2}}{4q_{\tau}^{2}}\right]\,-\frac{4\pi^{3}}{q_{\tau}}\ln\left[\frac{M_{t{\bar{t}}}^{2}}{4q_{\tau}^{2}}\right]\,+\frac{64\pi\zeta_{3}}{q_{\tau}},\\\
\vdots\;&\\\
\mathcal{F}^{(k)}_{\tau}(q_{\tau},M_{t{\bar{t}}})\,=&\;-\frac{k\pi}{q_{\tau}}\ln^{k-1}\left[\frac{M_{t{\bar{t}}}^{2}}{4q_{\tau}^{2}}\right]+\ldots\,.\end{split}$
(3.5)

Here we have focussed on the regime where $q_{\tau}$ is small, but always
larger than 0. Otherwise, contributions such as $\delta[q_{\tau}]$ will enter
the expressions above. From these results, it is observed that the
coefficients $s_{m,n}$ and $a_{m,n}^{\pm}$ cannot induce any divergent
behaviour by themselves as exhibited in Eq. (3.5), whilst the logarithmic
terms $L_{\mathrm{M}}^{k}$ produce the singular series up to
$\ln^{k-1}[{M_{t{\bar{t}}}^{2}}/{(4q_{\tau}^{2})}]/q_{\tau}$, i.e. the entire
LP asymptotic behaviors of Eq. (2.6). As a consequence, the resummation of the
$q_{\tau}$ spectrum can now be expressed as an exponentiation of the
$L_{\mathrm{M}}$, resembling the corresponding resummation in the Drell-Yan
processes and Higgs production [91, 92, 93, 94, 98, 99, 100, 87, 88, 89, 90].
In Sec. 3.2, we will use the solutions of RaGE and RGE to accomplish this
exponentiation.

It is interesting to note that, by analogy to the $q_{\tau}$ resummation, the
leading singular behaviour of the azimuthally averaged $q_{\mathrm{T}}$
distribution is also governed by the characteristic logarithmic terms from the
impact space [55, 56], such that the R(a)GE framework is in principle
applicable therein as well to accomplish the resummation. However, aside from
those two observables, the transverse momentum resummation is generally more
involved on the $pp\to{t{\bar{t}}}+X$ process. For instance, the singular
behaviour of the double differential observable ${\rm
d}\sigma_{{t{\bar{t}}}}/{\rm d}\vec{q}_{\mathrm{T}}$ is not only contained in
the logarithmic terms, but the AATs can also make up in part the asymptotic
series [58]. The appearance of these asymmetric divergent terms can have non-
trivial impacts on the pattern of the LP singularities and the choice of
resummation scheme. In App. A, we will deliver a comparative study on this
issue.

###  3.2 Evolution equations

We will now elaborate on the evolution equations of the hard, beam, and soft
functions introduced in Eqs. (2.29)-(2.30).

The hard function, see Eqs. (2.31)-(2.32), contains the squared UV-
renormalised amplitudes multiplied by the IRC regulator. Considering that the
scale dependences induced by the UV renormaliation cancel fully within the
amplitudes, the evolution equation of the hard sector is governed solely by
the IRC counterterm $\mathcal{Z}^{h}_{[\kappa]}$. According to the
parametrisation of [128, 161, 42], we have,

$\displaystyle\frac{{\rm d}}{{\rm
d}\ln\mu^{2}}\,\mathcal{H}^{[\kappa]}_{\alpha\beta}\big{(}M_{t{\bar{t}}}^{2},\beta_{t{\bar{t}}},x_{t},\mu\big{)}\,=$
$\displaystyle\;-C_{[\kappa]}\,\Gamma_{\mathrm{cusp}}\,\ln\left[\frac{\mu^{2}}{M_{t{\bar{t}}}^{2}}\right]\,\mathcal{H}^{[\kappa]}_{\alpha\beta}\big{(}M_{t{\bar{t}}}^{2},\beta_{t{\bar{t}}},x_{t},\mu\big{)}\,$
$\displaystyle\;{}+\sum_{\gamma}\left[\frac{\mathbf{\gamma}^{[\kappa]}_{h,\beta\gamma}}{2}\,\mathcal{H}^{[\kappa]}_{\alpha\gamma}\big{(}M_{t{\bar{t}}}^{2},\beta_{t{\bar{t}}},x_{t},\mu\big{)}+\left(\frac{\mathbf{\gamma}^{[\kappa]}_{h,\alpha\gamma}}{2}\right)^{\\!\\!\\!*}\mathcal{H}^{[\kappa]}_{\gamma\beta}\big{(}M_{t{\bar{t}}}^{2},\beta_{t{\bar{t}}},x_{t},\mu\big{)}\right]\,,$
(3.6)

where the subscripts $\\{\alpha,\beta,\gamma\\}$ represent the colour indices
in the set of colour basis in Eqs. (2.33)-(2.34). The helicity indices
involved in the gluonic channel have been omitted for brevity. $\kappa$ runs
over
$\\{g_{n}g_{{\bar{n}}},q^{i}_{n}\bar{q}^{j}_{{\bar{n}}},q^{i}_{{\bar{n}}}\bar{q}^{j}_{n}\\}$,
indicating the partonic channel. $C_{[\kappa]}$ denotes the colour factor in
QCD with

$\displaystyle\kappa\in\\{g_{n}g_{{\bar{n}}}\\}\,:\;C_{[\kappa]}=C_{A}\;,\qquad\qquad\kappa\in\\{q^{i}_{n}\bar{q}^{j}_{{\bar{n}}},q^{i}_{{\bar{n}}}\bar{q}^{j}_{n}\\}\,:\;C_{[\kappa]}=C_{F}\,.$
(3.7)

$\Gamma_{\mathrm{cusp}}$ is the cusp anomalous dimension. It is needed to
three-loop accuracy [162] in this paper, but is available to four-loop
precision in the literature [163, 164]. A numeric estimation of the five-loop
contribution is addressed in [165]. Analogously, $\gamma^{[\kappa]}_{h}$ is
the non-cusp anomalous dimension for the hard contribution. Their analytic
expressions up to N2LO can be found in [128, 161] and progress towards the
three-loop result has been made in [166].

The quark and gluon beam functions are given in Eqs. (2.35)-(2.36) in terms of
the SCETII field operators. They are same as those participating into the
Drell-Yan processes and Higgs production. They admit the RGEs [89],

$\displaystyle\frac{\partial}{\partial\ln\mu^{2}}\ln\mathcal{B}_{n}^{[q_{n}^{i}]}(\eta_{n},b_{\mathrm{T}},\mu,\nu)\,=$
$\displaystyle\;C_{F}\,\Gamma_{\mathrm{cusp}}\,\ln\bigg{[}\frac{\nu}{\eta_{n}\,\sqrt{s}}\bigg{]}+\gamma^{[q]}_{b}\;,$
(3.8)
$\displaystyle\frac{\partial}{\partial\ln\mu^{2}}\ln\mathcal{B}^{[g_{n}]}_{n,h^{\prime}h}(\eta_{n},\vec{b}_{\mathrm{T}},\mu,\nu)\,=$
$\displaystyle\;C_{A}\,\Gamma_{\mathrm{cusp}}\,\ln\bigg{[}\frac{\nu}{\eta_{n}\,\sqrt{s}}\bigg{]}+\gamma^{[g]}_{b}\;,$
(3.9)

as well as the RaGEs [89],

$\displaystyle-\frac{2}{C_{F}}\,\frac{\partial}{\partial\ln\nu^{2}}\ln\mathcal{B}_{n}^{[q_{i}]}(\eta_{n},b_{\mathrm{T}},\mu,\nu)\,=$
$\displaystyle\;-\frac{2}{C_{A}}\,\frac{\partial}{\partial\ln\nu^{2}}\ln\mathcal{B}^{[g_{n}]}_{n,h^{\prime}h}(\eta_{n},\vec{b}_{\mathrm{T}},\mu,\nu)$
$\displaystyle\,=$
$\displaystyle\;\gamma_{r}\left[\alpha_{s}\left(\frac{b_{0}}{b_{\mathrm{T}}}\right)\right]\,+\int_{\mu^{2}}^{\frac{b_{0}^{2}}{b_{\mathrm{T}}^{2}}}\,\frac{{\rm
d}\bar{\mu}^{2}}{\bar{\mu}^{2}}\,\Gamma_{\mathrm{cusp}}\left[\alpha_{s}(\bar{\mu})\right]\,.$
(3.10)

Here, $\eta_{n}$ stands for the momentum fraction along the $n$ direction.
$\gamma^{[q,g]}_{b}$ and $\gamma_{r}$ are the non-cusp anomalous dimensions
brought about in the virtuality and rapidity renormalisations, respectively.
Their specific expressions are dependent on the choice of regularisation
prescription, and this work will use those that correspond with our choice of
using the exponential rapidity regulator [89]. They are known at N2LO accuracy
[136, 135], which we use in the following. N3LO results are also available in
the literature [90, 167, 137, 168, 138], and N4LO corrections [169, 170, 171,
172] have appeared recently.

In order to derive the evolution equations for the soft sector, we utilise the
scale invariance condition in Eq. (2.81) and the R(a)GEs above. It follows
that

$\displaystyle\frac{\partial}{\partial\ln\mu^{2}}\mathcal{S}_{[\kappa]}^{\alpha\beta}\big{(}\vec{b}_{\mathrm{T}},\mu,\nu\big{)}\,=$
$\displaystyle\;-\bigg{\\{}C_{[\kappa]}\,\Gamma_{\mathrm{cusp}}\,\ln\left[\frac{\nu^{2}}{\mu^{2}}\right]+2\gamma_{b}^{[\kappa]}\bigg{\\}}\,\mathcal{S}_{[\kappa]}^{\alpha\beta}\big{(}\vec{b}_{\mathrm{T}},\mu,\nu\big{)}$
$\displaystyle\;-\sum_{\gamma}\,\left[\mathcal{S}_{[\kappa]}^{\alpha\gamma}\big{(}\vec{b}_{\mathrm{T}},\mu,\nu\big{)}\,\frac{\mathbf{\gamma}^{[\kappa]}_{h,\gamma\beta}}{2}+\mathcal{S}_{[\kappa]}^{\gamma\beta}\big{(}\vec{b}_{\mathrm{T}},\mu,\nu\big{)}\,\left(\frac{\mathbf{\gamma}^{[\kappa]}_{h,\gamma\alpha}}{2}\right)^{\\!\\!\\!*}\,\right]\,,$
(3.11)

and

$\displaystyle\frac{1}{C_{[\kappa]}}\,\frac{\partial}{\partial\ln\nu^{2}}\ln\mathcal{S}^{\alpha\beta}_{[\kappa]}\big{(}\vec{b}_{\mathrm{T}},\mu,\nu\big{)}\,=$
$\displaystyle\;\gamma_{r}\left[\alpha_{s}\left(\frac{b_{0}}{b_{\mathrm{T}}}\right)\right]\,+\int_{\mu^{2}}^{\frac{b_{0}^{2}}{b_{\mathrm{T}}^{2}}}\,\frac{{\rm
d}\bar{\mu}^{2}}{\bar{\mu}^{2}}\,\Gamma_{\mathrm{cusp}}\left[\alpha_{s}(\bar{\mu})\right]\,.$
(3.12)

For brevity, we make use of $\gamma_{b}^{[\kappa]}$ here to represent the
virtuality anomalous dimensions in Eqs. (3.8)-(3.9), more specifically,

$\kappa\in\\{g_{n}g_{{\bar{n}}}\\}\,:\;\gamma_{b}^{[\kappa]}=\gamma_{b}^{[g]}\,,\qquad\qquad\kappa\in\\{q^{i}_{n}\bar{q}^{j}_{{\bar{n}}},q^{i}_{{\bar{n}}}\bar{q}^{j}_{n}\\}\,:\;\gamma_{b}^{[\kappa]}=\gamma_{b}^{[q]}\,.$
(3.13)

As observed in Eqs. (3.11)-(3.12), since the anomalous dimensions herein are
all independently extracted from the hard and beam functions, these evolution
equations provide a non-trivial opportunity to examine our soft function of
Sec. 2.5 and in turn the factorisation in Eqs. (2.29)-(2.30)). Substituting
the expressions of Eq. (2.85) into Eqs. (3.11)-(3.12), we have checked that
our results indeed satisfy the criteria above on the scale dependences.

Solving those RGEs and RaGEs permits us to bridge the intrinsic scales of the
hard, beam, and soft ingredients, thereby exponentiating the characteristic
logarithmic terms in Eq. (3.1). Substituting these solutions into the master
formula of Eq. (2.43), we arrive at the resummed $q_{\tau}$ spectrum,

$\displaystyle\frac{{\rm d}{\sigma^{\mathrm{res}}_{{t{\bar{t}}}}}}{{\rm
d}M_{t{\bar{t}}}^{2}\,{\rm d}^{2}\vec{P}_{t}^{\perp}\,{\rm
d}Y_{t{\bar{t}}}\,{\rm d}q_{\tau}}\,=$
$\displaystyle\;\sum_{\mathrm{sign}[x_{t}]}\frac{\Theta_{\mathrm{kin}}^{(0)}}{16\pi^{3}\,\beta_{t{\bar{t}}}\,|x_{t}|\,M^{4}_{{t{\bar{t}}}}\,s}\sum_{\kappa}\int^{\infty}_{-\infty}{\rm
d}{b}_{\tau_{\|}}\,\cos\left({b}_{\tau_{\|}}{q}_{\tau}\right)\,\widetilde{\Sigma}_{{t{\bar{t}}}}^{\mathrm{res},[\kappa]}(b_{\tau_{\|}}\vec{\tau},M_{t{\bar{t}}},\beta_{t{\bar{t}}},x_{t},Y_{t{\bar{t}}})\,,$
(3.14)

where

$\displaystyle\widetilde{\Sigma}_{{t{\bar{t}}}}^{\mathrm{res},[q^{i}_{n}\bar{q}^{j}_{{\bar{n}}}]}(\vec{b}_{\mathrm{T}},M_{t{\bar{t}}},\beta_{t{\bar{t}}},x_{t},Y_{t{\bar{t}}})$
(3.15) $\displaystyle\,=$
$\displaystyle\;\left(\frac{1}{2N_{c}}\right)^{2}\,\mathcal{D}^{\mathrm{res}}_{[q^{i}_{n}\bar{q}^{j}_{{\bar{n}}}]}(b_{\mathrm{T}},M_{t{\bar{t}}},\mu_{h},\mu_{b},\mu_{s},\nu_{b},\nu_{s})\,\times\mathcal{B}_{n}^{[q_{n}^{i}]}(\eta_{n},b_{\mathrm{T}},\mu_{b},\nu_{b})\,\mathcal{B}_{{\bar{n}}}^{[\bar{q}_{\bar{n}}^{j}]}(\eta_{{\bar{n}}},b_{\mathrm{T}},\mu_{b},\nu_{b})\,$
$\displaystyle\sum_{\\{\alpha,\beta\\}}\Bigg{\\{}\mathcal{S}^{\alpha_{1}\beta_{1}}_{[q_{n}\bar{q}_{{\bar{n}}}]}(\vec{b}_{\mathrm{T}},\mu_{s},\nu_{s})\left[\mathcal{V}^{[q_{n}\bar{q}_{{\bar{n}}}]}_{\alpha_{1}\alpha_{2}}(\beta_{t{\bar{t}}},x_{t},\mu_{s},\mu_{h})\right]^{*}\,\mathcal{V}^{[q_{n}\bar{q}_{{\bar{n}}}]}_{\beta_{1}\beta_{2}}(\beta_{t{\bar{t}}},x_{t},\mu_{s},\mu_{h})\,\mathcal{H}_{\alpha_{2}\beta_{2}}^{[q^{i}_{n}\bar{q}^{j}_{{\bar{n}}}]}(M_{t{\bar{t}}},\beta_{t{\bar{t}}},x_{t},\mu_{h})\,\Bigg{\\}}\,,$
$\displaystyle\widetilde{\Sigma}_{{t{\bar{t}}}}^{\mathrm{res},[g_{n}g_{{\bar{n}}}]}(\vec{b}_{\mathrm{T}},M_{t{\bar{t}}},\beta_{t{\bar{t}}},x_{t},Y_{t{\bar{t}}})$
(3.16) $\displaystyle\,=$
$\displaystyle\;\left(\frac{1}{N^{2}_{c}-1}\right)^{2}\,\mathcal{D}^{\mathrm{res}}_{[g_{n}g_{{\bar{n}}}]}(b_{\mathrm{T}},M_{t{\bar{t}}},\mu_{h},\mu_{b},\mu_{s},\nu_{b},\nu_{s})\,\sum_{\\{\alpha,\beta,h,h^{\prime}\\}}\Bigg{\\{}\mathcal{S}^{\alpha_{1}\beta_{1}}_{[g_{n}g_{{\bar{n}}}]}(\vec{b}_{\mathrm{T}},\mu_{s},\nu_{s})\,$
$\displaystyle\times\,\mathcal{B}^{[g_{n}]}_{n,h_{n}^{\prime}{h_{n}}}(\eta_{n},\vec{b}_{\mathrm{T}},\mu_{b},\nu_{b})\,\mathcal{B}^{[g_{{\bar{n}}}]}_{{\bar{n}},h_{{\bar{n}}}^{\prime}{h_{{\bar{n}}}}}(\eta_{{\bar{n}}},\vec{b}_{\mathrm{T}},\mu_{b},\nu_{b})\,\left[\mathcal{V}^{[g_{n}g_{{\bar{n}}}]}_{\alpha_{1}\alpha_{2}}(\beta_{t{\bar{t}}},x_{t},\mu_{s},\mu_{h})\right]^{*}\,\mathcal{V}^{[g_{n}g_{{\bar{n}}}]}_{\beta_{1}\beta_{2}}(\beta_{t{\bar{t}}},x_{t},\mu_{s},\mu_{h})\,$
$\displaystyle\times\mathcal{H}_{\alpha_{2}\beta_{2};{h_{{\bar{n}}}^{\prime}}h_{{\bar{n}}};{h_{n}^{\prime}}h_{n}}^{[g_{n}g_{{\bar{n}}}]}(M_{t{\bar{t}}},\beta_{t{\bar{t}}},x_{t},\mu_{h})\,\Bigg{\\}}\,.$

The expression of
$\widetilde{\Sigma}_{{t{\bar{t}}}}^{\mathrm{res},[q^{i}_{{\bar{n}}}\bar{q}^{j}_{n}]}$
can be obtained from Eq. (3.15) by adjusting the label momenta of the beam-
collinear modes as appropriate. As is apparent in Eqs. (3.15-3.16), two sets
of auxiliary scales $\\{\mu_{h},\mu_{b},\mu_{s}\\}$ and
$\\{\nu_{b},\nu_{s}\\}$ have been introduced to define the initial conditions
utilised in solving the RaGEs and RGEs above. An appropriate choice of their
values minimises the missing higher-order corrections, and in this paper, in a
bid to minimise the logarithmic dependences on the respective sectors, the
following values will be taken as defaults in this paper [88, 173],

$\displaystyle\mu^{\mathrm{def}}_{h}=\nu^{\mathrm{def}}_{b}=M_{t{\bar{t}}}\,,\quad\mu^{\mathrm{def}}_{b}=\mu^{\mathrm{def}}_{s}=\nu^{\mathrm{def}}_{s}=b_{0}/|b_{\tau_{\|}}|\,.$
(3.17)

With the choices in Eq. (3.17), the impact space integral in Eq. (3.14) may
approach or cross the Landau singularity in the large $|b_{\tau_{\|}}|$
regime. In order to avoid the divergence, we impose upper and lower boundaries
$|b_{\tau_{\|}}|\leq b^{\mathrm{cut}}_{\tau_{\|}}=2\,\text{GeV}^{-1}$ in Eq.
(3.14) [173]. Alternative schemes have been discussed in [98, 174].

Further, Eqs. (3.15)-(3.16) also include the kernels
$\mathcal{D}^{\mathrm{res}}_{[\kappa]}$ and
$\mathcal{V}^{[\kappa]}_{\alpha\beta}$ to evolve the intrinsic scales amongst
the fixed-order functions. $\mathcal{D}^{\mathrm{res}}_{[\kappa]}$ is induced
by the diagonal anomalous dimensions, such as $\gamma_{\mathrm{cusp}}$,
$\gamma_{b}^{[\kappa]}$, and $\gamma_{r}$ in Eqs. (3.2)-(3.10). Its definition
reads

$\displaystyle\ln\mathcal{D}^{\mathrm{res}}_{[\kappa]}(b_{\mathrm{T}},M_{t{\bar{t}}},\mu_{h},\mu_{b},\mu_{s},\nu_{b},\nu_{s})$
$\displaystyle\,=\;\int^{\mu^{2}_{s}}_{\mu^{2}_{b}}\,\frac{{\rm
d}\bar{\mu}^{2}}{\bar{\mu}^{2}}\,\Bigg{\\{}C_{[\kappa]}\,\Gamma_{\mathrm{cusp}}\big{[}\alpha_{s}(\bar{\mu})\big{]}\,\ln\bigg{[}\frac{\nu_{b}^{2}}{M_{t{\bar{t}}}^{2}}\bigg{]}+2\,\gamma^{[\kappa]}_{b}\big{[}\alpha_{s}(\bar{\mu})\big{]}\,\Bigg{\\}}-\int^{\mu_{s}^{2}}_{\mu_{h}^{2}}\frac{{\rm
d}\bar{\mu}^{2}}{\bar{\mu}^{2}}\,\bigg{\\{}{C_{[\kappa]}}\,\Gamma_{\mathrm{cusp}}\big{[}\alpha_{s}(\bar{\mu})\big{]}\,\ln\left[\frac{\bar{\mu}^{2}}{M_{t{\bar{t}}}^{2}}\right]\,\bigg{\\}}$
$\displaystyle\,\phantom{=}\;{}+C_{[\kappa]}\,\ln\bigg{[}\frac{\nu_{s}^{2}}{\nu_{b}^{2}}\bigg{]}\,\int_{\frac{b_{0}^{2}}{b_{\mathrm{T}}^{2}}}^{\mu_{s}^{2}}\,\frac{{\rm
d}\bar{\mu}^{2}}{\bar{\mu}^{2}}\,\Gamma_{\mathrm{cusp}}\big{[}\alpha_{s}(\bar{\mu})\big{]}\,-C_{[\kappa]}\,\ln\bigg{[}\frac{\nu_{s}^{2}}{\nu_{b}^{2}}\bigg{]}\,\gamma_{r}\bigg{[}\alpha_{s}\left(\frac{b_{0}}{b_{\mathrm{T}}}\right)\bigg{]}\,.$
(3.18)

$\mathcal{V}^{[\kappa]}_{\alpha\beta}$ accounts for the contributions from the
non-diagonal anomalous dimension $\gamma_{h,\alpha\beta}^{[\kappa]}$, which,
in principle, can be extracted from the solutions of the RGE of the hard
function in Eq. (3.2). However, in the presence of the non-diagonal elements
in $\gamma_{h}^{[\kappa]}$, it is somewhat challenging to achieve a closed
expression. Hence, in this work we adopt the perturbative approaches suggested
in [175, 176, 42]. The details on the implementation are collected in App. B.

Logarithmic accuracy$\int_{A}^{B}$ |  $\mathcal{H}$, $\mathcal{S}$, $\mathcal{B}$ | $\Gamma_{\text{cusp}}$ | $\gamma_{h,s,b,r}$
---|---|---|---
NLL′ $\int_{A}^{B}$ | $\mathcal{O}(\alpha^{0}_{s})$ | $\mathcal{O}(\alpha^{2}_{s})$ | $\mathcal{O}(\alpha_{s})$
N2LL $\int_{A}^{B}$ | $\mathcal{O}(\alpha_{s})$ | $\mathcal{O}(\alpha^{3}_{s})$ | $\mathcal{O}(\alpha^{2}_{s})$
N2LL′ $\int_{A}^{B}$ | $\mathcal{O}(\alpha^{2}_{s})$ | $\mathcal{O}(\alpha^{3}_{s})$ | $\mathcal{O}(\alpha^{2}_{s})$
Table 1:  Precision prerequisites on the anomalous dimensions and the fixed-