# NNLL Resummation of Sudakov Shoulder Logarithms in the Heavy Jet Mass
Distribution

Arindam Bhattacharya  Johannes K. L. Michel  Matthew D. Schwartz
Iain W. Stewart  and Xiaoyuan Zhang

(June 13, 2023)

###### Abstract

The heavy jet mass event shape has large perturbative logarithms near the
leading order kinematic threshold at $\rho=\frac{1}{3}$. Catani and Webber
named these logarithms Sudakov shoulders and resummed them at double-
logarithmic level. A resummation to next-to-leading logarithmic level was
achieved recently. Here, we extend the resummation using an effective field
theory framework to next-to-next-to-leading logarithmic order and show how to
combine it with the resummation of dijet logarithms. We also solve the open
problem of an unphysical singularity in the resummed momentum space
distribution, in a way similar to how it is resolved in the Drell-Yan $q_{T}$
spectrum: through a careful analysis of the kinematics and scale-setting in
position space. The heavy jet mass Sudakov shoulder is the first observable
that does not involve transverse momentum for which position space resummation
is critical. These advances may lead to a more precise extraction of the
strong coupling constant from $e^{+}e^{-}$ data.

##  1 Introduction

Event shapes in electron-positron collisions have long been valuable testing
grounds for quantum field theory and quantum chromodynamics. High quality
event shape data predominantly from the Large Electron-Positron Collider (LEP)
have allowed for a detailed comparison between theory and experiment and
engendered improvements on both fronts. On the theory side, inspiration from
data has lead to higher-precision fixed-order computations, up to the current
state of the art of NNLO ($\alpha_{s}^{3}$) [1, 2, 3, 4, 5, 6]. Data has also
inspired exponential progress in resummation, up to N3LL precision using Soft-
Collinear Effective Theory (SCET) [7, 8, 9, 10, 11] in the dijet region for
thrust, heavy jet mass, and the $C$ parameter [12, 13, 14, 15], and NNLL′
precision on other event shapes [16]. The comparison between theory and data
has also motivated improvements in the treatment of non-perturbative
corrections from hadronization, including using matrix elements defined in
quantum field theory [17, 18, 19], an improved understanding of universality
relations [20, 21], and the inclusion of hadron mass corrections [22, 23].

Until recently, however, there has been very little understanding of the
trijet region, with 3 hard jets, in $e^{+}e^{-}$ events. This region has
traditionally been assumed to be well-modeled in fixed-order perturbation
theory, since it is well-separated from the dijet region where large threshold
logarithms spoil the convergence of the perturbation expansion in
$\alpha_{s}$. However, fixed-order perturbation theory fails spectacularly
right at the 3-jet threshold where it predicts kinks and discontinuities. The
failure of perturbation theory in the trijet region was understood by Catani
and Webber [24], and also Seymour [25], as due to the appearance of Sudakov
shoulder logarithms. Although the resummation of the Sudakov shoulder
logarithms at the double-logarithmic level was achieved analytically for the
$C$ parameter, little additional progress was made in the last 25 years.
Recently, Sudakov shoulder resummation was revisited using SCET where a
factorization theorem was derived and next-to-leading logarithmic (NLL)
resummation was achieved for thrust and heavy jet mass [26]. In addition,
there has been recent progress in understanding power corrections near the
3-jet region using renormalon methods [27, 28] and the dispersive approach
[29]. In this paper, we focus on understanding and improving perturbative
Sudakov shoulder resummation for heavy jet mass. Heavy jet mass is special in
that it has a Sudakov shoulder on the left of the trijet threshold, so the
shoulder logarithms propagate into the region usually used for fits to
$\alpha_{s}$. In contrast, thrust only has a right shoulder, so that large
shoulder logarithms do not impact the region used for such fits.

Thrust is defined in the center-of-mass frame of an $e^{+}e^{-}$ collision as

$T\equiv\max_{\vec{n}}\frac{\sum_{j}|\vec{p}_{j}\cdot\vec{n}|}{\sum_{j}|\vec{p}_{j}|}\,,$
(1.1)

where the sum is over all particles in the event and the maximum is over
3-vectors $\vec{n}$ of unit norm. The vector $\vec{n}$ that maximizes thrust
is known as the thrust axis. The thrust axis splits the event into two
hemispheres. One can then add all the momenta in each hemisphere and compute
the hemisphere invariant masses $M_{L}$ and $M_{R}$. Heavy jet mass is defined
as

$\rho=\frac{1}{Q^{2}}\max(M_{L}^{2},M_{R}^{2})\,.$ (1.2)

Conventionally one also defines $\tau=1-T$. So at zeroth order with two
massless partons $\tau=\rho=0$. We call the region near $\tau=\rho=0$ the
dijet region. At leading order in $\alpha_{s}$ the thrust and heavy jet mass
distributions are the same (for massless particles) with a kinematic maximum
at $\tau=\rho=\frac{1}{3}$.

Both thrust and heavy jet mass have non-analytic behavior near
$\tau=\rho=\frac{1}{3}$. To be concrete, the leading non-analytic behavior of
thrust near $\tau=\frac{1}{3}$ at NLO was computed in [26] and found to be in

$\frac{1}{\sigma_{0}}\frac{d\sigma}{d\tau}\approx\frac{\alpha_{s}}{4\pi}48C_{F}\,\left[3-\left(\frac{\alpha_{s}}{4\pi}\right)\frac{3}{2}(2C_{F}+C_{A})\Gamma_{0}\ln^{2}t+\cdots\right]t\theta(t)\,,$
(1.3)

where $t\equiv\tau-\frac{1}{3}$, with $\sigma_{0}$ denoting the Born cross
section. Thrust only has a right shoulder: its behavior for
$\tau\lesssim\frac{1}{3}$ is smooth. To get the full leading order thrust
distribution in this region one must add to Eq. (1.3) the analytic linear
function $-\frac{\alpha_{s}}{4\pi}144C_{F}t$ which effectively turns
$t\theta(t)$ into $(-t)\theta(-t)$ with support only for $\tau<\frac{1}{3}$.
For heavy jet mass near $\rho=\frac{1}{3}$ we define $r=\frac{1}{3}-\rho$. The
leading non-analytic behavior of its distribution near $r=0$ is

$\displaystyle\frac{1}{\sigma_{0}}\frac{d\sigma}{d\rho}\approx\frac{\alpha_{s}}{4\pi}\,48C_{F}\,\Bigl{[}3$
$\displaystyle-\Bigl{(}\frac{\alpha_{s}}{4\pi}\Bigr{)}\frac{1}{2}(2C_{F}+C_{A})\Gamma_{0}\ln^{2}r+\cdots\Bigr{]}\,r\,\theta(r)$
$\displaystyle+\Bigl{(}\frac{\alpha_{s}}{4\pi}\Bigr{)}^{2}\,48C_{F}\Bigl{[}-(2C_{F}+C_{A})\Gamma_{0}\ln^{2}(-r)+\cdots\Bigr{]}(-r)\,\theta(-r)\,.$
(1.4)

Heavy jet mass has large logarithms on both sides of $\rho=\frac{1}{3}$,
giving it a left $(r>0)$ and right ($r<0$) Sudakov shoulder.

Physically, the origin of the large logarithms can be understood as follows.
At leading order (three partons) there is only a single phase space point with
$\rho=\frac{1}{3}$, namely the “trijet configuration” where the partons have
energy $\frac{Q}{3}$ and are spaced $\frac{2\pi}{3}$ apart; two partons are
clustered into one hemisphere and one parton remains in the other. The thrust
axis can be aligned with any of the partons. As these three partons emit soft
and collinear gluons, the cross section becomes logarithmically enhanced.
These soft and collinear logarithms in the cross section in the trijet region
become Sudakov shoulder logarithms after the thrust or heavy jet mass
measurement function is applied.

A careful analysis of soft and collinear kinematics near the trijet region
leads to a factorization formula for Sudakov shoulders. We can treat collinear
emissions as transforming the three massless partons into jets with masses
$M_{1,2,3}$. The thrust axis will align with one of the jets, but which one it
aligns with depends on the exact kinematics. If $M_{1}$ is the mass of the jet
in the light hemisphere, and $M_{2}$ and $M_{3}$ are the masses of the jets in
the heavy hemisphere, then, defining $m_{i}=\frac{M_{i}}{Q}$, the bound is

$m_{1}^{2}<r+m_{2}^{2}+m_{3}^{2}\,.$ (1.5)

As $m_{1}$ gets larger or as $m_{2}$ or $m_{3}$ gets smaller, the thrust axis
flips and a different parton is in the light hemisphere. Soft contributions
can be analyzed similarly. The soft and collinear emissions into the light and
heavy hemispheres can be grouped into effective light-hemisphere and heavy-
hemisphere mass-shift variables $m_{\ell}^{2}\approx m^{2}_{1}$ and
$m_{h}^{2}\approx m_{2}^{2}+m_{3}^{2}$ and the factorization formula for heavy
jet mass written as

$\frac{d\sigma}{d\rho}=\int_{0}^{\infty}dm_{h}^{2}\int_{0}^{\infty}dm_{\ell}^{2}\frac{d^{2}\sigma}{dm_{\ell}^{2}dm_{h}^{2}}(r+m_{h}^{2}-m_{\ell}^{2})\theta(r+m_{h}^{2}-m_{\ell}^{2})\,.$
(1.6)

The double differential distribution
$\frac{d^{2}\sigma}{dm_{\ell}^{2}dm_{h}^{2}}$ can then be resummed using
standard techniques, see [30].

A key feature of the heavy jet mass Sudakov shoulder factorization formula in
Eq. (1.6) that distinguishes it from typical event shape factorization
formulas is the relative minus sign between $m_{h}^{2}$ and $m_{\ell}^{2}$. In
particular, for small $r$, the region where $m_{\ell}^{2}\gg 1$ and
$m_{h}^{2}\gg 1$ contributes, even though this region is not described by soft
and collinear physics. Such contributions could generically lead to non-global
logarithms, or worse, a failure of the factorization theorem to describe the
threshold behavior. However, in this case, as explained in [26], we still
expect all the shoulder logs to be described correctly. The reason is that
when $m_{h}$ and $m_{\ell}$ are large, the point $r=0$ is not special:
contributions to this region should go smoothly from $r<0$ to $r>0$. The kinks
and large logarithms in the shoulder region emerge only when $m_{h}$ and
$m_{\ell}$ are small so that the lower bounds $m_{h}>0$ and $m_{\ell}>0$
become relevant. Despite this, when one uses standard resummation methods with
canonical scale choices, spurious “Sudakov Landau pole” singularities result.
Similar poles have been seen also in the small $q_{T}$ Drell-Yan spectrum [31,
32, 33, 34, 35, 36]. A main result of this paper is a demonstration that these
Sudakov Landau poles are avoided by carrying out the resummation with
canonical scales in position space. Additional results are the extension of
the Sudakov shoulder resummation to NNLL, and smoothly combining the dijet
resummation, shoulder resummation, and fixed-order corrections for the heavy
jet mass.

An outline of the paper is as follows. In Section 2, we review the
factorization theorem for the Sudakov shoulder region, and describe the
necessary ingredients for NNLL resummation. Following that, in Section 3, we
review the appearance of spurious Sudakov Landau poles when setting scales in
momentum ($r$) space, and address how these can be carefully removed by
setting scales in position (Fourier) space. In Section 4, we describe our
resummation procedure for Sudakov shoulder logs, and set up a matching
procedure to account for dijet logs as well, thereby forming a consistent
picture of resummation across both the dijet ($\rho\gtrsim 0$) and trijet
region ($\rho\sim 1/3$) in heavy jet mass. In Section 5, we provide our
numerical results which combine NNLL shoulder resummation and NNLL′ dijet
resummation matched to NLO. Finally, we conclude in Section 6.

##  2 Factorization formula

In this section we review the factorization formula from [26]. In the trijet
limit, there are three independent channels associated with a gluon jet in the
light hemisphere, a quark jet in the light hemisphere or an anti-quark jet in
the light hemisphere. The contribution of these channels to the total cross
section add incoherently, leading to a form

$\frac{d\sigma_{\text{sh}}}{d\rho}=\sum_{i}\frac{d\sigma_{i}}{d\rho}=\frac{d\sigma_{g}}{d\rho}+2\frac{d\sigma_{q}}{d\rho}\,,$
(2.1)

where we have used symmetry under charge-conjugation to identify the quark and
antiquark channels. In our notation, when $\sigma$ has a $g$ or $q$ subscript
it automatically identifies it as a leading power factorized contribution near
the shoulder.

###  2.1 Factorization formula

We first focus on the gluon channel for concreteness. In [26], it was shown
that the factorization formula for heavy jet mass in the shoulder region has
the form

$\displaystyle\frac{1}{\sigma_{\text{LO}}}\frac{d\sigma_{g}}{d\rho}$
$\displaystyle=\Pi_{g}\,\int_{0}^{\infty}dm_{\ell}^{2}\int_{0}^{\infty}dm_{h}^{2}\,\left.\widetilde{\phantom{\,}sj}_{g}\cdot\left(\frac{Qm_{l}^{2}}{\mu_{s\ell}}\right)^{a}\left(\frac{Qm_{h}^{2}}{\mu_{sh}}\right)^{b}\frac{1}{m_{l}^{2}m_{h}^{2}}\frac{e^{-\gamma_{E}(a+b)}}{\Gamma(a)\Gamma(b)}\right|_{\overset{\scriptstyle
a=\eta_{\ell}\mathstrut}{\scriptstyle b=\eta_{h}}}$
$\displaystyle\qquad\times(r+m_{h}^{2}-m_{\ell}^{2})\theta(r+m_{h}^{2}-m_{\ell}^{2})\,,$
(2.2)

with $\sigma_{\text{LO}}=\frac{\alpha_{s}}{4\pi}48C_{F}\sigma_{0}$, where
$\sigma_{0}$ is the Born cross-section for $e^{+}e^{-}\to q\bar{q}$. Here
$m_{h}^{2}$ and $m_{\ell}^{2}$ are the (dimensionless) mass-shift variables,
representing the change in the heavy or light jet masses-squared due to soft
and collinear emissions relative to $Q^{2}$; the evolution kernel is,

$\displaystyle\Pi_{g}$
$\displaystyle=\exp\Big{[}4C_{F}\left(S(\mu_{h},\mu_{jh})+S(\mu_{sh},\mu_{jh})\right)+2C_{A}\left(S(\mu_{h},\mu_{j\ell})+S(\mu_{sl},\mu_{j\ell})\right)\Big{]}$
$\displaystyle\quad\times\exp\Big{[}2A_{\gamma
sg}(\mu_{s\ell},\mu_{h})+2A_{\gamma sqq}(\mu_{sh},\mu_{h})+2A_{\gamma
jg}(\mu_{j\ell},\mu_{h})+4A_{\gamma jq}(\mu_{jh},\mu_{h})\Big{]}$
$\displaystyle\quad\times
H(Q,\mu_{h})\left(\frac{Q^{2}}{\mu_{h}^{2}}\right)^{-2C_{F}A_{\Gamma}(\mu_{h},\mu_{jh})-C_{A}A_{\Gamma}(\mu_{h},\mu_{j\ell})}\,,$
(2.3)

the fixed-order matching operator is

$\widetilde{\phantom{\,}sj}_{g}=\widetilde{j}_{q}\left(\partial_{b}+\ln\frac{Q\mu_{sh}}{\mu_{jh}^{2}}\right)\widetilde{j}_{\bar{q}}\left(\partial_{b}+\ln\frac{Q\mu_{sh}}{\mu_{jh}^{2}}\right)\widetilde{j}_{g}\left(\partial_{a}+\ln\frac{Q\mu_{s\ell}}{\mu_{j\ell}^{2}}\right)\widetilde{s}_{q\bar{q};g}(\partial_{a},\partial_{b})\,,$
(2.4)

and after taking derivatives, $a$ and $b$ are set to

$a=\eta_{\ell}=2C_{A}A_{\Gamma}(\mu_{j\ell},\mu_{s\ell}),\qquad
b=\eta_{h}=4C_{F}A_{\Gamma}(\mu_{jh},\mu_{sh})\,.$ (2.5)

Here we have allowed for different jet and soft factorization scales in the
heavy hemisphere ($\mu_{jh}$, $\mu_{sh}$) and light hemisphere ($\mu_{j\ell}$,
$\mu_{s\ell}$), slightly generalizing the expression in [26]. As we discuss in
Section 3.4 below, this separation is not sensible beyond NLL, and ultimately
we will set $\mu_{s\ell}$ and $\mu_{sh}$ to a common scale $\mu_{s}$. For the
quark channel the evolution kernel $\Pi_{g}$ is replaced by

$\displaystyle\Pi_{q}=\exp\Big{[}(2C_{F}+2C_{A})\left(S(\mu_{h},\mu_{jh})+S(\mu_{sh},\mu_{jh})\right)+2C_{F}\left(S(\mu_{h},\mu_{j\ell})+S(\mu_{sl},\mu_{j\ell})\right)\Big{]}$
$\displaystyle\times\exp\Big{[}2A_{\gamma sq}(\mu_{s\ell},\mu_{h})+2A_{\gamma
sqg}(\mu_{sh},\mu_{h})+2A_{\gamma jq}(\mu_{j\ell},\mu_{h})+2A_{\gamma
jq}(\mu_{jh},\mu_{h})+2A_{\gamma jg}(\mu_{jh},\mu_{h})\Big{]}$
$\displaystyle\times
H(Q,\mu_{h})\left(\frac{Q^{2}}{\mu_{h}^{2}}\right)^{-(C_{F}+C_{A})A_{\Gamma}(\mu_{h},\mu_{jh})-C_{F}A_{\Gamma}(\mu_{h},\mu_{j\ell})}\,,$
(2.6)

and $\widetilde{\phantom{\,}sj}_{g}$ is replaced by

$\widetilde{\phantom{\,}sj}_{q}=\widetilde{j}_{q}\left(\partial_{b}+\ln\frac{Q\mu_{sh}}{\mu_{jh}^{2}}\right)\widetilde{j}_{g}\left(\partial_{b}+\ln\frac{Q\mu_{sh}}{\mu_{jh}^{2}}\right)\widetilde{j}_{q}\left(\partial_{a}+\ln\frac{Q\mu_{s\ell}}{\mu_{j\ell}^{2}}\right)\widetilde{s}_{qg;q}(\partial_{a},\partial_{b})\,,$
(2.7)

with

$a=\eta_{\ell}=2C_{F}A_{\Gamma}(\mu_{j\ell},\mu_{s\ell}),\qquad
b=\eta_{h}=2(C_{F}+C_{A})A_{\Gamma}(\mu_{jh},\mu_{sh})\,.$ (2.8)

For both channels, the resummation is expressed in terms of [37]

$S_{\Gamma}(\nu,\mu)=-\int_{\alpha_{s}(\nu)}^{\alpha_{s}(\mu)}d\alpha\
\frac{\Gamma_{\text{cusp}}(\alpha)}{\beta(\alpha)}\int_{\alpha_{s}(\nu)}^{\alpha}\frac{d\alpha^{\prime}}{\beta(\alpha^{\prime})}\,,$
(2.9)

and

$\displaystyle A_{\Gamma}(\nu,\mu)$
$\displaystyle=-\int_{\alpha_{s}(\nu)}^{\alpha_{s}(\mu)}d\alpha\
\frac{\Gamma_{\text{cusp}}(\alpha)}{\beta(\alpha)},\quad$ $\displaystyle
A_{\gamma}(\nu,\mu)$
$\displaystyle=-\int_{\alpha_{s}(\nu)}^{\alpha_{s}(\mu)}d\alpha\
\frac{\gamma(\alpha)}{\beta(\alpha)}\,.$ (2.10)

The heavy jet mass distribution near the shoulder has the form
$\frac{d\sigma}{d\rho}=\sum c_{n,m}\alpha_{s}^{n}r\ln^{m}r$. Because the
integrals over $m_{\ell}$ and $m_{h}$ in Eq. (2.1) are UV divergent (as
discussed in [26]), we consider first the third derivative
$\frac{d^{3}\sigma}{d\rho^{3}}$. Taking an extra two derivatives turns the
measurement function
$(r+m_{h}^{2}-m_{\ell}^{2})\theta(r+m_{h}^{2}-m_{\ell}^{2})$ in Eq. (2.1) into
a simpler $\delta(r+m_{h}^{2}-m_{\ell}^{2})$. Using additional boundary
conditions obtained from fixed-order matching, this third derivative can be
integrated back up to the physical spectrum, as we will demonstrate in Section
4.1. The expansion of $\frac{d^{3}\sigma}{d\rho^{3}}$ near the shoulder yields
a series in distributions of the form $\left[\frac{\ln^{n}r}{r}\right]_{+}$,
analogous to how the expansion of $\frac{d\sigma}{d\rho}$ in the dijet limit
yields a series of distributions of the form
$\left[\frac{\ln^{n}\rho}{\rho}\right]_{+}$. The resummed form of
$\frac{d^{3}\sigma}{d\rho^{3}}$ is

$\displaystyle\frac{1}{\sigma_{\text{LO}}}\frac{d^{3}\sigma_{g}}{d\rho^{3}}$
$\displaystyle=\Pi_{g}\int_{0}^{\infty}dm_{\ell}^{2}\int_{0}^{\infty}dm_{h}^{2}\
\widetilde{\phantom{\,}sj}_{g}\cdot\left(\frac{Qm_{\ell}^{2}}{\mu_{s\ell}}\right)^{{{\eta_{\ell}}}}\left(\frac{Qm_{h}^{2}}{\mu_{sh}}\right)^{{\eta_{h}}}\frac{1}{m_{\ell}^{2}m_{h}^{2}}\frac{e^{-\gamma_{E}(\eta_{\ell}+\eta_{h})}}{\Gamma(\eta_{\ell})\Gamma(\eta_{h})}$
$\displaystyle\qquad\times\delta(r+m_{h}^{2}-m_{\ell}^{2})\,,$ (2.11)

for the gluon channel, and similarly for the quark channel.

###  2.2 NNLL resummation

For NNLL resummation, we need the three-loop cusp anomalous dimension, the
three-loop $\beta$-function coefficients, the two-loop regular (non-cusp)
anomalous dimensions for the hard, jet and soft functions, and the one-loop
finite parts of the hard, jet and soft functions. The two-loop gluon and quark
jet functions [38, 39] are the same as for other processes, see e.g. [40, 41],
and we summarize all perturbative ingredients in Appendix B.

The Sudakov shoulder hard function is obtained by matching the QCD vector
current onto the corresponding SCET current operator. The virtual diagrams in
SCET are all scaleless in dimensional regularization, so the hard function can
be extracted from the virtual contribution in QCD. For the case of matching
onto three hard partons in the final state, the QCD results can be found in
[42]. The result is summarized in Appendix B. The hard function is the same up
to crossing as for any $q\bar{q}g\gamma^{\star}$ process, such as direct
photon production [40] or 3-jettiness [41], restricted to the trijet
configuration. We have checked the hard function against these results as
well.

The soft function for the heavy jet mass Sudakov shoulder is defined as a
matrix element of Wilson lines in the trijet configuration. The soft function
involves collecting emissions into 6 separate sextant regions shaped like
wedges of a 6-slice orange. Then the soft gluons in the sextants within the
light or heavy hemispheres are combined separately with appropriate
projections. This leads to a two-argument soft function, which in momentum
space is defined as

$S_{g}(k_{\ell},k_{h})=\frac{1}{NC_{F}}\sum_{X_{s}}\bigg{|}\langle
X_{s}|\mathbf{Tr}(Y^{\dagger}_{n_{2}}Y_{n_{1}}T^{a}Y_{n_{1}}^{\dagger}Y_{n_{3}})|0\rangle\bigg{|}^{2}\
\hat{\mathscr{M}}(k_{\ell},k_{h})\ .$ (2.12)

Here, $Y_{n}$ denotes a path-ordered exponential of gauge fields (a Wilson
line) extending from $x=0$ to $x=\infty$ along the $n^{\mu}$ direction. As
written, all of these Wilson lines are in the fundamental representation of
$SU(3)$. We conventionally choose $n_{1}^{\mu}$ to point along the thrust axis
towards the light hemisphere so that Eq. (2.12) corresponds to a soft function
with a gluon jet in the light hemisphere. The other channels arise from
permuting the Wilson line indices. For heavy jet mass, the trijet measurement
function $\hat{\mathscr{M}}(k_{\ell},k_{h})$ is defined as

$\displaystyle\hat{\mathscr{M}}(k_{\ell},k_{h})$ (2.13) $\displaystyle\
=\delta\left(k_{\ell}-\sum_{k_{i}\in|X_{s}\rangle}\left[\theta_{\mathbf{1}}(k_{i})\left(\frac{2}{3}n_{1}\cdot
k_{i}\right)+\theta_{\mathbf{\bar{2}}}(k_{i})\left(\frac{2}{3}N_{2}\cdot
k_{i}\right)+\theta_{\mathbf{\bar{3}}}(k_{i})\left(\frac{2}{3}N_{3}\cdot
k_{i}\right)\right]\right)$ $\displaystyle\ \
\times\delta\left(k_{h}-\sum_{k_{m}\in|X_{s}\rangle}\left[\theta_{\mathbf{\bar{1}}}(k_{m})\left(\frac{2}{3}\bar{n}_{1}\cdot
k_{m}\right)+\theta_{\mathbf{2}}(k_{m})\left(\frac{2}{3}n_{2}\cdot
k_{m}\right)+\theta_{\mathbf{3}}(k_{m})\left(\frac{2}{3}n_{3}\cdot
k_{m}\right)\right]\right).$

where the sum runs over the momenta of all particles in the state
$|X_{s}\rangle$. The notation $\theta_{\mathbf{1}}(k)$ refers to a projection
onto the sextant wedge centered on the direction $n_{1}^{\mu}$ (see [26]),
i.e. $\theta_{\mathbf{1}}(k)=\theta(n_{2}\cdot k-\bar{n}_{2}\cdot
k)\theta(n_{3}\cdot k-\bar{n}_{3}\cdot k)$, and analogously for
$\theta_{\mathbf{2}}(k)$ and $\theta_{\mathbf{3}}(k)$. Similarly,
$\theta_{\mathbf{\bar{1}}}$ refers to a projection onto the wedge in the
direction $\bar{n}_{1}^{\mu}$, which is opposite to the one in the direction
$n_{1}^{\mu}$, $\theta_{\overline{\mathbf{1}}}(k)=\theta(\bar{n}_{2}\cdot
k-n_{2}\cdot k)\theta(\bar{n}_{3}\cdot k-n_{3}\cdot k)$, and similarly for
$\theta_{\mathbf{\bar{2}}}(k)$ and $\theta_{\mathbf{\bar{3}}}(k)$.

In each region only one component of the momentum $k_{i}$ is retained. Among
the six projections, four of them are $\frac{2}{3}n\cdot k$ with $n^{\mu}$ a
lightlike vector pointing to the center of the sextant. However, for the two
sextants adjacent to the light-hemisphere direction $n_{1}^{\mu}$, the
projections involve two spacelike directions

$N_{2}^{\mu}={\bar{n}}_{3}^{\mu}+\frac{1}{2}(n_{1}^{\mu}-{\bar{n}}_{1}^{\mu}),\quad
N_{3}^{\mu}={\bar{n}}_{2}^{\mu}+\frac{1}{2}(n_{1}^{\mu}-{\bar{n}}_{1}^{\mu})\,.$
(2.14)

Ref. [26] computed the one-loop soft anomalous dimension needed for NLL
resummation. The complete one-loop soft function including the soft constant
is calculated in Appendix A. For NNLL resummation, we also need the 2-loop
regular soft anomalous dimension $\gamma_{s}$. This can be fixed by RG
invariance and is expressible in terms of the hard and jet anomalous
dimensions. The result is also summarized in Appendix A.

##  3 Scale setting

In this Section, we show that by choosing scales in position space, the
Sudakov Landau pole is avoided. We also explore reasons why this should be
true.

###  3.1 Momentum space vs. position space

One can do the integrals over $m_{h}^{2}$ and $m_{\ell}^{2}$ exactly before
setting scales. The relevant integral has the form

$\displaystyle f(r)$
$\displaystyle\equiv\frac{1}{\Gamma(a)}\frac{1}{\Gamma(b)}\int_{0}^{\infty}dx\int_{0}^{\infty}dy\
x^{a-1}y^{b-1}\delta(r+y-x)$
$\displaystyle=\frac{1}{\Gamma(a+b)}\left[r^{a+b-1}\frac{\sin(\pi
a)}{\sin(\pi(a+b))}\theta(r)+(-r)^{a+b-1}\frac{\sin(\pi
b)}{\sin(\pi(a+b))}\theta(-r)\right]\,,$ (3.1)

where we shorten notation by denoting $x\equiv m_{\ell}^{2}\ ,\ y\equiv
m_{h}^{2}$. This expression has poles at $a+b\in{\mathbb{Z}}$. The pole at
$a+b=1$ is the Sudakov Landau pole frustrating the numerical resummation of
the Sudakov shoulder logarithms. It comes from the region where $x,y\gg 1$ and
$x-y=r$.

Although the function $f(r)$ in Eq.(3.1) seems complicated, its Fourier
transform is not. It is easiest to compute the Fourier transform before
integrating over $x$ and $y$, doing the integrals in a different order in the
$r>0$ and $r<0$ region:

$\displaystyle\tilde{f}(z)$
$\displaystyle=\frac{1}{\Gamma(a)}\frac{1}{\Gamma(b)}\left[\mathcal{}\int_{0}^{\infty}dr\
e^{izr}\int_{0}^{\infty}dy\ (r+y)^{a-1}y^{b-1}+\int_{-\infty}^{0}dr\
e^{izr}\int dx\ x^{a-1}(x-r)^{b-1}\right]$
$\displaystyle=\frac{\Gamma(1-a-b)}{\pi}\left[\mathcal{}\sin(\pi
a)\int_{0}^{\infty}dr\ e^{izr}\mathcal{}r^{-1+a+b}+\sin(\pi
b)\int_{-\infty}^{0}dr\ e^{izr}(-r)^{-1+a+b}\right]$
$\displaystyle=(-iz)^{-a-b}\frac{\sin(\pi
a)}{\sin(\pi(a+b))}+(iz)^{-a-b}\frac{\sin(\pi b)}{\sin(\pi(a+b))}$
$\displaystyle=(-iz)^{-a}(iz)^{-b}\,.$ (3.2)

So the Fourier transform is remarkably simple.

The simple form of the Fourier transform belies the complicated non-analytic
behavior of $(-iz)^{-a}(iz)^{-b}$. In taking the inverse Fourier transform,
one must carefully adjust the contour to avoid the singularity at $z=0$. For
example, to take the Fourier transform of a single power $z^{-c}$ we must
adjust the contour above or below the branch point. Then closing the contour
either encircles the branch point or gives zero, depending on the sign of the
argument. The result is

$\int_{-\infty}^{\infty}dz\ e^{-izr}(z\pm i\varepsilon)^{-c}=\pm 2\pi
i\frac{(ir)^{c-1}}{\Gamma(c)}\theta(\pm r)\,.$ (3.3)

In a similar manner, the branch point in the position space function
$(-iz)^{-a}(iz)^{-b}$ is ultimately responsible for the complicated piecewise
form in Eq. (3.1).111 As an aside, it is worth noting that the result of
applying the commonly used one-sided Laplace transform to $f(r)$ is not so
simple. The Laplace transform is defined with a semi-infinite integration
$\int_{0}^{\infty}dre^{-\nu r}$ and so for fixed $\nu$ cannot be convergent
for both positive and negative $r$. The best one can do is flip the sign of
the exponent by hand between the $r>0$ and $r<0$ regions leading to
$\mathcal{L}[f](\nu)=\int_{0}^{\infty}dr\ e^{-\nu r}f(r)+\int_{-\infty}^{0}dr\
e^{\nu r}f(r)=\nu^{-a-b}\frac{\sin(\pi a)+\sin(\pi b)}{\sin(\pi(a+b))}\,.$
(3.4) This form is not simple at all, and moreover still has the Sudakov
Landau pole singularities at $a+b\in\mathbb{Z}$.

The pole at $a+b=1$ comes from the large $y$ region:

$f(r)\approx\int_{0}^{\Lambda}dy\ y^{a+b-2}=\frac{\Lambda^{a+b-1}}{a+b-1}\,.$
(3.5)

This is the region of large jet mass shifts $m_{\ell}^{2}=x\gg 1$ and
$m_{h}^{2}=y\gg 1$ where their difference $r=x-y$ is small. If one sets scales
in momentum space, so $a\sim b\sim\alpha_{s}\ln r$ then the pole at $a+b=1$
corresponds to one value of $r$. Large $x$ and $y$ in momentum space
corresponds to small $|z|$ in position space. For small $|z|$ the integral is
approximately

$f(r)\approx\int\frac{dz}{2\pi}(-iz)^{-a}(iz)^{-b}\,.$ (3.6)

This would indeed be singular for $a+b=1$ at small $|z|$. However, when scales
are set in position space, so $a\sim b\sim\alpha_{s}\ln|z|$, then $a$ and $b$
change over the integration region, in contrast to momentum space scale
setting where $a$ and $b$ are fixed. Moreover, the region where $a+b=1$ does
not correspond to small $z$. In this way, the singular region is simply
avoided if one sets scales in position space.

###  3.2 Scale setting

Above we showed how setting scales in position space should avoid the Sudakov
Landau poles. Now we will explore the differences between momentum and
position space-setting in more detail.

The idea behind resummation is to include systematically terms in the
perturbation expansion of the form $c_{n,m}\alpha^{n}L^{m}$ for some infinite
series of $c_{n,m}$ where $L$ is a relevant logarithm and $\alpha$ the
relevant coupling. It can be helpful to distinguish logarithmic counting
schemes. All these schemes agree to some order in resummed perturbation theory
and differ only by higher-order terms. However, these higher-order terms can
be important, and in this case, are the source of the spurious Sudakov Landau
poles as we will see.

A common scheme counts logs in the log of the cumulant [43]

$R(\alpha,L)=\exp[Lg_{1}(\alpha_{s}L)+g_{2}(\alpha_{s}L)+\alpha_{s}g_{3}(\alpha_{s}L)+\cdots]\,.$
(3.7)

Here $R(\alpha,L)$ is the cumulant of some observable. This is in contrast to
the differential distribution, which is the derivative of $R(\alpha,L)$ with
respect to the observable and which generally has an expansion in terms of
plus functions and other distributions. With this organization, $g_{1}(\alpha
L)$ is the leading-logarithmic series, $g_{2}(\alpha L)$ is the next-to-
leading logarithmic series, and so on.

An alternative set of schemes is based on the perturbative expansion of
anomalous dimensions. Here, the leading cusp anomalous dimension $\Gamma_{0}$
and the one-loop $\beta$ function are included for LL resummation,
$\Gamma_{1}$, the two-loop $\beta$ function, and the one-loop regular
anomalous dimensions $\gamma_{0}$ are included for NLL, and one higher order
for each to go to NNLL, and so on. In addition, the finite parts of the hard,
jet, and soft functions are required to one-loop (two-loop) order to achieve
NNLL (N3LL) resummation. Another commonly used combination is that of NNLL RG
evolution (i.e., two-loop regular anomalous dimensions) with two-loop fixed-
order finite parts, which we refer to as NNLL′ and which typically features
improved perturbative uncertainties over NNLL. Within this set of schemes one
still has the freedom to set scales in distribution space, cumulant space, or
position space. Ultimately, we will be setting scales in position space for
our final results.

In [44] different log-counting schemes are compared for a number of dijet-
observables (we review the case of dijet thrust below) and it is shown that
results in all spaces coincide up to subleading logarithmic terms. However, it
is precisely these subleading terms which are responsible for the spurious
divergences in the Sudakov shoulder case. What we will show is that there is
an exact correspondence between the terms associated with the leading cusp
anomalous dimension and the double logarithmic series $(\alpha_{s}L^{2})^{n}$
only if scales are set in position space. In contrast, when scales are set in
momentum space, problematic subleading terms are produced.

### Dijet thrust

Before discussing the shoulder case, let us review scale setting in thrust in
the dijet region. Let us work in the simplest possible setting, the
“$\Gamma_{0}$ approximation” where we take the leading cusp anomalous
dimension $\Gamma_{0}\neq 0$ but we set $\beta_{j}=0$ and $\gamma_{j}=0$ for
$j\geq 0$ and set $\Gamma_{j}=0$ for $j>0$. We also take $H=1$ and set
$\mu_{h}=Q$ throughout as the hard scale is not of particular interest. This
$\Gamma_{0}$-approximation suffices to show how one obtains different results
depending on in which space scales are set. All choices correctly predict the
result at the fixed-scale double-logarithmic order where
$R\sim\exp(\alpha_{s}L^{2})$, but give different results for higher order
terms. It is exactly this issue of how to treat terms beyond the logarithmic
accuracy one it targeting which we wish to address, since for some observables
the higher order terms can induce spurious singularities. In this section, we
define $\hat{\alpha}=\frac{\alpha_{s}}{4\pi}\Gamma_{0}$.

In the $\Gamma_{0}$-approximation, the thrust cumulant is

$R(\tau)=e^{-2\hat{\alpha}C_{F}\ln^{2}\frac{\mu_{h}}{\mu_{j}}-2\hat{\alpha}C_{F}\ln^{2}\frac{\mu_{s}}{\mu_{j}}}\left[\tilde{j}\left(\partial_{\eta}+\ln\frac{Q\mu_{s}}{\mu_{j}^{2}}\right)\right]^{2}\tilde{s}(\partial_{\eta})\left(\frac{\tau
Q}{\mu_{s}}\right)^{\eta}\frac{e^{-\gamma_{E}\eta}}{\Gamma(1+\eta)}\theta(\tau)\,,$
(3.8)

with $\eta=4C_{F}\hat{\alpha}\ln\frac{\mu_{j}}{\mu_{s}}$. From this, we can
read off the canonical scales $\mu_{s}=Q\tau$ $\mu_{j}=\sqrt{Q\mu_{s}}$. Using
these values and setting $\tilde{j}=\tilde{s}=1$ the resummed cumulant in
momentum space takes the form

$\left[R(\tau)\right]_{\text{mom}}=e^{-\hat{\alpha}C_{F}\ln^{2}\tau}\frac{e^{-\gamma_{E}\eta}}{\Gamma(1+\eta)}\theta(\tau)\,,$
(3.9)

with $\eta=-2C_{F}\hat{\alpha}\ln\tau$. $R(\tau)$ has the form in Eq. (3.7)
with $g_{1}=-\hat{\alpha}C_{F}\ln\tau$ and
$g_{2}=\ln\frac{e^{-\gamma_{E}\eta}}{\Gamma(1+\eta)}$. On the other hand if we
Laplace transform before setting scales, then

${\mathcal{L}}[R](\nu)=\frac{1}{\nu}e^{-2\hat{\alpha}C_{F}\ln^{2}\frac{\mu_{h}}{\mu_{j}}-2\hat{\alpha}C_{F}\ln^{2}\frac{\mu_{s}}{\mu_{j}}}\left(\frac{Qe^{-\gamma_{E}}}{\nu\mu_{s}}\right)^{\eta}\,.$
(3.10)

Here we can read off canonical scales as
$\mu_{s}=\frac{Qe^{-\gamma_{E}}}{\nu}$ and $\mu_{j}=\sqrt{Q\mu_{s}}$. With the
Laplace-space scale choice

${\mathcal{L}}[R](\nu)=\frac{1}{\nu}e^{-\hat{\alpha}C_{F}\ln^{2}(\nu
e^{\gamma_{E}})}\,,$ (3.11)

which is a purely LL expansion. The inverse Laplace transform gives

$\left[R(\tau)\right]_{\text{Lap.}}=e^{-\hat{\alpha}C_{F}\ln^{2}\tau}e^{-\hat{\alpha}C_{F}\partial_{\eta}^{2}}\frac{e^{-\gamma_{E}\eta}}{\Gamma(1+\eta)}\theta(\tau)\,.$
(3.12)

This agrees with the momentum space form up to terms formally of order NNLL
and higher, so the results of setting scales in either space are equivalent.
(In practice, setting scales in Laplace space requires some care because one
has to numerically integrate over the Landau pole in QCD [45, 46].)

In summary, we saw from the example of thrust in the dijet limit that setting
scales in position space allows for exact agreement between the logarithmic
counting and the solution of the RG equations. In momentum space, some
subleading terms are generated even if only $\Gamma_{0}$ is kept as in this
simplified example. It is these subleading terms which are problematic in the
Sudakov shoulder case, as we will see next.

### Heavy jet mass shoulder

Now we return to the Sudakov shoulder of heavy jet mass. We work again in the
$\Gamma_{0}$ approximation with $\mu_{h}=Q$, for which the factorization
formula in Eq. (2.11) for the gluon channel becomes

$\displaystyle\frac{1}{\sigma_{\text{LO}}}\frac{d^{3}\sigma_{g}}{d\rho^{3}}=e^{-2\hat{\alpha}C_{F}\ln^{2}\frac{\mu_{h}}{\mu_{jh}}-2\hat{\alpha}C_{F}\ln^{2}\frac{\mu_{sh}}{\mu_{jh}}-\hat{\alpha}C_{A}\ln^{2}\frac{\mu_{h}}{\mu_{j\ell}}-\hat{\alpha}C_{A}\ln^{2}\frac{\mu_{s\ell}}{\mu_{j\ell}}}\times\frac{e^{-\gamma_{E}(a+b)}}{\Gamma(a+b)}$
(3.13) $\displaystyle\
\times\left[\frac{1}{r}\left(\frac{rQ}{\mu_{s\ell}}\right)^{a}\left(\frac{rQ}{\mu_{sh}}\right)^{b}\frac{\sin(\pi
a)}{\sin(\pi(a+b))}\theta(r)+\frac{1}{(-r)}\left(\frac{-rQ}{\mu_{s\ell}}\right)^{a}\left(\frac{-rQ}{\mu_{sh}}\right)^{b}\frac{\sin(\pi
b)}{\sin(\pi(a+b))}\theta(-r)\right].$

Here, $a=2C_{A}A_{\Gamma}(\mu_{j},\mu_{s\ell})$ and
$b=4C_{F}A_{\Gamma}(\mu_{jh},\mu_{s})$. Canonical scales in momentum space are
then

$\mu_{s\ell}=\mu_{sh}=|r|Q,\quad\mu_{j\ell}^{2}=\mu_{jh}^{2}=Q\mu_{s\ell}=Q\mu_{sh}\,.$
(3.14)

Note that $|r|$ must be used to obtain canonical scales for both the $r<0$ and
$r>0$ regions, and hence the scales are continuous but not smooth functions.
With these choices, the momentum space distribution becomes

$\displaystyle\left(\frac{1}{\sigma_{\mathrm{LO}}}\frac{d^{3}\sigma_{g}}{d\rho^{3}}\right)_{\text{mom}}$
(3.15) $\displaystyle\ \
=\frac{1}{r}e^{-\Gamma_{0}\frac{1}{2}C_{A}\hat{\alpha}\ln^{2}|r|-\Gamma_{0}C_{F}\ln^{2}|r|}\frac{e^{-\gamma_{E}(a+b)}}{\Gamma(a+b)}\left[\frac{\sin(\pi
a)}{\sin(\pi(a+b))}\theta(r)-\frac{\sin(\pi
b)}{\sin(\pi(a+b))}\theta(-r)\right],$

with $a=-C_{A}\hat{\alpha}\Gamma_{0}\ln|r|$ and
$b=-2C_{F}\hat{\alpha}\Gamma_{0}\ln|r|$. This has the form of a NLL
distribution $\exp[Lg_{1}(\hat{\alpha}L)+g_{2}(\hat{\alpha}L)]$. The Sudakov
Landau pole is in the NLL contribution, $g_{2}(\hat{\alpha}L)$.

Now let us look at the distribution in position space. Before setting scales,
it has the form

$\displaystyle\widetilde{\sigma}_{g}(z)$
$\displaystyle=\int_{-\infty}^{\infty}dr\ e^{izr}\
\frac{1}{\sigma_{\mathrm{LO}}}\frac{d^{3}\sigma_{g}}{d\rho^{3}}$ (3.16)
$\displaystyle=e^{-2\hat{\alpha}C_{F}\ln^{2}\frac{\mu_{h}}{\mu_{jh}}-2\hat{\alpha}C_{F}\ln^{2}\frac{\mu_{sh}}{\mu_{jh}}-\hat{\alpha}C_{A}\ln^{2}\frac{\mu_{h}}{\mu_{j\ell}}-\hat{\alpha}C_{A}\ln^{2}\frac{\mu_{s\ell}}{\mu_{j\ell}}}\left(-iz\frac{\mu_{s\ell}e^{\gamma_{E}}}{Q}\right)^{-a}\left(iz\frac{\mu_{sh}e^{\gamma_{E}}}{Q}\right)^{-b}.$

We can read off that the canonical scales are

$\mu_{s\ell}=i\frac{Qe^{-\gamma_{E}}}{z},\quad\mu_{sh}=-i\frac{Qe^{-\gamma_{E}}}{z},\quad\mu_{j\ell}^{2}=Q\mu_{s\ell}\quad\text{and}\quad\mu_{jh}^{2}=Q\mu_{sh}\,.$
(3.17)

With these scale choices,

$\widetilde{\sigma}_{g}(z)=\exp\left[-\frac{1}{2}\hat{\alpha}C_{A}\Gamma_{0}\ln^{2}(-ize^{\gamma_{E}})-\hat{\alpha}C_{F}\Gamma_{0}\ln^{2}(ize^{\gamma_{E}})\right]\,.$
(3.18)

Thus in position space this is exactly an LL expansion with no higher terms.

Figure -359: Comparison of the resummed $\frac{d^{3}\sigma_{g}}{d\rho^{3}}$
distribution in the $\Gamma_{0}$ approximation when scales are set in position
or momentum space for $\hat{\alpha}=\frac{\alpha_{s}}{4\pi}=0.01$. The blue
curve shows the result with canonical scales in momentum space from Eq.
(3.14), the orange curve with canonical complex scales from Eq. (3.17) and the
green curve with real scales from Eq. (3.19). The Sudakov Landau poles are
divergences in blue curve at $a+b=1$ ($r=0.012$) and $a+b=2$ ($r=0.00015$) are
seen as an artifact of momentum-space scale setting.

Although we cannot Fourier transform back to momentum space analytically, we
can do so numerically. A comparison between the inverse Fourier transform of
$\widetilde{\sigma}_{g}(z)$ with $\mu_{\rm Can}(z)$, and
$\left(\frac{d^{3}\sigma_{g}}{d\rho^{3}}\right)_{\text{mom}}$ with $\mu_{\rm
Can}(r)$, is shown in Fig. -359. We see that when canonical scales are chosen
in position space, $\mu_{\rm Can}(z)$, the Sudakov Landau pole is gone.

When specifying the canonical scales, we are always free to make alternate
choices that agree up to ${\cal O}(1)$ factors, since such choices still
minimize large logs. An alternative is to use real scales, such as

$\mu_{s\ell}=\mu_{sh}=\frac{Qe^{-\gamma_{E}}}{|z|},\quad\mu_{j\ell}^{2}=\mu_{jh}^{2}=Q\mu_{s\ell}=Q\mu_{sh}\,.$
(3.19)

Choosing real scales like this in position space will still remove the Sudakov
Landau pole. Real scales are additionally easier to implement since they avoid
analytically continuing $\alpha_{s}(\mu)$ to complex scales, as e.g. done for
the resummation of so-called timelike logarithms in gluon-fusion Higgs
production [47, 48, 49, 50, 51, 52]. Moreover, real scales are also preferable
since they guarantee that the distribution will be real when Fourier
transformed back into momentum space. For generic complex scales the higher
order terms, beyond the logarithmic order one is working, can in general be
complex. With the scales in Eq. (3.19) the position space distribution is

$\widetilde{\sigma}_{g}(z)=\exp\left[-\frac{C_{A}+2C_{F}}{2}\hat{\alpha}\Gamma_{0}\ln^{2}(|z|e^{\gamma_{E}})\right]\left(-i\frac{z}{|z|}\right)^{C_{A}\hat{\alpha}\Gamma_{0}\ln(|z|e^{\gamma_{E}})}\left(i\frac{z}{|z|}\right)^{2C_{F}\hat{\alpha}\Gamma_{0}\ln(|z|e^{\gamma_{E}})}\,.$
(3.20)

The numerical inverse Fourier transform of this function is also shown in Fig.
-359.

When setting scales in position space and transforming back, one might be
concerned that an entirely new expression is generated, with no relation to
the one with scales set in momentum space. First of all, we know that if the
computations and matching are done to common fixed order $\alpha_{s}^{n}$,
then the momentum and position space expressions must both be scale-
independent to that order. Thus when we Fourier transform the position space
expression back to momentum space it must agree to that order with the
momentum space expression. Although it is harder to show directly, the series
of large logarithms must also agree. In the $\Gamma_{0}$ approximation, we can
derive an analytic relation between the resummed expressions that also gives
insight into how the Sudakov Landau pole is removed. To Fourier transform the
distribution using canonical position space scales, we can write formally that

$\displaystyle\left(\frac{1}{\sigma_{\mathrm{LO}}}\frac{d^{3}\sigma_{g}}{d\rho^{3}}\right)_{\text{pos}}=\int_{-\infty}^{\infty}\frac{dz}{2\pi}e^{-izr}\tilde{\sigma}_{g}(z)$
$\displaystyle=e^{-\hat{\alpha}\Gamma_{0}\frac{1}{2}C_{A}\partial_{a}^{2}-\hat{\alpha}\Gamma_{0}C_{F}\partial_{b}^{2}}r^{a+b-1}\frac{e^{-\gamma_{E}(a+b)}}{\Gamma(a+b)}\left[\frac{\sin(\pi
a)}{\sin(\pi(a+b))}\theta(r)-\frac{\sin(\pi
b)}{\sin(\pi(a+b))}\theta(-r)\right]_{a=b=0}$
$\displaystyle=e^{-\hat{\alpha}\Gamma_{0}\frac{1}{2}C_{A}\partial_{a}^{2}-\hat{\alpha}\Gamma_{0}C_{F}\partial_{b}^{2}}\left(\frac{1}{\sigma_{\mathrm{LO}}}\frac{d^{3}\sigma_{g}}{d\rho^{3}}\right)_{\text{mom}}\,.$
(3.21)

To get to the last line, we used $f(\partial_{a})r^{a}=r^{a}f(\partial_{a}+\ln
r)$ and then expanded $(\partial_{a}+L)^{2}$ and used
$e^{x\partial_{a}}f(a)=f(a+x)$. In this final form, we must set
$a=-C_{A}\hat{\alpha}\Gamma_{0}\ln|r|$ and
$b=-2C_{F}\hat{\alpha}\Gamma_{0}\ln|r|$ after taking the derivatives as in Eq.
(3). Since $a\sim b\sim\hat{\alpha}\ln r$ each application of
$\hat{\alpha}\partial_{a}^{2}$ or $\hat{\alpha}\partial_{b}^{2}$ gives a
contribution suppressed by $\hat{\alpha}$ which makes the position and
momentum space distributions differ at the NNLL level and beyond.
Nevertheless, these derivatives cure the Sudakov Landau pole: near the pole
$\partial_{a}^{2}\sim\partial_{r}^{2}\to\infty$ and the divergence gets
exponentially suppressed.

###  3.3 Comparison to $q_{T}$ resummation

It has been known for a long time that color-singlet transverse-momentum
($q_{T}$) spectra in hadron collisions feature poles at finite values of
$q_{T}$ when strictly truncating the logarithmic expansion in momentum space
[31], or equivalently when naively setting scales [36] in momentum space.222
We refer the reader to ref. [36] for a detailed discussion of the issues of
naive momentum-space scale setting, as well as for a formal momentum-space
construction that avoids them, essentially by employing cumulant-space scales
for every individual boundary term in the factorization theorem. See also ref.
[35], which avoids the issues of momentum-space scale setting by approximating
logarithms of $q_{T}$ by logarithms of the hardest emission’s $p_{T}$, and
ref. [33, 34, 53, 54], where only a subset of scales is set in momentum space.
In both $q_{T}$ and the Sudakov shoulder, these spurious poles are due to
cancellations between energetic emissions that contribute to the observable:
through momenta with opposite direction in the transverse plane for
$q_{T}=\lvert\vec{q}_{T}\rvert$, and through mass shifts with opposite sign
for the Sudakov shoulder, as in Eq. (1.6). For the $q_{T}$ distribution, it is
in fact standard practice to perform the resummation by setting scales (or
truncating the logarithmic expansion) as a function of the Fourier-conjugate
variable $b_{T}=\lvert\vec{b}_{T}\rvert$ of $q_{T}$ [45, 46, 55], and then
performing a numerical [56] or approximate analytic [57] inverse Fourier
integral. This procedure is well known to lift the spurious poles by retaining
an appropriate set of terms that would be subleading when counting logarithms
directly in momentum space, in close analogy to our discussion in Section 3.
It is thus interesting to ask in what ways our results here differ from the
$q_{T}$ case. Two intriguing differences are

  * •

In the $\vec{q}_{T}$ case, each individual beam and soft function (or TMD PDF)
is an azimuthally symmetric function of the magnitude $k_{T}=\vec{k}_{T}$ of
the total transverse momentum of soft or collinear radiation. Accordingly, the
resulting spectrum is azimuthally symmetric and only a function of
$\lvert\vec{q}_{T}\rvert$. This is distinct from the case of the heavy jet
mass shoulder, where individual jet functions always contribute with a given
fixed sign and the soft function itself is not symmetric under
$m_{h}\leftrightarrow m_{\ell}$. In position space, this physical asymmetry
between hemispheres is encoded in an imaginary part of the Fourier transform
of the cross section, cf. Fig. -356, which is an interesting novel feature of
the Sudakov shoulder case.

  * •

To our knowledge, the heavy-jet mass distribution in the shoulder region is
the first invariant mass-like observable which features these spurious poles,
and for which position-space resummation is required to arrive at physical
results. Examples of transverse momentum-like observables exist whose
resummation can be carried out in momentum space, either because they are
strictly additive scalar quantities (transverse energy $E_{T}$) or because
they involve an ordering that vetoes harder emissions (leading jet $p_{T}$).
Our analysis shows that an observable being transverse is neither sufficient
_nor_ necessary for the appearance of spurious poles, and that position-space
resummation is a crucial tool beyond the case of $q_{T}$. In other language
[58], we can say that this is the first SCETI observable where position space
scale-setting is critical, while earlier examples are in SCETII.

###  3.4 Complex scales and non-global logarithms

In the $\Gamma_{0}$ approximation, the HJM shoulder distribution receives
distinct contributions from the light and heavy hemispheres. These separate
contributions are well-defined since the leading Sudakov double logarithms
come from the soft-collinear limit where the soft emissions are always near a
jet and therefore uniquely associated with a hemisphere. This means we can
choose distinct complex soft scales as in Eq. (3.17).

Beyond the leading-log approximation, when single logs come into play, it is
no longer possible to cleanly separate left and right hemisphere
contributions. In the dijet limit, the hemisphere soft function factorizes as
[30, 59, 13, 14, 60, 61]

$\tilde{s}_{\text{hemi}}(\nu_{1},\nu_{2},\mu)=\tilde{s}_{\mu}(\nu_{1},\mu),\tilde{s}_{\mu}(\nu_{2},\mu)\,\tilde{s}_{f}(\nu_{1},\nu_{2})\,,$
(3.22)

where $\nu_{1}$ and $\nu_{2}$ are the Laplace conjugates of the left and right
hemisphere masses. This factorization is non-trivial. It follows from RG
invariance and symmetry of the hemisphere soft function under the interchange
of the two hemisphere masses. It is also not uniquely defined – one can
shuffle constants among the three factors. We conventionally define
$\tilde{s}_{\mu}(\nu,\mu)$ to be all the terms determined by RG invariance.
The factor $\tilde{s}_{f}(\nu_{1},\nu_{2})$ describes the non-global structure
of the hemisphere soft function. This function is known to 2-loops [60, 61]
and contains the leading non-global logarithm [62].

The position-space trijet hemisphere soft function $S(z_{\ell},z_{h},\mu)$
relevant for the heavy jet mass Sudakov shoulder resummation is superficially
similar. We might like to write

$\tilde{s}(z_{\ell},z_{h},\mu)=\tilde{s}_{\ell}(z_{\ell},\mu)\,,\tilde{s}_{h}(z_{h},\mu)\,\tilde{s}_{\text{ng}}(z_{\ell},z_{h})\,.$
(3.23)

Unfortunately, although we can still use RG invariance to determine its
$\mu$-dependence we no longer have a symmetry in the interchange of the two
hemispheres to give a precise meaning to two factors
$\tilde{s}_{\ell}(z_{\ell},\mu)$ and $\tilde{s}_{h}(z_{h},\mu)$ (beyond the
cusp anomalous dimension). To be explicit, in the gluon channel, the
coefficient of the overall soft anomalous dimension at ${\cal O}(\alpha_{s})$
can be written as $\gamma_{0}^{s}=\gamma_{0}^{sg}+\gamma_{0}^{sqq}$ where

$\gamma_{0}^{sg}=-2C_{A}\ln 3+4C_{F}\ln 2,\quad\gamma_{0}^{sqq}=-4C_{F}\ln
6\,.$ (3.24)

This separation is motivated by the contributions to the soft function from
radiation into the light hemisphere containing the single gluon jet, or the
heavy hemisphere containing the quark and antiquark jets. This suggests that
we might identify $\gamma_{0}^{\ell}=\gamma_{0}^{sg}$ as the anomalous
dimension governing $\tilde{s}_{\ell}(z_{\ell},\mu)$ and
$\gamma_{0}^{h}=\gamma_{0}^{sqq}$ as the anomalous dimension governing
$\tilde{s}_{h}(z_{h},\mu)$. Also, there is no guarantee that resummation
achieved by assuming separate anomalous dimensions will predict any correct
logarithms at higher orders.

Figure -358: Problem with complex scales, using NNLL resummation and no
matching. This figure shows the effect on the resummed distribution if the
two-loop soft anomalous dimension is artificially separated into light and
heavy hemisphere components so that separate complex scales can be chosen.
Different choices of the split introduce unphysical logarithms which distort
the distribution. We see that complex scales introduce spurious uncertainty.

If we cannot define single-scale objects $\tilde{s}_{\ell}(z,\mu)$ and
$\tilde{s}_{h}(z,\mu)$ there is no sense to choosing different scales for the
heavy and light hemispheres. Moreover, it can be problematic to choose
different complex scales for the two hemispheres as in Eq. (3.17). In Fig.
-358 we show the prediction (using the full NNLL setup we describe below) with
canonical complex scales for different choices of how the two-loop soft
anomalous dimension splits into heavy and light components:

$\gamma_{1}^{\ell}=\lambda\gamma_{1}^{s},\quad\gamma_{1}^{h}=(1-\lambda)\gamma_{1}^{s}\,.$
(3.25)

Curves for $\lambda=-\frac{1}{2}\ldots 2$ are shown in Fig. -358 and compared
with the real scale choice in Eq. (3.19). We see that artificially separating
the anomalous dimension can generate large distortions in the spectrum. It
effectively induces large logarithms of the ratio of these scales
$\ln\frac{\mu_{\ell}}{\mu_{h}}\sim\pi$. These logarithms are formally of the
same order as non-global logarithms of $k_{\ell}/k_{h}$ in the
$\tilde{s}_{\text{ng}}$ part of the trijet soft function in Eq. (3.23), which
we necessarily treat at fixed order when evaluating $\tilde{s}_{\text{ng}}$ at
any given individual scale. We thus conclude that there is no sensible way to
split the soft anomalous dimension into light and heavy components, and that a
common scale should be used for all parts of the soft function.

###  3.5 Summary

In this section, we examined scale setting in position, momentum, and Laplace
space for the Sudakov shoulder region of heavy jet mass in the $\Gamma_{0}$
approximation. From the momentum-space point of view, the regions with $r<0$
and $r>0$ are distinct, and the canonical scale choice $\mu_{s}\sim|r|Q$ is
not smooth across $r=0$. However, $r=0$ is not special from the point of QCD:
the shoulder at $r=0$ is an artifact of perturbation theory. In position
space, the natural complex scales are $\mu_{s\ell}\sim-\frac{i}{z}Q$ and
$\mu_{sh}\sim\frac{i}{z}Q$. Choosing scales in position space makes the
distribution smooth across $r=0$ and removes the spurious Sudakov Landau
singularity caused by momentum-space scale setting.

To choose canonical complex scales requires a separation into contributions
from heavy or light hemispheres, which is a sensible criterion only in the LL
approximation. Beyond that, artificially separating scales in the heavy and
light hemispheres is problematic. Instead, one should choose equal soft scales
such as $\mu_{s\ell}=\mu_{sh}\sim\frac{Q}{|z|}$. The Sudakov Landau pole is
still removed as long as scales are chosen in position rather than momentum
space.

##  4 Integration and matching

With the factorization formula and the insights from Section 2 and 3 in hand,
we can now turn to full NNLL resummation. As discussed, choosing scales in
momentum space is problematic: it generates spurious Sudakov Landau poles from
the region where the heavy and light jet mass shifts are both large, but there
difference is still small.

In position space this region corresponds to small $z$. When we set scales in
position space the distribution is heuristically $z^{-\hat{\alpha}\ln
z}=e^{-\hat{\alpha}\ln^{2}z}$ so the small $z$ region becomes Sudakov
suppressed. However, despite this suppression, one must still be careful when
performing the inverse Fourier transform back to momentum space in actual QCD
because of the Landau pole in $\alpha_{s}(\mu)$. In addition, the
factorization formula we have been using so far holds for
$\frac{d^{3}\sigma}{d\rho^{3}}$, so one must integrate over $\rho$ twice to
get back the physical spectrum $\frac{d\sigma}{d\rho}$. This integration
requires boundary conditions generating a linear function $c_{0}+c_{1}\rho$.
These two constants $c_{0}$ and $c_{1}$ correspond to finite remainders after
adding counterterms to remove the logarithmic and linear UV divergences in the
original $\frac{d\sigma}{d\rho}$ distribution in Eq. (2.1). More physically,
they correspond to the fact that the offset and slope of the distribution at
$r=0$ are not predicted by the shoulder factorization formula alone: there can
be contributions from four-parton and higher-order configurations which have
$\rho\sim\frac{1}{3}$ but do not arise from configurations close to the trijet
configuration. The integration constants $c_{0}$ and $c_{1}$ must be
determined by some boundary conditions for the double integral from
$\frac{d^{3}\sigma}{d\rho^{3}}\to\frac{d\sigma}{d\rho}$ as we discuss in this
section.

###  4.1 Integration

In terms of the position space distribution, the factorization formula has the
form

$\frac{1}{\sigma_{\text{LO}}}\frac{d^{3}\sigma_{g}}{d\rho^{3}}=\int_{-\infty}^{\infty}\frac{dz}{2\pi}\,e^{-izr}\,\widetilde{\sigma}_{g}(z)\,,$
(4.1)

where $r=\frac{1}{3}-\rho$ and

$\widetilde{\sigma}_{g}(z)=\Pi_{g}\
\cdot\widetilde{\phantom{\,}sj}_{g}\cdot\left(-iz\frac{e^{\gamma_{E}}\mu_{s\ell}}{Q}\right)^{-a}\left(iz\frac{e^{\gamma_{E}}\mu_{sh}}{Q}\right)^{-b}\,,$
(4.2)

with $\Pi_{g}$ the same as in Eq. (2.3), but now
$\widetilde{\phantom{\,}sj}_{g}$ is a function rather than an operator

$\widetilde{\phantom{\,}sj}_{g}=\widetilde{j}_{q}\left(\ln\frac{Q^{2}e^{-\gamma_{E}}}{iz\mu_{jh}^{2}}\right)\widetilde{j}_{\bar{q}}\left(\ln\frac{Q^{2}e^{-\gamma_{E}}}{iz\mu_{jh}^{2}}\right)\widetilde{j}_{g}\left(\ln\frac{Q^{2}e^{-\gamma_{E}}}{-iz\mu_{jl}^{2}}\right)\widetilde{s}_{q\bar{q};g}\left(\ln\frac{Qe^{-\gamma_{E}}}{-iz\mu_{s\ell}},\ln\frac{Qe^{-\gamma_{E}}}{iz\mu_{sh}}\right).$
(4.3)

For NLL′ or NNLL resummation we expand $\widetilde{\phantom{\,}sj}_{g}$ to
order $\alpha_{s}$. The quark channel is analogous.

To obtain the physical spectrum, we need to inverse Fourier transform and
integrate over $r$ twice. A sufficient condition for $\sigma(\rho)$ to be real
is that $\widetilde{\sigma}^{\star}(-z)=\widetilde{\sigma}(z)$. This condition
on $\widetilde{\sigma}$ holds automatically for real scale choices, since
$(i(-z))^{\star}=iz$. We can thus write the inverse Fourier integral as

$\frac{1}{\sigma_{\text{LO}}}\frac{d^{3}\sigma_{i}}{d\rho^{3}}=\int_{-\infty}^{+\infty}\frac{dz}{2\pi}\,e^{-izr}\,\widetilde{\sigma_{i}}(z)=2\
\text{Re}\left[\int_{0}^{\infty}\frac{dz}{2\pi}e^{-izr}\widetilde{\sigma_{i}}(z)\right]\,.$
(4.4)

Integrating twice with respect to $r$ generates a factor of $\frac{1}{z^{2}}$
in the integrand along with two integration constants. We can then write the
integral as333In practice, we also introduce $u\equiv\frac{1}{1+z}$ to map the
Fourier integration region into a finite region $u\in[0,1]$. And to improve
the convergence, we set a cutoff $\delta=10^{-12}$ in the upper limit of $u$,
and by varying $\delta$ from $10^{-8}$ to $10^{-20}$, we find that the error
is very tiny.

$\displaystyle\frac{1}{\sigma_{\text{LO}}}\frac{d\sigma_{i}}{d\rho}=2\
\text{Re}\left[\int_{0}^{\infty}\frac{dz}{2\pi}{K(z,r)}\widetilde{\sigma_{i}}(z)\right]\,,$
(4.5)

with a kernel function

$K(z,r)=-\frac{e^{-irz}}{z^{2}}+K_{0}(z)+K_{1}(z)\,r\,,$ (4.6)

where $K_{0}$ and $K_{1}$ serve to fix the constants in the boundary
conditions. For example, if we performed the integrals from a lower cutoff
$r_{L}$ so that

$\frac{1}{\sigma_{\text{LO}}}\frac{d\sigma_{i}}{d\rho}=\int_{r_{L}}^{r}dr^{\prime}\int_{r_{L}}^{r^{\prime}}dr^{\prime\prime}\frac{1}{\sigma_{\text{LO}}}\frac{d^{3}\sigma}{d\rho^{3}}\,,$
(4.7)

then we would find

$K^{r_{L}}(z,r)=\frac{e^{-ir_{L}z}-e^{-irz}}{z^{2}}+i\frac{e^{-ir_{L}z}(r_{L}-r)}{z}\,.$
(4.8)

Unlike for dijet observables, where the limit of integration $\rho=0$ or
$\tau=0$ is natural, there is no natural choice for $r_{L}$ for the Sudakov
shoulder. Thus, this boundary condition forcing the distribution to vanish at
$r=r_{L}$ is not ideal.

A more natural boundary condition is to ask that the unmatched distribution
when expanded to fixed order reproduces the leading behavior of the full
theory distribution near $r=0$. It does this automatically for all the
$r\ln^{m}r$ terms, which are true predictions of the effective theory
factorization formula. It also reproduces the difference in slopes in the left
and right shoulders. This difference comes from integrating over the
$\delta$-function terms in the triple-derivative. For example,

$\int_{r_{L}}^{r}dr^{\prime}\int_{r_{L}}^{r^{\prime}}dr^{\prime\prime}\delta(r^{\prime\prime})=r\theta(r)-r\theta(r_{L})=\frac{1}{2}r\Big{[}\theta(r)-\theta(-r)\Big{]}+\left[\frac{1}{2}-\theta(r_{L})\right]r\,.$
(4.9)

The first term on the right is independent of $r_{L}$ and therefore a
prediction of the theory. The value for the second term is determined by one
of the boundary conditions.

More generally, consider the unmatched resummed distribution, expanded in
$\alpha_{s}$. The distribution
$g(r)=\frac{1}{\sigma_{0}}\frac{d\sigma}{d\rho}$ has the generic form

$\displaystyle\hskip 28.45274ptg^{\text{res}}(r)$
$\displaystyle=\theta(r)\sum_{m\geq 1}d_{m}^{+}r\ln^{m}r+\theta(r)\sum_{m\geq
1}d_{m}^{-}r\ln^{m}(-r)$ $\displaystyle\ \
+br[\theta(r)-\theta(-r)]+c_{0}+c_{1}\,r\,,\hskip 28.45274pt$ (4.10)

with $c_{1},b,d_{m}^{\pm}$ all functions of $\alpha_{s}$, but only $c_{0}$ and
$c_{1}$ depend on the choice of boundary conditions. That is, $d^{\pm}$ and
$b$ are predicted by the factorized expression while $c_{0}$ and $c_{1}$ are
not. We can then demand that the offset at $r=0$ and the sum of the slopes on
the left and right shoulders at $r=0$ agree with a fixed-order calculation to
fix $c_{0}$ and $c_{1}$. Moreover, using perturbation theory to determine the
boundary conditions implies that the uncertainty on these constants will be
included in the uncertainty associated with the fixed-order contribution. Thus
we will not have to account for it separately.

At leading order, the distribution in the two mass shifts is

$\frac{1}{\sigma_{0}}\frac{d^{2}\sigma}{dm_{\ell}^{2}dm_{h}^{2}}=\delta(m_{\ell}^{2})\,\delta(m_{h}^{2})\,,$
(4.11)

so the integrals over $m_{\ell}^{2}$ and $m_{h}^{2}$ in Eq. (1.6) are finite
and the unmatched distribution is

$g(r)=\int_{0}^{\infty}dm_{h}^{2}\int_{0}^{\infty}dm_{\ell}^{2}\,\frac{d^{2}\sigma}{dm_{\ell}^{2}dm_{h}^{2}}(r+m_{h}^{2}-m_{\ell}^{2})\,\theta(r+m_{h}^{2}-m_{\ell}^{2})=r\theta(r)\,.$
(4.12)

We read off $d_{m}^{\pm}=0$, $b=c_{1}=\frac{1}{2}$ and $c_{0}=0$. Note that at
leading order there is no UV divergence in the integrals over $m_{\ell}^{2}$
and $m_{h}^{2}$ so the integration constants are calculable without using the
full theory. Since $g^{\prime\prime}(r)=\delta(r)$ we can implement these
integration constants with an integral kernel

$K(z,r)=\frac{1}{z^{2}}\left[1-e^{-izr}+r(1-e^{iz})\right]\,.$ (4.13)

This kernel when integrated against $\tilde{\sigma}(z)=1$ gives $r\theta(r)$
as desired.

At higher orders we need to modify the boundary conditions by the inclusion of
$\alpha_{s}$ corrections from fixed-order perturbation theory. Furthermore, we
also want to add any fixed order contributions that are not enhanced by large
logarithms, and hence not part of the resummed cross-section (so called “non-
singular” terms). This can be accomplished by writing

$\displaystyle\frac{d\sigma^{\text{match,sh}}}{d\rho}$
$\displaystyle=\frac{d\sigma^{\text{FO}}(\mu_{\text{FO}})}{d\rho}+2\text{Re}\biggl{\\{}\int_{0}^{\infty}\frac{dz}{2\pi}K(z,r)\left[\tilde{\sigma}^{\text{sh}}(z,\mu_{\text{sh}})-\widetilde{\sigma}^{\text{sh}}(z,\mu^{\text{FO}})\right]\biggr{\\}}$
$\displaystyle=\overbrace{2\text{Re}\biggl{\\{}\int_{0}^{\infty}\frac{dz}{2\pi}K(z,r)\,\tilde{\sigma}^{\text{sh}}(z,\mu_{\text{sh}})\biggr{\\}}}^{\text{resummed}}+\underbrace{\biggl{[}\frac{d\sigma^{\text{FO}}(\mu_{\text{FO}})}{d\rho}-\overbrace{\frac{d\sigma^{\text{sh}}(\mu_{\text{FO}})}{d\rho}}^{\text{singular}}\biggr{]}}_{\text{non-
singular}}\,,$ (4.14)

where $\mu^{\text{res}}$ stands for the scale choices leading to resummation
and $\mu^{\text{FO}}$ refers to setting all the scales equal to a common
fixed-order scale. When writing the matched cross section as on the first
line, the integration constants will contribute only beyond the order to which
the resummed cross section is matched, since the latter is scale-invariant to
that order. Thus the uncertainty on the constants will be within the fixed-
order (hard scale) uncertainty. We therefore use the leading-order kernel in
Eq. (4.13) for simplicity. On the second line, we have rearranged terms and
performed the analytic $z$ integral over the fixed-order singular cross
section to isolate the “non-singular” term in square brackets, which is
analytic near the shoulder.

###  4.2 Profile functions

Although choosing canonical real scales resums the correct logarithms to NNLL
order, in practice it is important to choose scales which avoid the Landau
pole at $\mu=\Lambda_{\rm QCD}$, and which allow for resummation to turn off
in regions that are more accurately described by an alternate description.
These features can be achieved with profile functions [63, 19]. In our case,
we want to turn off the resummation for the shoulder region as we approach the
dijet region, where conventional dijet resummation applies, and to switch to
fixed-order perturbation theory when we are far enough above the shoulder
threshold.

To construct profile functions, we follow the procedure used for $q_{T}$
resummation in [64]. We start with the canonical (real) scales, as justified
in Section 2. To ensure that the soft and jet scales do not hit the Landau
pole of QCD, we adjust these canonical scales to be

$\displaystyle\mu_{h}^{\text{can}}=Q\,,\qquad\mu_{s}^{\text{can}}(z)=\sqrt{\Big{(}\frac{Qe^{-\gamma_{E}}}{|z|}\Big{)}^{2}+\Big{(}\mu_{s}^{\text{min}}\Big{)}^{2}}\,,\qquad\mu_{j}^{\text{can}}(z)=\sqrt{\mu_{h}^{\text{can}}\mu_{s}^{\text{can}}}\,,$
(4.15)

and take $\mu_{s}^{\text{min}}=2$ GeV. We would like the canonical scales to
be used for a certain range of $\rho$ and then smoothly turn off away from the
resummation region. To do so, we choose central jet and soft profile scales
based on the canonical scales as

$\displaystyle\mu_{j,s}^{\text{pro}}(z,\rho)=\mu_{h}^{1-g(\rho)}[\mu^{\mathrm{can}}_{j,s}(z)]^{g(\rho)}\,,$
(4.16)

where $0\leq g(\rho)\leq 1$ in the exponent determines whether the resummation
is fully on ($g=1$) or turned off to recover the fixed-order singular result
($g=0$), and smoothly interpolates between the two cases through a weighted
geometric mean. Since the RG evolution resums logarithms of $\mu_{s}/\mu_{j}$
and $\mu_{j}/\mu_{h}$, this effectively amounts to having $g(\rho)$ multiply
logarithms of $z$ in the resummation exponent.

To achieve a matched result, the shoulder resummation needs to turn off
outside of the resummation region. To implement this, we take $g(\rho)$ to be
a product of profile functions governing the turn-off towards the left and the
right of the shoulder,

$\displaystyle g(\rho)=g_{L}(\rho)\,g_{R}(\rho)\,.$ (4.17)

We use two different functional forms for the individual profile functions,
each with their own advantages. For the left profile, we use a sigmoid
function:

$\displaystyle g_{L}(\rho)=1-\frac{1}{1+e^{(\rho-\rho_{L})/\sigma_{L}}}\ =\
\raisebox{-0.5pt}{\includegraphics[width=113.81102pt]{sigmoid.pdf}}\,.$ (4.18)

The advantage of using an infinitely differentiable function like the sigmoid
is that it does not introduce unphysical kinks in the matched distribution
which may bias comparisons to data. A disadvantage is that the sigmoid
function turns the resummation exactly on or off only in the asymptotic limits
$\rho\to\pm\infty$. (In practice, $g_{L}(\rho)\geq 0.9$ for
$\rho>\rho_{L}+2\sigma_{L}$, which leads to scales that are equivalent to
canonical ones within the variations we consider. Similarly, we find that
$g_{L}(\rho)\leq 0.1$ for $\rho<\rho_{L}-2\sigma_{L}$ is close enough to the
strict fixed-order limit to the left of the shoulder for our purposes.) For
the right profile, we use a piecewise quadratic function [64]:

$\displaystyle g_{R}(\rho)=\begin{cases}1\,,\quad&\rho<\rho_{R1}\,,\\\
1-\frac{2(\rho-\rho_{R1})^{2}}{(\rho_{R1}-\rho_{R2})^{2}}\,,\quad&\rho_{R1}<\rho<\tfrac{\rho_{R1}+\rho_{R2}}{2}\,,\\\
\frac{2(\rho-\rho_{R2})^{2}}{(\rho_{R1}-\rho_{R2})^{2}}\,,\quad&\tfrac{\rho_{R1}+\rho_{R2}}{2}<\rho<\rho_{R2}\,,\\\
0\,,\quad&\rho>\rho_{R2}\,,\end{cases}\quad=\
\raisebox{-0.5pt}{\includegraphics[width=113.81102pt]{quad.pdf}}\,.$ (4.19)

This form has the advantage that the resummation is exactly off at
$\rho\geq\rho_{R2}$. Since the cross section on the right shoulder falls
rapidly, one must turn off resummation rather quickly to avoid unphysical
negative values from the resummed distribution [19]. This is easiest to
achieve with a quadratic profile.

Figure -357: The left panel shows the weight function used in the shoulder
profiles to interpolate between canonical scales ($g=1$) and fixed scales
($g=0$). The dashed line indicates the central transition point for the left-
shoulder profile, at $\rho_{L}=0.2$. The right panel shows the position space
scales we use for canonical scales ($g\approx 1$; red dot on left, red curves
on right) and for scales which are part way on the transition to fixed scales
($g=0.3$; green dot on left, green curves on right). We take
$Q=100~{}\text{GeV}$ and $\mu_{s}^{\min}=2~{}\text{GeV}$,
$\rho_{L}=0.2,\sigma_{L}=0.025,\rho_{R1}=\frac{1}{3}$ and $\rho_{R2}=0.342$.
The region where $|z|<e^{-\gamma_{E}}=0.56$ has $\mu_{j},\mu_{s}>Q$ but gives
a negligible contribution to the momentum space distribution.

The left panel of Fig. -357 shows the combined interpolating function
$g(\rho)$ we use for the full shoulder region. The right panel of this figure
shows the scales as a function of $z$ for a point in the canonical region
where $g=1$, and a point with $g=0.3$ in the transition region, illustrating
the weighted geometric mean of scales. Although the region with $|z|\lesssim
0.56$ where the jet and soft scales are larger than $Q$ is not excluded from
the Fourier integral, this region contributes negligibly.444In the
$|z|\lesssim 0.56$ region, since $r$ is small, we can approximate $e^{izr}=1$
in the inverse Fourier transform. Thus not only does the region $|z|\lesssim
0.56$ not contribute much overall (because the integrand is bounded), but the
amount it contributes is also independent of $r$ and therefore its leading
contribution is fixed by the $c_{0}+c_{1}r$ boundary condition. We can
understand this from looking directly at $\tilde{\sigma}(z)$ as shown in Fig.
-356.

Figure -356: The Fourier transform $\widetilde{\sigma}(z)$ of the cross
section with NNLL resummation and real scale choices. Right panel zooms in for
the small $z$ region. The region within the two grey dashed lines have
canonical scales above $Q$ but contributes negligibly to the momentum space
distribution.

To determine the profile midpoint $\rho_{L}$ and the transition width
$\sigma_{L}$ on the left shoulder, we compare the size of the singular cross
section predicted by the shoulder factorization at fixed order to the total
cross section in the top left panel of Fig. -355. The non-singular remainder,
which is given by the difference of the two, has the same size as the singular
cross section around $\rho\approx 0.15$, suggesting that the resummation
should be turned off around this point. Taking $\rho\lesssim 0.15$ to be the
region where the resummation is (nearly) off, this suggests
$\rho_{L}-2\sigma_{L}=0.15$. On the other hand, the singular cross section
dominates over the remainder by a factor $\geq 5$ down to $\rho\approx 0.25$,
which motivates keeping the resummation on until $\rho_{L}+2\sigma_{L}\approx
0.25$. Thus the $g_{L}(\rho)$ sigmoid profile is configured to make the bulk
of the transition between $\rho=0.25$ and $\rho=0.15$, as evident in Fig.
-357.

For the right shoulder, a comparison of singular and non-singular
contributions is not a useful tool to determine the resummation region. The
problem is that the cross section dies off very fast for
$\rho\gtrsim\frac{1}{3}$ and moreover there is a second shoulder determined by
the 4-parton threshold at $\rho=0.42$. In fact, the resummed distribution with
canonical scale choices, although smooth at $\rho=\frac{1}{3}$, falls rapidly
for larger $\rho$ and becomes negative by $\rho=0.345$ (see Fig. -354 below).
To obtain a physical, positive cross section, the resummation must be turned
off at least by then, so we take $\rho_{R2}=0.342$. We take
$\rho_{R1}=\frac{1}{3}$ to at least extend the right shoulder resummation as
far as possible. Part of the reason for using a quadratic profile in this
region, which has $g_{R}(\rho)=0$ exactly for $\rho\geq\rho_{R2}$ is that even
the very small difference between the sigmoid and quadratic profiles has big
effects relative to the value of the cross section. When we match the resummed
distribution to fixed order, we add the singular and nonsingular contributions
(see Eq. (4.1)). In the fully fixed order region the singular and nonsingular
terms are individually large with opposite sign, and must sum to give a cross
section which is close to zero. Even a small amount of leftover resummation
from not turning off higher order logarithms can spoil this sensitive
cancellation. The cancellation in the fixed order region is guaranteed by the
quadratic profile in Eq. (4.19), which imposes $g(\rho\geq\rho_{R2})=0$
exactly, as can be seen in Fig. -357. The same delicate cancellations do not
take place in the region controlled by $g_{L}(\rho)$, allowing us to use the
smoother sigmoid for the left shoulder transition.

###  4.3 Dijet matching

Figure -355:  Comparison of fixed order contributions from various sources to
determine which terms dominate in various regions. In the top-left panel we
compare the should singular terms $\sigma_{\text{sh}}^{\text{LO}}$ with the
full LO and their difference. In the top right panel we compare the dijet
singular terms $\sigma_{\text{dij}}^{\text{LO}}$ with the full LO and their
difference. The bottom left shows the negative of the remainder when both
expansions are subtracted from LO. The bottom right panel shows that taking a
linear combination of the two LO singular contributions weighted by the
function $g(\rho)$ used for the profile functions can nearly completely
reproduce the full LO distribution. In our actual matched results, the
transition between different regions is implemented through scale choices, not
a simple linear combination.

To obtain the best possible prediction valid across phase space, we should
match the resummed shoulder distribution to both the fixed-order result and
the resummed dijet distribution. The resummation for heavy jet mass in the
dijet limit has been done up to the N3LL level [14]. The general form for the
resummed distribution in the dijet region is

$\displaystyle\frac{1}{\sigma_{0}}\frac{d\sigma^{\text{dij}}}{d\rho}$
$\displaystyle=\exp\left[4C_{F}S(\mu_{h},\mu_{j})+4C_{F}S(\mu_{s},\mu_{j})-2A_{H}(\mu_{h},\mu_{s})+4A_{J}(\mu_{j},\mu_{s})\right]$
$\displaystyle\quad\times\left(\frac{Q^{2}}{\mu_{h}^{2}}\right)^{-2C_{F}A_{\Gamma}(\mu_{h},\mu_{j})}H^{\text{dij}}(Q^{2},\mu_{h}^{2})\
\tilde{j}_{q}\left(\partial_{\eta_{1}}+\ln\frac{Q\mu_{s}}{\mu_{j}^{2}}\right)\tilde{j}_{q}\left(\partial_{\eta_{2}}+\ln\frac{Q\mu_{s}}{\mu_{j}^{2}}\right)$
$\displaystyle\quad\times\tilde{s}_{\mu}(\partial_{\eta_{1}})\tilde{s}_{\mu}(\partial_{\eta_{2}})\tilde{s}_{f}\left(\partial_{\eta_{1}}-\partial_{\eta_{2}}\right)\frac{1}{\rho}\left(\frac{\rho
Q}{\mu_{s}}\right)^{\eta_{1}+\eta_{2}}\frac{e^{-\gamma_{E}\eta_{1}}}{\Gamma(\eta_{1})}\frac{e^{-\gamma_{E}\eta_{2}}}{\Gamma(\eta_{2})}\,,$
(4.20)

where after taking derivatives one sets

$\displaystyle\eta_{1}=\eta_{2}=2C_{F}A_{\Gamma}(\mu_{j},\mu_{s})\,.$ (4.21)

As discussed in Section 3, for the dijet region it is not necessary to set
scales in position space (although one could). For simplicity, for the
canonical scales we use conventional canonical scales in momentum space,

$\displaystyle\mu_{h}=Q\,,\qquad\mu_{j,{\rm dij}}^{\rm
can}=Q\sqrt{\rho}\,,\qquad\mu_{s,{\rm dij}}^{\rm can}=Q\rho\,.$ (4.22)

Before we present the complete matched result including shoulder resummation,
we first discuss how to match the resummed dijet distribution to fixed order
using standard procedures [19] in a setup where the shoulder is not accounted
for as a special region. Here one adds a non-singular remainder with a common
scale $\mu_{\rm FO}$ to the resummed dijet cross section evaluated with soft
and jet profile scales (denoted together by $\mu^{\text{pro}}_{\rm dij}$)

$\displaystyle\frac{d\sigma^{\text{match,dij}}}{d\rho}=\overbrace{\frac{d\sigma^{\text{dij}}(\mu^{\text{pro}}_{\rm
dij})}{d\rho}}^{\text{resummed}}+\underbrace{\biggl{[}\frac{d\sigma^{\text{FO}}(\mu_{\rm
FO})}{d\rho}-\overbrace{\frac{d\sigma^{\text{dij}}(\mu_{\rm
FO})}{d\rho}}^{\text{singular}}\biggr{]}}_{\text{non-singular}}\,.$ (4.23)

To ensure that the plain fixed-order result is recovered in the tail, it is
common to turn off the dijet resummation with profile functions that
interpolate between canonical scales in Eq. (4.22) and equal scales
$\mu_{s}=\mu_{j}=\mu_{h}$ directly as a function of $\rho$. When matching to
the plain fixed-order result, appropriate parameters for the dijet profile can
be inferred by the top right panel of Fig. -355 in analogy to the discussion
for the shoulder. When the LO cross section vanishes at $\rho=\frac{1}{3}$,
the dijet singular cross section exactly cancels the non-singular. Any small
mismatch between the singular and non-singular, resulting for example from
scales that are close but not equal in the profiles, can lead to large
relative errors on the prediction. This suggests that the dijet resummation
should be turned off before reaching this point. In practice, when comparing
to partially matched results based on Eq. (4.23), we will make use of the
profile functional form developed in ref. [15] for the $C$ parameter
distribution, adjusting the point $t_{s}$ (defined in that reference) where
the resummation is exactly turned off to be $t_{s}=\frac{1}{3}$.

We now turn to the complete matched result incorporating both dijet and
shoulder resummation. The matched cross section takes the following form,

$\displaystyle\frac{d\sigma^{\text{match}}}{d\rho}=\frac{d\sigma^{\text{dij}}(\mu_{\rm
dij}^{\text{pro}})}{d\rho}+\frac{d\sigma^{\text{sh}}(\mu_{\rm
sh}^{\text{pro}})}{d\rho}+\left[\frac{d\sigma^{\text{FO}}(\mu_{\rm
FO})}{d\rho}-\frac{d\sigma^{\text{dij}}(\mu_{\rm
FO})}{d\rho}-\frac{d\sigma^{\text{sh}}(\mu_{\rm FO})}{d\rho}\right]\,.$ (4.24)

This matched distribution has the formal properties that far from the
shoulder, when the shoulder scales are equal to the fixed-order scale, it
reduces to the dijet distribution matched to fixed order as in Eq. (4.23).
Similarly, parametrically far from the dijet region, where the dijet scales
reduce to fixed order scales, the dijet resummation turns off and it reduces
to a matched distribution of shoulder and fixed order only. We note that
matching different regions of a one-dimensional distribution is in a sense
easier than the joint resummation required for double-differential
distributions as in [64], where generically there is a region where two
separate types of logarithms can be simultaneously large.

For Eq. (4.24) to have the right limits, we moderate the dijet resummation
through profile functions similar to Eq. (4.16). We define the profile scales
with a function $g_{\text{dij}}(\rho)$ playing the role that $g(\rho)$ plays
for shoulder resummation:

$\displaystyle\mu_{j,\text{dij}}^{\text{pro}}(\rho)$
$\displaystyle=\mu_{h}^{1-g_{\text{dij}}(\rho)}\big{[}\mu^{\mathrm{can}}_{j,\text{dij}}(\rho)\big{]}^{g_{\text{dij}}(\rho)}\,,$
$\displaystyle\mu_{s,\text{dij}}^{\text{pro}}(\rho)$
$\displaystyle=\mu_{h}^{1-g_{\text{dij}}(\rho)}\big{[}\mu^{\mathrm{can}}_{s,\text{dij}}(\rho)\big{]}^{g_{\text{dij}}(\rho)}\,.$
(4.25)

Note that for, dijet resummation, the profile scales only depend on $\rho$, in
contrast to Eq. (4.16) where they also depend on $z$, because we perform the
resummation directly in momentum space.

The weight function $g_{\text{dij}}(\rho)$ should interpolate between the
region where the dijet resummation is fully on ($g_{\text{dij}}=1$) and the
region where it is off ($g_{\text{dij}}=0$). To determine the parameters for
$g_{\text{dij}}$, we return to the comparison of the singular and non-singular
contributions, but now look at the non-singular with both singular shoulder
and singular dijet subtracted. This overall non-singular contribution is given
by the expression in square brackets in Eq. (4.24). At leading order it is
shown (multiplied by $-1$) in the bottom-left panel of Fig. -355. From this
plot, we conclude there is no reason to adjust the shoulder profiles when
matching to dijet – the shoulder singular still accounts for nearly the
complete LO distribution for $\rho\gtrsim 0.25$ and the shoulder singular is
the same size as the overall non-singular at $\rho\approx 0.15$. However, for
the dijet, when matching to the shoulder resummation it is no longer sensible
to maintain dijet resummation up to $\rho\approx 0.33$. If dijet resummation
is allowed past the point where the dijet singular and the full nonsingular
curves are nearly equal in magnitude, $\rho\approx 0.25$, the resummation
would upset the large cancellations between them, biasing the result. If dijet
resummation is turned off at $\rho\gtrsim 0.25$, these cancellations are
guaranteed.

A profile function for the matched dijet resummation consistent with these
constraints is simply the same profile as for the shoulder, but growing in the
opposite direction. That is, we take

$\displaystyle
g_{\mathrm{dij}}(\rho)=\frac{1}{1+e^{(\rho-\rho_{\mathrm{dij}})/\sigma_{\mathrm{dij}}}}\,g_{R}(\rho)\,,$
(4.26)

with $\rho_{\mathrm{dij}}=\rho_{L}$ and $\sigma_{\mathrm{dij}}=\sigma_{L}$
equal to those of $g_{L}$ in Eq. (4.18). The additional factor of
$g_{R}(\rho)$ in Eq. (4.26) ensures that all resummation (dijet and shoulder)
is turned off in the far tail $\rho\gtrsim 0.342$, maintaining consistency
with our choice of profiles for the right shoulder.

###  4.4 Uncertainties

A common way to assess theoretical uncertainties due to missing higher orders
is to vary the factorization scale. We consider three scale variations in the
shoulder region. A hard variation involves changing $Q\to 2^{v_{h}}Q$ in all
the scales in Eq. (4.15). A soft variation involves varying the soft scale in
the same was as for the hard variation, but not varying the hard and
maintaining $\mu_{j}=\sqrt{\mu_{s}\mu_{h}}$. For the jet scale variation, we
take $\mu_{j}=(\mu_{h}\mu_{s})^{v_{j}}$ and allow $v_{j}$ to differ from
$\frac{1}{2}$. That is, these variations correspond to

$\displaystyle\mu_{h}$ $\displaystyle=2^{v_{h}}Q\,,$
$\displaystyle\mu_{s,\text{sh}}(z)$
$\displaystyle=\sqrt{\Big{(}2^{v_{h}}2^{v_{s}}\frac{Qe^{-\gamma_{E}}}{|z|}\Big{)}^{2}+\Big{(}\mu_{s}^{\text{min}}\Big{)}^{2}}\,,$
$\displaystyle\mu_{j,\text{sh}}(z)$
$\displaystyle=\bigl{[}\mu_{s,\text{sh}(z)}\bigr{]}^{v_{j}}\mu_{h}^{1-v_{j}}\,.$
(4.27)

These varied scales are then inserted into the profiles using Eq. (4.16). To
vary the scales by a factor of 2 we separately vary $v_{h}$ and $v_{s}$ from
$-1$ to $1$. This variation is inherited by the jet scale, maintaining its
canonical relation with the soft and hard jet scale. In addition, we vary the
jet scale away from its canonical value by varying the power
$v_{j}=\\{0.4,0.6\\}$ in the exponent.

For the dijet resummed distribution we consider scale variations of the form:

$\displaystyle\mu_{h}$ $\displaystyle=2^{v_{h}}Q\,,$
$\displaystyle\mu_{s,\text{dij}}(\rho)$
$\displaystyle=2^{v_{h}}2^{v^{\prime}_{s}}Q\rho\,,$
$\displaystyle\mu_{j,\text{dij}}(\rho)$
$\displaystyle=\bigl{[}\mu_{s,\text{dij}(\rho)}\bigr{]}^{v_{j}^{\prime}}\mu_{h}^{1-v_{j}^{\prime}}\,,$
(4.28)

i.e., we vary the dijet scales in a fashion that is similar to the shoulder
scales. We use a single hard scale variation $v_{h}$, which is the same for
both dijet and shoulder logs and appears in the nonsingular cross section in
Eq. (4.24) as well, where we take $\mu_{\mathrm{FO}}=\mu_{h}$.

We independently vary the profile parameters $\rho_{L}$, $\sigma_{L}$,
$\rho_{\mathrm{dij}}$, and $\sigma_{\mathrm{dij}}$. We do not consider
variations of the profile functions in the right shoulder since resummation
must turn off quite quickly on the right shoulder, so any reasonable
variations which do not make the cross section negative are in any case
smaller than the hard scale variation.

A summary of the profile parameters along with their variation ranges is shown
in Table 1. As a preliminary estimate of the overall perturbative uncertainty,
we consider an envelope of all the scale variations and the profile parameter
variations.

Parameter | Central value | Variation Range
---|---|---
$v_{h}$ | $0$ |  $-1$ to $1$
$v_{s},v_{s}^{\prime}$ | $0$ |  $-1$ to $1$
$v_{j},v_{j}^{\prime}$ | $0.5$ |  $0.4$ to $0.6$
$\rho_{L},\rho_{\rm dij}$ | 0.2 |  $0.18$ to $0.22$
$\sigma_{L},\sigma_{\rm dij}$ | $0.025$ |  $0.02$ to $0.03$
Table 1: Parameter ranges used for the uncertainty estimates in both the
shoulder profile function and dijet profile function. We fix the profile
function for the right shoulder to have $\rho_{R1}=\frac{1}{3}$ and
$\rho_{R2}=0.342$, and fix $\mu_{s}^{\text{min}}=2\,\mathrm{GeV}$. Varying
these parameters has negligible effect on the distribution. For the dijet
resummation (when matching to the fixed-order result only), we use the profile
function from [15], with the following set of parameters defined in that
reference: $n_{0}=2$, $n_{1}=10$, $t_{2}=0.2\pm 0.025$, $t_{s}=\frac{1}{3}\pm
0.025$, $r_{s}=1^{\times 2}_{\div 2}$, $e_{j}=0\pm 1$, $e_{h}=1^{\times
2}_{\div 2}$.

##  5 Numerical results

Figure -354:  Plots of Sudakov shoulder resummation in heavy jet mass
distribution, with matching to FO distributions. Profiled scales are indicated
by solid curves, and interpolate between canonical scales near
$\rho=\frac{1}{3}$ (dot-dashed curves) and fixed order scales
($\mu_{h}=\mu_{j,\text{sh}}=\mu_{s,\text{sh}}$, dashed curves). The right
panel indicates a zoomed in view of the distribution near $\rho=\frac{1}{3}$.

With all the ingredients in place, we can now look at numerical predictions.
All our results are given for $e^{+}e^{-}$ hadron events at $Q=91.2$ GeV with
$\alpha_{s}=0.1179$ [65].

We first examine the prediction using Sudakov shoulder resummation without
dijet matching in Fig. -354. We compare canonical scales choices, as given in
Eq. (4.15), to profile scales, as given in Eq. (4.16), and compare to fixed
order distributions at LO and NLO [66, 67]. One can also see from this figure
that the profile functions smoothly transition out of the region where
shoulder logarithms are large ($\rho\sim\frac{1}{3}$) into the fixed order
region at smaller $\rho$. This figure also confirms that setting scales in
position space eliminates the Sudakov Landau poles encountered in Ref. [26].
The right panel in Fig. -354 zooms in on the region around $\rho=\frac{1}{3}$.
It shows that in the matched curves, the non-analytic behavior (kink) of the
fixed-order distributions is completely resolved by the resummation of
shoulder logarithms, leaving behind a distribution with continuous first
derivative. Although the NLL+LO matched curve goes negative for
$\rho>\frac{1}{3}$, this is due to our boundary conditions on integrating the
third derivative which makes it match exactly to LO at $\rho=\frac{1}{3}$. For
the NNLL+NLO curve, the offset is positive and the profile functions of the
right shoulder are chosen to turn off the shoulder resummation quickly enough
to avoid a negative cross section. Comparing the NLO curve to those at
NNLL+NLO, we see that the shoulder resummation has a considerable impact for
values of $\rho\gtrsim 0.25$.

Figure -353: Relative contributions dijet and shoulder resummation to the
matched distribution. Right panel shows the ratio with respect to
$r=\frac{1}{3}-\rho$, i.e.
$\frac{\rho}{r}\frac{1}{\sigma_{0}}\frac{d\sigma}{d\rho}$. This is done to
highlight and expand the transition region without distorting the curves. For
the green (blue) curve, we vary the shoulder (dijet) profile only, and for the
red curve, we present the envelope of both profiles.

Figure -352: Plots showing various uncertainties associated with the scale
variations. Top left panel shows the fixed order distributions for reference
only (these uncertainties are obtained by varying the renormalization scale by
a factor of 2, which are part of our hard uncertainties and not included in
the overall envelope). The top-right panel shows the hard variation, where
$\mu_{h}$ is varied in the dijet, shoulder and matched contributions
simultaneously. The remaining panels show soft, jet and profile variations for
the shoulder and dijet contributions separately.

Next, we consider results with both shoulder and dijet resummation, and show
in Fig. -353 how the profile functions interpolate between the regions where
shoulder resummation and dijet resummation are important. As discussed in
Section 4.3, we use the same profile functions for both dijet resummation and
shoulder resummation, such that the dijet resummation is active in the range
$0\leq\rho\lesssim 0.2$, and the shoulder resummation is active in the range
$\rho\gtrsim 0.25$, with a smooth transition between them. The parameters used
and their variations are given in Table 1. We use NNLL resummation for the
shoulder, NNLL′ resummation for the dijet, and match to the fixed NLO
($\mathcal{O}(\alpha_{s}^{2})$) distribution. These form a consistent set of
resummed orders since they both contain the complete
$\mathcal{O}(\alpha_{s}^{2})$ information in the respective singular limit. We
see that the matched resummed distribution smoothly merges into the dijet
resummation near $\rho\sim 0.2$, and into the shoulder resummation near
$\rho\sim 0.25$, as desired. In the right panel of Fig. -353, we highlight the
transition region by dividing the distribution by $r=\rho-\frac{1}{3}$. In
this panel we include profile uncertainty estimates obtained by varying the
profile parameters within ranges given in Table 1. In all cases the band shown
is the envelope of the relevant individual variations.

Figure -351: Envelope of all perturbative scale variations in the matched
predictions shown in Fig -352.

We next plot the various uncertainties from individual hard, jet, and soft
scale variations and profile variations in Fig. -352. The first panel shows
the fixed-order predictions. We use Event2 [66, 67] for NLO and CoLoRFulNNLO
[5, 6] for NNLO. The hard scale variation for the shoulder and dijet
resummation is correlated, since there is only a single common hard scale for
both factorizations. We vary the hard scale by a factor of 2 above and below
its central value at 91.2 GeV. The soft scales
($\mu_{s,\text{sh}},\mu_{s,\text{dij}}$) are both independently varied for the
shoulder and the dijet. The range of variation is approximately a factor of 2
relative to their pointwise central value, while simultaneously varying the
jet scale using the relation $\mu_{j}=\sqrt{\mu_{h}\mu_{s}}$. The jet scale
variation is done as per the prescription in Eqs. (4.4) and (4.4). Note that
all variations maintain the monotonicity condition $\mu_{s}<\mu_{j}<\mu_{h}$
(up to the numerically negligible region where $|z|<0.56$), and are thus
sensibly defined. Finally, the profile parameter variations are performed as
summarized in Table 1. We recall that the shoulder profile variation is
performed while fixing the profile for dijet scales and varying only the
shoulder profiles, and vice-versa for the dijet profile variation.

Our final prediction, with uncertainties given by an overall envelope of all
variations, is shown in Fig. -351. Increasing the resummation order leads to a
visible decrease in perturbative uncertainty. At the same time, the NNLL
result is fully covered by the larger uncertainty estimate at NLL. The
dominant source of perturbative uncertainty in both the shoulder and dijet
region is the variation of the respective soft scale.

Figure -350:  Plots comparing our dijet + shoulder resummation to both “full”
dijet resummation (matched to NLO, but no shoulder) and fixed order results at
NLO and NNLO. In the left panels we show results over a wider range of $\rho$
where both shoulder and dijet resummation are active, while in the right
panels we zoom in on values of $\rho$ near the Sudakov shoulder. In the upper
panels we show absolute cross sections, while in the lower panels we show
variations relative to the central curve for our highest order shoulder+dijet
resummed result, namely
$\Delta(\text{resum})\equiv\frac{\text{resum}}{\text{NNLL}^{\prime\text{pro}}_{\text{dij}}+\text{NNLL}^{\text{pro}}_{\text{sh}}+\text{NLO}}-1$.

Finally, we compare our prediction combining shoulder and dijet resummation to
the previous approach where only dijet resummation was included [14] in Fig.
-350. For reference, we also plot the NLO and NNLO distributions in the range
$0.10\lesssim\rho\lesssim 0.4$. Fig. -350 also shows the associated
perturbative uncertainty with each resummation, which we obtain by considering
an envelope of all scale and profile variations. (Note that the resummed
results with and without shoulder resummation use different profile functions,
as described in Section 4.2.) The uncertainty associated with the FO curves is
given by a scale variation by a standard factor of 2 around the central value
$\mu_{\text{FO}}=m_{Z}$. The sharp spikes in the dijet resummation and FO
distributions are indicative of the presence of large shoulder logs, which are
exponentiated by shoulder resummation, leaving behind a smooth result. The
addition of shoulder resummation has a numerically significant impact on the
value of the cross section in the region $\rho\gtrsim 0.2$, which includes
both the transition region between dijet and shoulder resummation
$0.2\lesssim\rho\lesssim 0.25$ and the region where shoulder resummation
dominates, for $\rho\gtrsim 0.25$. These differences will be relevant for
$\alpha_{s}(m_{Z})$ fits using the heavy jet mass distribution. Next, turning
to the region of very large $\rho>\frac{1}{3}$ to the right of the shoulder,
we can see that the fully matched prediction merges onto the fixed NLO result
past $\rho\gtrsim 0.34$, as desired for both the central value and uncertainty
band. Since the NLO is the first non-zero order to contribute in this region,
its $\mathcal{O}(20\%)$ uncertainty is substantially larger than any of the
predictions below the shoulder. Finally, we observe that the fully matched
prediction with shoulder resummation provides a realistic uncertainty estimate
across the shoulder, smoothly connecting the NNLL uncertainties on the left
with the NLO (effectively LO) uncertainties on the right. In contrast, both
fixed order and dijet-resummation-only predictions are discontinuous at the
shoulder and their uncertainty bands can be seen to be unreliable, since they
do not capture the impact of the large logarithmic contributions. From the
NNLO curve in Fig -350 we see that this higher order contribution does not
improve the ability of fixed order results to capture the large resummation
effects in the $0.1<\rho<0.34$ region. Including NNLO results can be expected
to reduce the perturbative uncertainty for $\rho>0.34$, where we note that the
NNLO central value is within the NLO uncertainties.

##  6 Conclusions

In this paper we have considered the resummation of Sudakov shoulder
logarithms in the heavy jet mass distribution in $e^{+}e^{-}$ events. These
logarithms were identified by Catani and Webber [24] in 1997, and a
systematically improvable method to compute and resum these logarithms was
developed in 2022 [26]. When the factorized distribution is resummed in
momentum space, that is the space where the heavy jet mass itself lives,
spurious “Sudakov Landau poles” result. These poles are divergences in the
distribution arising from scale setting and have no physical significance.
Similar poles have been found in factorized expressions at hadron colliders,
such as the $q_{T}$ spectrum in Drell-Yan. This is the first example we know
of for Sudakov Landau poles in a non-transverse momentum distribution.

The Sudakov Landau poles in the heavy jet mass distribution were found in [26]
and prohibited the quantitative application of the resummation formula. In
this paper we have identified the origin of these poles as arising from an
inconsistent solution of the evolution equations: when scales are set in
momentum space, subleading terms in logarithmic power counting are included
which give rise to these singularities. In contrast, when the resummation is
performed in position space and the momentum-space distribution is only
generated by a final Fourier transform after scales have been set, then the
Sudakov Landau poles are removed. We verified this analytically and
numerically.

In addition to the Sudakov Landau pole, another complication in the HJM
shoulder is the existence of a linear UV divergence in $\frac{d\sigma}{d\rho}$
associated with the region where the light and jet mass are both very large,
but their difference is finite. We regulate this divergence by considering the
third-derivative $\frac{d^{3}\sigma}{d\rho^{3}}$, and specifying appropriate
boundary conditions to carry out the double integral, which fix terms of the
form $c_{0}+c_{1}\rho$. In particular we demanded that the resummed and
matched distribution and its first derivative agree with the fixed order
distribution at $r=\frac{1}{3}-\rho=0$, where perturbative uncertainties
provide natural variations in this constraint.

The trijet hemisphere soft function for the Sudakov shoulder is a two-scale
function, like the one for heavy jet mass in the dijet region. Thus it has
non-global structure which is worth exploring. Unlike the hemisphere soft
function, the trijet hemisphere soft function is not symmetric in its
arguments, and there is not a natural way to separate its global and non-
global parts. While one might like to try to choose different complex scales
$\mu_{s\ell}$ and $\mu_{sh}$ for the light and heavy hemispheres, we find such
complex scales are problematic, since it is difficult to ensure logs of the
form $\ln\frac{\mu_{s\ell}}{\mu_{sh}}$ do not induce spurious effects. If
instead, one uses a single real soft scale, then the only missing logarithms
are non-global in nature and integrate to a fixed-order contribution to the
trijet soft function starting at order $\alpha_{s}^{3}$, similar to the
leading non-global logarithmic contribution to the dijet heavy jet mass
distribution. Such contributions are formally N3LL, the same order as the
unknown 2-loop soft function constant. Thus our NNLL calculation is self-
consistent with real scales and dropping the non-global contributions.

To make a quantitative prediction, we needed to match between the dijet and
shoulder regions. This requires us to modulate the hard, jet and soft scales
with profile functions, where as $\rho$ grows we turn off the dijet
resummation and turn on the shoulder resummation. We included such profiles
and provided the first quantitative predictions that include both types of
resummation. We find that at the NNLL level, the Sudakov shoulder resummation
provides as much as a 20% increase over the NLO cross section for values of
$\rho\gtrsim 0.25$, and captures effects that may be numerically important for
$\alpha_{s}(m_{Z})$ fits for $\rho\gtrsim 0.2$.

For a full phenomenological analysis and comparison to data, the main
remaining ingredient needed is the incorporation of non-perturbative power
corrections due to hadronization. In the dijet region, much is understood
about how the power corrections affect the distribution. In the trijet region,
much less is known. It will be important to understand how power corrections
coming from factorization formula in the trijet region can be combined with
power corrections from the dijet region in a manner that goes beyond
hadronization models. Combining the NNLL resummation of Sudakov shoulder
logarithms that we have developed in this paper with a rigorous field
theoretic treatment of power corrections could be critical in resolving long-
standing discrepancies between the values of $\alpha_{s}$ extracted from the
heavy jet mass distribution and the world average.

## Acknowledgments

The authors would like to thank T. Becher, P. Monni, A. Gao and Y. Li for
useful conversations. M.D.S., XY.Z. and A.B. are supported in part by the U.S.
Department of Energy under contract DE-SC0013607. I.S. and J.M. were supported
in part by the U.S. Department of Energy contract DE-SC0011090. I.S. was also
supported in part by the Simons Foundation through the Investigator grant
327942.

##  Appendix A One-loop trijet hemisphere soft function

In this section, we review the calculation of trijet hemisphere soft function
needed for both thrust and hemisphere. Based on the analysis of trijet
kinematics in [26], we define the six-directional differential soft function
via

$S_{6,g/q}(q_{i})=2g_{s}^{2}\mu^{2\epsilon}\iota^{\epsilon}\int\frac{d^{d}k}{(2\pi)^{d-1}}\delta^{+}(k^{2})\hat{\mathscr{M}}(k,q_{i})|\mathcal{M}_{g/q}|^{2},\quad
i=1,2,3,\bar{1},\bar{2},\bar{3}$ (A.1)

with the conventional factor $\iota=\frac{e^{\gamma_{E}}}{4\pi}$ and the
eikonal matrix elements

$\displaystyle|\mathcal{M}_{g}|^{2}=\left[\left(C_{F}-\frac{1}{2}C_{A}\right)\frac{n_{2}\cdot
n_{3}}{(n_{2}\cdot k)(n_{3}\cdot k)}+\frac{1}{2}C_{A}\frac{n_{1}\cdot
n_{2}}{(n_{1}\cdot k)(n_{2}\cdot k)}+\frac{1}{2}C_{A}\frac{n_{1}\cdot
n_{3}}{(n_{1}\cdot k)(n_{3}\cdot k)}\right],$
$\displaystyle|\mathcal{M}_{q}|^{2}=\left[\left(C_{F}-\frac{1}{2}C_{A}\right)\frac{n_{1}\cdot
n_{3}}{(n_{1}\cdot k)(n_{3}\cdot k)}+\frac{1}{2}C_{A}\frac{n_{1}\cdot
n_{2}}{(n_{1}\cdot k)(n_{2}\cdot k)}+\frac{1}{2}C_{A}\frac{n_{2}\cdot
n_{3}}{(n_{2}\cdot k)(n_{3}\cdot k)}\right],$ (A.2)

and the soft measurement function ($\hat{\mathscr{M}}$) at one-loop

$\displaystyle\hat{\mathscr{M}}(k,q_{i})$ $\displaystyle=\theta(n_{2}\cdot
k-\bar{n}_{2}\cdot k)\theta(n_{3}\cdot k-\bar{n}_{3}\cdot
k)\delta\left(q_{j\neq 1}\right)\delta\left(q_{1}-\frac{2}{3}n_{1}\cdot
k\right)$ $\displaystyle+\theta(n_{1}\cdot k-\bar{n}_{1}\cdot
k)\theta(n_{3}\cdot k-\bar{n}_{3}\cdot k)\delta\left(q_{j\neq
2}\right)\delta\left(q_{2}-\frac{2}{3}n_{2}\cdot k\right)$
$\displaystyle+\theta(n_{1}\cdot k-\bar{n}_{1}\cdot k)\theta(n_{2}\cdot
k-\bar{n}_{2}\cdot k)\delta\left(q_{j\neq
3}\right)\delta\left(q_{3}-\frac{2}{3}n_{3}\cdot k\right)$
$\displaystyle+\theta(\bar{n}_{3}\cdot k-n_{3}\cdot k)\theta(\bar{n}_{2}\cdot
k-n_{2}\cdot
k)\delta\left(q_{j\neq\bar{1}}\right)\delta\left(q_{\bar{1}}-\frac{2}{3}\bar{n}_{1}\cdot
k\right)$ $\displaystyle+\theta(\bar{n}_{1}\cdot k-n_{1}\cdot
k)\theta(\bar{n}_{3}\cdot k-n_{3}\cdot
k)\delta\left(q_{j\neq\bar{2}}\right)\delta\left(q_{\bar{2}}-\frac{2}{3}N_{2}\cdot
k\right)$ $\displaystyle+\theta(\bar{n}_{1}\cdot k-n_{1}\cdot
k)\theta(\bar{n}_{2}\cdot k-n_{2}\cdot
k)\delta\left(q_{j\neq\bar{3}}\right)\delta\left(q_{\bar{3}}-\frac{2}{3}N_{3}\cdot
k\right)\,.$ (A.3)

Here $n_{1}$, $n_{2}$ and $n_{3}$ are three lightlike vectors that describe
the trijet configuration. For concreteness, we align

$\displaystyle n_{1}$ $\displaystyle=\left(1,0,0,1\right),$ $\displaystyle
n_{2}$ $\displaystyle=\left(1,0,\frac{\sqrt{3}}{2},-\frac{1}{2}\right),$
$\displaystyle n_{3}$
$\displaystyle=\left(1,0,-\frac{\sqrt{3}}{2},-\frac{1}{2}\right)\ .$ (A.4)

For thrust, the last two projections $N_{2,3}$ are lightlike while for HJM, it
is spacelike. Explicitly,

$\displaystyle\text{Thrust}:$ $\displaystyle\quad N_{2}=\bar{n}_{2},\quad
N_{3}=\bar{n}_{3}$ $\displaystyle\text{HJM}:$ $\displaystyle\quad
N_{2}=\left(1,0,\frac{\sqrt{3}}{2},\frac{3}{2}\right),\quad
N_{3}=\left(1,0,-\frac{\sqrt{3}}{2},\frac{3}{2}\right)\,.$ (A.5)

The short-hand notation $\delta(q_{j\neq 1})$ stands for
$\delta(q_{2})\delta(q_{3})\delta(q_{\bar{1}})\delta(q_{\bar{2}})\delta(q_{\bar{3}})$.

In the following, we will focus on the HJM case. To simplify the calculation,
we perform the standard topology identification. First of all, let us
introduce the notation

$I_{n_{a},n_{b},n_{c},n_{d},N}=\int d^{d}k\frac{n_{a}\cdot n_{b}}{(n_{a}\cdot
k)(n_{b}\cdot k)}\delta^{+}(k^{2})\delta\left(q-\frac{2}{3}N\cdot
k\right)\theta\left(n_{c}\cdot k-\bar{n}_{c}\cdot
k\right)\theta\left(n_{d}\cdot k-\bar{n}_{d}\cdot k\right)$ (A.6)

then our soft function can be written as a linear combination of
$I_{n_{a},n_{b},n_{c},n_{d},N}$. For example, the first term in Eq. (A) can be
mapped to

$S_{6,g}(q_{i})\supset\left(C_{F}-\frac{1}{2}C_{A}\right)I_{n_{2},n_{3},n_{2},n_{3},n_{1}}+\frac{1}{2}C_{A}I_{n_{1},n_{2},n_{2},n_{3},n_{1}}+\frac{1}{2}C_{A}I_{n_{1},n_{3},n_{2},n_{3},n_{1}}\,.$
(A.7)

Using rotations and reflection, we can rename the projection vectors and
reduce the soft function into 7 master integrals:

$\displaystyle I_{0}=I_{n_{1},n_{2},n_{2},n_{3},n_{1}},\quad
I_{1}=I_{n_{2},n_{3},n_{2},n_{3},n_{1}},\quad
I_{2}=I_{n_{1},n_{3},\bar{n}_{2},\bar{n}_{3},\bar{n}_{1}},\quad
I_{3}=I_{n_{2},n_{3},\bar{n}_{2},\bar{n}_{3},\bar{n}_{1}},\quad$
$\displaystyle I_{4}=I_{n_{2},n_{3},\bar{n}_{1},\bar{n}_{2},N_{3}},\quad
I_{5}=I_{n_{1},n_{2},\bar{n}_{1},\bar{n}_{2},N_{3}},\quad
I_{6}=I_{n_{1},n_{3},\bar{n}_{1},\bar{n}_{2},N_{3}}\,.$ (A.8)

To evaluate the integrals, we use the lightcone parameterization,
appropriately expand to the first order in $d=4-2\epsilon$ and compute each
piece numerically. It turns out some of the numbers can be reconstructed with
the pentagon table [68] via PSLQ [69, 70]. $I_{0}$ is the only divergent
integral, which reads

$\displaystyle I_{0}$
$\displaystyle=\mathcal{N}\bigg{[}\frac{4\pi}{3\epsilon}+\frac{1}{3}\pi\ln\frac{81}{64}-2\kappa+\epsilon\bigg{(}\frac{24c_{1}}{5}+\frac{32c_{2}}{3}+\frac{10}{9}\pi\text{Li}_{2}\left(-\frac{1}{2}\right)-\frac{91\pi^{3}}{270}+\frac{17}{30}\pi\ln^{2}3$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad-\frac{25}{9}\pi\ln^{2}2-2\pi\ln
2\ln 3\bigg{)}\bigg{]}\,.$ (A.9)

For others, we get

$\displaystyle I_{1}$
$\displaystyle=\mathcal{N}\bigg{[}4\kappa-\frac{4}{3}\pi\ln 2+\epsilon
c_{3}\bigg{]},\quad I_{2}=\mathcal{N}\bigg{[}-2\kappa+2\pi\ln 2+\epsilon
c_{4}\bigg{]}\,,$ $\displaystyle I_{3}$
$\displaystyle=\mathcal{N}\bigg{[}4\kappa+\frac{4}{3}\pi\ln 2+\epsilon
c_{5}\bigg{]},\quad I_{4}=\mathcal{N}\bigg{[}-2\kappa+2\pi\ln 2+\epsilon
c_{6}\bigg{]}\,,$ $\displaystyle I_{5}$
$\displaystyle=\mathcal{N}\bigg{[}4\kappa+\frac{4}{3}\pi\ln 2+\epsilon
c_{7}\bigg{]},\quad I_{6}=\mathcal{N}\bigg{[}-2\kappa+2\pi\ln 2+\epsilon
c_{8}\bigg{]}\,,$ (A.10)

where the constants $\kappa$ and $c_{i}$ are

$\displaystyle\kappa=\Im\left[\text{Li}_{2}\left(e^{\frac{i\pi}{3}}\right)\right],\quad
c_{1}=\Im\left[\text{Li}_{3}\left(\frac{i}{\sqrt{3}}\right)\right],\quad
c_{2}=\Im\left[\text{Li}_{3}\left(1+i\sqrt{3}\right)\right],\quad
c_{3}=-1.89958,$ $\displaystyle c_{4}=-3.83452,\quad c_{5}=-12.3488,\quad
c_{6}=5.56704,\quad c_{7}=15.9482,\quad c_{8}=3.52263,$ (A.11)

and the overall normalization is

$\mathcal{N}=\frac{3}{8}\Omega_{d-3}\frac{1}{q^{1+2\epsilon}}\left(\frac{2}{3}\right)^{2\epsilon},\quad\Omega_{d}\equiv\frac{2\pi^{d/2}}{\Gamma(d/2)}\,.$
(A.12)

Putting back the color factors and renormalizing in $\overline{\text{MS}}$, we
obtain the differential soft function defined in Eq. (A.1). To get the trijet
hemisphere soft function, we need to integrate the soft momenta over the
hemispheres, i.e.

$\displaystyle S_{i}(q_{L},q_{H},\mu)$ $\displaystyle=\int
d^{6}q_{i}S_{6i}(q_{i},\mu)\delta\left(q_{L}-q_{1}-q_{\bar{2}}-q_{\bar{3}}\right)\delta\left(q_{H}-q_{\bar{1}}-q_{2}-q_{3}\right)$
$\displaystyle\equiv S_{iL}(q_{L},\mu)S_{iH}(q_{H},\mu)S_{f}(q_{L}/q_{H})\,.$
(A.13)

Note that for convenience, we split the soft constant into $S_{iL}(q_{L},\mu)$
and $S_{iH}(q_{H},\mu)$ and leave only non-global logarithmic structure in
$S_{f}(q_{L}/q_{H})$. Since non-global logarithms are the same order as the
finite part of the 2-loop soft function, which is relevant only for N3LL
resummation, we simply set $S_{f}(q_{L}/q_{H})=1$. The final result is

$\displaystyle S_{gL}(q_{L},\mu)$
$\displaystyle=\delta(q_{L})\left(1+\frac{\alpha_{s}}{4\pi}c_{gL}\right)+\frac{\alpha_{s}}{4\pi}\bigg{[}\frac{-2C_{A}\Gamma_{0}\ln\frac{q_{L}}{\mu}+2\gamma_{sg}}{q_{L}}\bigg{]}_{\star}\,,$
$\displaystyle S_{gH}(q_{H},\mu)$
$\displaystyle=\delta(q_{H})\left(1+\frac{\alpha_{s}}{4\pi}c_{gH}\right)+\frac{\alpha_{s}}{4\pi}\bigg{[}\frac{-4C_{F}\Gamma_{0}\ln\frac{q_{H}}{\mu}+2\gamma_{sqq}}{q_{H}}\bigg{]}_{\star}\,,$
$\displaystyle S_{qL}(q_{L},\mu)$
$\displaystyle=\delta(q_{L})\left(1+\frac{\alpha_{s}}{4\pi}c_{qL}\right)+\frac{\alpha_{s}}{4\pi}\bigg{[}\frac{-2C_{F}\Gamma_{0}\ln\frac{q_{L}}{\mu}+2\gamma_{sq}}{q_{L}}\bigg{]}_{\star}\,,$
$\displaystyle S_{qH}(q_{H},\mu)$
$\displaystyle=\delta(q_{H})\left(1+\frac{\alpha_{s}}{4\pi}c_{qH}\right)+\frac{\alpha_{s}}{4\pi}\bigg{[}\frac{-2(C_{F}+C_{A})\Gamma_{0}\ln\frac{q_{H}}{\mu}+2\gamma_{sqg}}{q_{H}}\bigg{]}_{\star}\,,$
(A.14)

with the non-cusp anomalous dimensions

$\displaystyle\gamma_{sg}$ $\displaystyle=8C_{F}\ln 2-4C_{A}\ln
3,\quad\gamma_{sqq}=-8C_{F}\ln 6\,,$ $\displaystyle\gamma_{sq}$
$\displaystyle=4C_{A}\ln
2+4C_{F}\ln\frac{2}{3},\quad\gamma_{sqg}=-4(C_{A}+C_{F})\ln 6\,,$ (A.15)

and the soft constants

$\displaystyle c_{gL}$ $\displaystyle=-4.44002C_{A}+1.68285C_{F},$
$\displaystyle c_{gH}$ $\displaystyle=-0.210218C_{A}+5.10882C_{F},$
$\displaystyle c_{qL}$ $\displaystyle=0.841426C_{A}-3.59860C_{F},$
$\displaystyle c_{qH}$ $\displaystyle=2.55411C_{A}+2.34389C_{F}\,.$ (A.16)

For completeness, we can also extract the exact 1-loop soft function for
thrust in the trijet region. The only difference between the thrust and HJM
trijet soft functions is that the soft constants in the light hemisphere are
different:

$\displaystyle c_{gL}^{\text{thrust}}$
$\displaystyle=8.09467C_{A}+10.6607C_{F},\quad
c_{gH}^{\text{thrust}}=-0.210218C_{A}+5.10882C_{F},$ $\displaystyle
c_{qL}^{\text{thrust}}$ $\displaystyle=5.33034C_{A}+13.4250C_{F},\quad
c_{qH}^{\text{thrust}}=2.55411C_{A}+2.34389C_{F}\,.$ (A.17)

##  Appendix B Hard, jet and soft functions

In this Appendix, we summarize the hard and jet functions. We write the cusp
anomalous dimension as

$\displaystyle\Gamma_{\text{cusp}}$
$\displaystyle=\left(\frac{\alpha_{s}}{4\pi}\right)\Gamma_{0}+\left(\frac{\alpha_{s}}{4\pi}\right)^{2}\Gamma_{1}+\left(\frac{\alpha_{s}}{4\pi}\right)^{3}\Gamma_{2}+\cdots\,,$

where [71, 72]

$\displaystyle\Gamma_{0}$ $\displaystyle=4\,,$ $\displaystyle\Gamma_{1}$
$\displaystyle=4\left[C_{A}\left(\frac{67}{9}-\frac{\pi^{2}}{3}\right)-\frac{20}{9}T_{F}n_{F}\right]\,,$
$\displaystyle\Gamma_{2}$
$\displaystyle=4\bigg{[}\left(-\frac{56\zeta_{3}}{3}-\frac{418}{27}+\frac{40\pi^{2}}{27}\right)C_{A}n_{f}T_{F}+\left(\frac{22\zeta_{3}}{3}+\frac{245}{6}-\frac{134\pi^{2}}{27}+\frac{11\pi^{4}}{45}\right)C_{A}^{2}$
$\displaystyle\qquad+\left(16\zeta_{3}-\frac{55}{3}\right)C_{F}n_{f}T_{F}-\frac{16}{27}n_{f}^{2}T_{F}^{2}\bigg{]}\,.$
(B.1)

#### Hard Functions

Hard functions for event shape factorizations in SCET are obtained by matching
QCD form factors onto corresponding SCET operators. In general, the hard
function RGE takes the form

$\displaystyle\frac{d\ln
H(Q^{2},\mu^{2})}{d\ln\mu}=2\Gamma^{H}\ln\frac{Q^{2}}{\mu^{2}}+2\gamma^{H}\,,$
(B.2)

where $\Gamma^{H}$ is the anomalous dimension proportional to the cusp
anomalous dimension, and $\gamma_{H}$ depends on the relevant factorization
and kinematics. The solution to the above RGE is given by

$\displaystyle
H(\mu^{2},Q^{2})=\exp\left[4S(\mu_{h},\mu)-2A_{\gamma_{H}}(\mu_{h},\mu)\right]\left(\frac{Q^{2}}{\mu_{h}^{2}}\right)^{-2A_{\Gamma}(\mu_{h},\mu)}\times
h(L)\,,$ (B.3)

where $L\equiv\ln\left(\frac{Q^{2}}{\mu_{h}^{2}}\right)$ is the fixed order
hard function at the scale $\mu_{h}$. Perturbatively, we have

$\displaystyle h(L)$
$\displaystyle=1+\left(\frac{\alpha_{s}}{4\pi}\right)\left[-\Gamma_{0}^{H}\frac{L^{2}}{2}-\gamma_{0}^{H}L+c_{1}^{H}\right]+\left(\frac{\alpha_{s}}{4\pi}\right)^{2}\bigg{[}(\Gamma_{0}^{H})^{2}\frac{L^{4}}{8}+(\beta_{0}+3\gamma_{0}^{H})\Gamma_{0}^{H}\frac{L^{3}}{6}$
$\displaystyle+\left(-\Gamma_{1}^{H}+(\gamma_{0}^{H})^{2}+\beta_{0}\gamma_{0}^{H}-c_{1}^{H}\Gamma_{0}^{H}\right)\frac{L^{2}}{2}+(-\gamma_{1}^{H}-\gamma_{0}^{H}c_{1}^{H}-\beta_{0}c_{1}^{H})L+c_{2}^{H}\bigg{]}\,.$
(B.4)

The constants $c_{i}^{H}$ encode the fixed order, non-log pieces in the hard
function. At our working order, the hard anomalous dimension has the form [73]

$\displaystyle\frac{d}{d\ln{\mu}}\ln{H(Q^{2},\mu^{2})}\
=-\Gamma_{\mathrm{cusp}}\sum_{i<j}\mathbf{T}_{i}\cdot\mathbf{T}_{j}\ln{\left[\frac{2q_{i}\cdot
q_{j}}{\mu^{2}}\right]}+2\gamma^{\prime}_{H}\,,$ (B.5)

where $\mathbf{T}_{j}^{a}$ denotes the color charge of the $j^{\mathrm{th}}$
external hard parton, written using color space formalism [74]. Here $q_{i}$
is the momentum of the outgoing $i^{\mathrm{th}}$ hard parton. For the
shoulder resummation, the hard kinematics is a trijet one, equally spaced in
angle given by

$\displaystyle q_{i}^{\mu}$ $\displaystyle=\frac{Q}{3}n^{i}_{\mu}\ ,$
$\displaystyle 2q_{i}\cdot q_{j}$ $\displaystyle=\frac{Q^{2}}{3}\ .$ (B.6)

with the relevant partons being $q,\bar{q},g$. The relevant color sums are

$\displaystyle\mathbf{T}_{q}\cdot\mathbf{T}_{\bar{q}}$
$\displaystyle=\frac{C_{A}}{2}-C_{F}\,,$
$\displaystyle\mathbf{T}_{q}\cdot\mathbf{T}_{g}$
$\displaystyle=\mathbf{T}_{\bar{q}}\cdot\mathbf{T}_{g}=-\frac{C_{A}}{2}\,.$
(B.7)

The non-cusp hard anomalous dimension $\gamma^{\prime}_{H}$ is also color
diagonal up to two loop order, and is thus a sum of quark and gluon pieces,
(i.e)

$\displaystyle\gamma^{\prime}_{H}=2\gamma^{q}+\gamma^{g}\,.$ (B.8)

These are extracted from the expressions in [75, 47] and are given by

$\displaystyle\gamma^{q}=\left(\frac{\alpha_{s}}{4\pi}\right)\gamma_{0}^{q}+\left(\frac{\alpha_{s}}{4\pi}\right)^{2}\gamma^{q}_{1}+\cdots\,,\qquad\gamma^{g}=\left(\frac{\alpha_{s}}{4\pi}\right)\gamma_{0}^{g}+\left(\frac{\alpha_{s}}{4\pi}\right)^{2}\gamma^{g}_{1}+\cdots\,,$
(B.9)

where [76, 77, 78, 79, 80]

$\displaystyle\gamma^{q}_{0}$ $\displaystyle=-6C_{F}\,,$ (B.10)
$\displaystyle\gamma^{g}_{0}$ $\displaystyle=-2\beta_{0}\,,$
$\displaystyle\gamma^{q}_{1}$
$\displaystyle=C_{F}^{2}\left(4\pi^{2}-3-48\zeta_{3}\right)+C_{F}C_{A}\left(-\frac{11\pi^{2}}{3}-\frac{961}{27}+52\zeta_{3}\right)+C_{F}T_{F}n_{f}\left(\frac{4\pi^{2}}{3}+\frac{260}{27}\right)\,,$
$\displaystyle\gamma^{g}_{1}$
$\displaystyle=C_{A}^{2}\left(\frac{11\pi^{2}}{9}-\frac{1384}{27}+4\zeta_{3}\right)+C_{A}T_{F}n_{f}\left(-\frac{4\pi^{2}}{9}+\frac{512}{27}\right)+8C_{F}T_{F}n_{f}\,.$

Using Eqns. (B.6), (B), (B.8) in (B.5), we find that the anomalous dimension
$\Gamma^{H}$ for the shoulder trijet hard function is

$\displaystyle\Gamma^{H}$
$\displaystyle=\left(C_{F}+\frac{1}{2}C_{A}\right)\Gamma_{\text{cusp}}\,,$
(B.11) $\displaystyle\gamma^{H}$
$\displaystyle=\left(\frac{\alpha_{s}}{4\pi}\right)\gamma_{0}^{H}+\left(\frac{\alpha_{s}}{4\pi}\right)^{2}\gamma_{1}^{H}+\cdots\
,$

where

$\displaystyle\gamma_{0}^{H}$
$\displaystyle=-\beta_{0}-6C_{F}-2(C_{A}+2C_{F})\ln 3\,,$
$\displaystyle\gamma_{1}^{H}$
$\displaystyle=C_{A}T_{F}n_{f}\left(\frac{256}{27}-\frac{2}{9}\pi^{2}+\frac{40}{9}\ln
3\right)+C_{A}^{2}\left(2\zeta_{3}+\frac{11}{18}\pi^{2}-\frac{692}{27}+\frac{2\pi^{2}}{3}\ln
3-\frac{134}{9}\ln 3\right)$
$\displaystyle\quad+C_{F}T_{F}n_{f}\left(\frac{368}{27}+\frac{4\pi^{2}}{3}+\frac{80}{9}\ln
3\right)+C_{F}^{2}\left(-48\zeta_{3}+4\pi^{2}-3\right)$
$\displaystyle\quad+C_{F}C_{A}\left(52\zeta_{3}-\frac{11}{3}\pi^{2}-\frac{961}{27}+\frac{4\pi^{2}}{3}\ln
3-\frac{268}{9}\ln 3\right)\,.$ (B.12)

For NNLL resummation, we also need the $c_{i}^{H}$ constants to one-loop,
namely $c_{1}^{H}$. It can be extracted from the real-virtual correction of
$\gamma^{\star}\rightarrow q\bar{q}$ squared matrix elements [42, 81],

$\displaystyle c_{1}^{H}$
$\displaystyle=C_{F}\left[-\frac{65}{4}+\frac{3\pi^{2}}{2}-\frac{21}{8}\ln
3-10\ln 2\ln 3+3\ln^{2}3+10\text{Li}_{2}\left(\frac{1}{3}\right)\right]$
$\displaystyle\quad+C_{A}\left[\frac{3}{4}+\frac{5\pi^{2}}{4}+\frac{3}{8}\ln
3+\ln 2\ln 3-\frac{3}{2}\ln^{2}3-\text{Li}_{2}\left(\frac{1}{3}\right)\right]\
.$ (B.13)

The hard anomalous dimensions are also common to other hard functions that
have appeared earlier in the SCET literature [41, 40], and we have verified
ours against the same.

#### Jet functions

The jet function RGE in the Laplace space is

$\frac{d\ln\tilde{j}(\nu)}{d\ln\mu}=-2\Gamma^{J}\ln\frac{1}{\nu\mu^{2}e^{\gamma_{E}}}-2\gamma^{J}\,,$
(B.14)

with perturbative solution

$\displaystyle\tilde{j}(L)$
$\displaystyle=1+\left(\frac{\alpha_{s}}{4\pi}\right)\left[\Gamma_{0}^{J}\frac{L^{2}}{2}+\gamma_{0}^{J}L+c_{1}^{J}\right]+\left(\frac{\alpha_{s}}{4\pi}\right)^{2}\bigg{[}(\Gamma_{0}^{J})^{2}\frac{L^{4}}{8}+(-\beta_{0}+3\gamma_{0}^{J})\Gamma_{0}^{J}\frac{L^{3}}{6}$
$\displaystyle+(\Gamma_{1}^{J}+(\gamma_{0}^{J})^{2}-\beta_{0}\gamma_{0}^{J}+c_{1}^{J}\Gamma_{0}^{J})\frac{L^{2}}{2}+(\gamma_{1}^{J}+\gamma_{0}^{J}c_{1}^{J}-\beta_{0}c_{1}^{J})+c_{2}^{J}\bigg{]}\,.$
(B.15)

We will need the $\mathcal{O}(\alpha_{s})$ piece for NNLL shoulder resummation
and $\mathcal{O}(\alpha_{s}^{2})$ piece for NNLL′ dijet resummation. For quark
jets [82, 83, 38], the anomalous dimension and constant are

$\displaystyle\Gamma^{J_{q}}$ $\displaystyle=C_{F}\Gamma_{\text{cusp}}\,,$
$\displaystyle\gamma_{0}^{J_{q}}$ $\displaystyle=-3C_{F}\,,$
$\displaystyle\gamma_{1}^{J_{q}}$
$\displaystyle=C_{F}^{2}\left(-\frac{3}{2}+2\pi^{2}-24\zeta_{3}\right)+C_{F}C_{A}\left(-\frac{1759}{54}-\frac{11}{9}\pi^{2}+40\zeta_{3}\right)+C_{F}T_{F}n_{f}\left(\frac{242}{27}+\frac{4}{9}\pi^{2}\right),$
$\displaystyle c_{1}^{J_{q}}$
$\displaystyle=C_{F}\left(7-\frac{2}{3}\pi^{2}\right)\,.$ (B.16)

For gluon jets, we have [40, 39]

$\displaystyle\Gamma^{J_{g}}$ $\displaystyle=C_{A}\Gamma_{\text{cusp}}\,,$
$\displaystyle\gamma_{0}^{J_{q}}$ $\displaystyle=-\beta_{0}\,,$
$\displaystyle\gamma_{1}^{J_{g}}$
$\displaystyle=C_{A}^{2}\left(-\frac{1096}{27}+\frac{11}{9}\pi^{2}+16\zeta_{3}\right)+C_{A}T_{F}n_{f}\left(\frac{368}{27}-\frac{4}{9}\pi^{2}\right)+4C_{F}T_{F}n_{f},$
$\displaystyle c_{1}^{J_{g}}$
$\displaystyle=C_{A}\left(\frac{67}{9}-\frac{2}{3}\pi^{2}\right)-\frac{20}{9}T_{F}n_{f}\,.$
(B.17)

#### Soft functions

The RGE for the trijet hemisphere soft function is

$\displaystyle\frac{d\ln\tilde{s}(\nu)}{d\ln\mu}=-2\Gamma^{S}\ln\frac{1}{\nu\mu
e^{\gamma_{E}}}-2\gamma^{S}\,,$ (B.18)

and its solution is

$\displaystyle\tilde{s}(L)$
$\displaystyle=1+\left(\frac{\alpha_{s}}{4\pi}\right)\left[\Gamma_{0}^{S}\frac{L^{2}}{2}+2\gamma_{0}^{S}L+c_{1}^{S}\right]+\left(\frac{\alpha_{s}}{4\pi}\right)^{2}\bigg{[}(\Gamma_{0}^{S})^{2}\frac{L^{4}}{8}+(-\beta_{0}+3\gamma_{0}^{S})\Gamma_{0}^{S}\frac{L^{3}}{3}$
$\displaystyle+(\Gamma_{1}^{S}+4(\gamma_{0}^{S})^{2}-4\beta_{0}\gamma_{0}^{S}+c_{1}^{S}\Gamma_{0}^{S})\frac{L^{2}}{2}+2(\gamma_{1}^{S}+\gamma_{0}^{S}c_{1}^{S}-\beta_{0}c_{1}^{S})L+c_{2}^{S}\bigg{]}\,.$
(B.19)

The soft anomalous dimension for each hemisphere is

$\displaystyle\Gamma^{sg}$ $\displaystyle=-2C_{A}\Gamma_{\text{cusp}}\,,$
$\displaystyle\quad\gamma_{0}^{sg}$ $\displaystyle=-2C_{A}\ln 3+4C_{F}\ln
2\,,$ $\displaystyle\Gamma^{sq}$
$\displaystyle=-2C_{F}\Gamma_{\text{cusp}}\,,$
$\displaystyle\quad\gamma_{0}^{sq}$
$\displaystyle=-2C_{F}\ln\frac{3}{2}+2C_{A}\ln 2\,,$
$\displaystyle\Gamma^{sqg}$
$\displaystyle=-2(C_{F}+C_{A})\Gamma_{\text{cusp}}\,,$
$\displaystyle\quad\gamma_{0}^{sqg}$ $\displaystyle=-2(C_{A}+C_{F})\ln 6\,,$
$\displaystyle\Gamma^{sqq}$ $\displaystyle=-4C_{F}\Gamma_{\text{cusp}}\,,$
$\displaystyle\quad\gamma_{0}^{sqq}$ $\displaystyle=-4C_{F}\ln 6\,.$ (B.20)

The two-loop non-cusp anomalous dimension can be obtained from RG invariance

$\displaystyle\gamma_{1}^{sg}+\gamma_{1}^{sqq}=\gamma_{1}^{sq}+\gamma_{1}^{sqg}$
$\displaystyle=C_{A}T_{F}n_{f}\left(-\frac{112}{27}+\frac{2\pi^{2}}{9}+\frac{40}{9}\ln
3\right)+C_{F}T_{F}n_{f}\left(-\frac{224}{27}+\frac{4\pi^{2}}{9}+\frac{80}{9}\ln
3\right)$
$\displaystyle\quad+C_{F}C_{A}\left(\frac{808}{27}-\frac{11}{9}\pi^{2}-\frac{268}{9}\ln
3+\frac{4\pi^{2}}{3}\ln 3-28\zeta_{3}\right)$
$\displaystyle\quad+C_{A}^{2}\left(\frac{404}{27}-\frac{11}{18}\pi^{2}-\frac{134}{9}\ln
3+\frac{2\pi^{2}}{3}\ln 3-14\zeta_{3}\right)\,.$ (B.21)

The soft constants were computed in Appendix A. In Laplace space and in the
above notation they translate to

$\displaystyle c_{1}^{sg}=\Gamma_{0}^{sg}\frac{\pi^{2}}{12}+c_{gL},\quad
c_{1}^{sqq}=\Gamma_{0}^{sqq}\frac{\pi^{2}}{12}+c_{gH},\quad
c_{1}^{sq}=\Gamma_{0}^{sq}\frac{\pi^{2}}{12}+c_{qL},\quad
c_{1}^{sqg}=\Gamma_{0}^{sqg}\frac{\pi^{2}}{12}+c_{qH}\,.$ (B.22)

##  Appendix C Running coupling and Sudakov kernels

We expand the QCD beta-function as

$\frac{\mathrm{d}\alpha_{s}(\mu)}{\mathrm{d}\ln\mu}=\beta(\alpha_{s}(\mu)),\quad\beta(\alpha)=-2\alpha\,\left[\left(\frac{\alpha}{4\pi}\right)\beta_{0}+\left(\frac{\alpha}{4\pi}\right)^{2}\beta_{1}+\left(\frac{\alpha}{4\pi}\right)^{3}\beta_{2}+\cdots\right]\,,$
(C.1)

where the coefficient up to three loops are given by [84, 85, 86, 87]

$\displaystyle\beta_{0}$
$\displaystyle=\frac{11}{3}C_{A}-\frac{4}{3}T_{F}n_{f}\,,$ (C.2)
$\displaystyle\beta_{1}$
$\displaystyle=\frac{34}{3}C_{A}^{2}-\frac{20}{3}C_{A}T_{F}n_{f}-4C_{F}T_{F}n_{f}\,,$
$\displaystyle\beta_{2}$
$\displaystyle=n_{f}^{2}T_{F}^{2}\left(\frac{158}{27}C_{A}+\frac{44}{9}C_{F}\right)+n_{f}T_{F}\left(2C_{F}^{2}-\frac{205}{9}C_{F}C_{A}-\frac{1415}{27}C_{A}^{2}\right)+\frac{2857}{54}C_{A}^{3}\,.$

It is also useful sometimes to have an expression for the coupling at one
scale in terms of the coupling at a different scale. For two-loop running,
such an expression is

$\alpha_{s}(\mu)=\alpha_{s}(Q)\left[X+\alpha_{s}(Q)\frac{\beta_{1}}{4\pi\beta_{0}}\ln
X\right]^{-1},\quad X\equiv
1+\frac{\alpha_{s}(Q)}{2\pi}\beta_{0}\ln\frac{\mu}{Q}\,,$ (C.3)

and at three loops

$\alpha_{s}(\mu)=\alpha_{s}(Q)\left\\{X+\alpha_{s}(Q)\frac{\beta_{1}}{4\pi\beta_{0}}\ln
X+\frac{\alpha_{s}^{2}(Q)}{16\pi^{2}}\left[\frac{\beta_{2}}{\beta_{0}}\left(1-\frac{1}{X}\right)+\frac{\beta_{1}^{2}}{\beta_{0}^{2}}\left(\frac{1}{X}-1+\frac{\ln
X}{X}\right)\right]\right\\}^{-1}\,.$ (C.4)

For reference, we also provide explicit expressions [88, 37] for the Sudakov
RG kernels appearing in the evolution factors. Their definitions are:

$\displaystyle S_{\Gamma}(\nu,\mu)$
$\displaystyle=-\int_{\alpha_{s}(\nu)}^{\alpha_{s}(\mu)}d\alpha\
\frac{\Gamma_{\text{cusp}}(\alpha)}{\beta(\alpha)}\int_{\alpha_{s}(\nu)}^{\alpha}\frac{d\alpha^{\prime}}{\beta(\alpha^{\prime})}\,,$
$\displaystyle A_{\Gamma}(\nu,\mu)$
$\displaystyle=-\int_{\alpha_{s}(\nu)}^{\alpha_{s}(\mu)}d\alpha\
\frac{\Gamma_{\text{cusp}}(\alpha)}{\beta(\alpha)}\,,$ $\displaystyle
A_{\gamma_{j}}(\nu,\mu)$
$\displaystyle=-\int_{\alpha_{s}(\nu)}^{\alpha_{s}(\mu)}d\alpha\
\frac{\gamma(\alpha)}{\beta(\alpha)}\,.$ (C.5)

For $S_{\Gamma}(\nu,\mu)$, we have

$S_{\Gamma}(\nu,\mu)=S_{\Gamma}^{\text{LL}}(\nu,\mu)+S_{\Gamma}^{\text{NLL}}(\nu,\mu)+S_{\Gamma}^{\text{NNLL}}(\nu,\mu)+\cdots$
(C.6)

with

$\displaystyle S_{\Gamma}^{\text{LL}}(\nu,\mu)$
$\displaystyle=\frac{\Gamma_{0}\pi}{\beta_{0}^{2}\alpha_{s}(\nu)}\left(1-\frac{1}{r}-\ln
r\right)\,,$ $\displaystyle S_{\Gamma}^{\text{NLL}}(\nu,\mu)$
$\displaystyle=\frac{\Gamma_{0}}{4\beta_{0}^{2}}\left[\left(\frac{\Gamma_{1}}{\Gamma_{0}}-\frac{\beta_{1}}{\beta_{0}}\right)(1-r+\ln
r)+\frac{\beta_{1}}{2\beta_{0}}\ln^{2}(r)\right]\,,$ $\displaystyle
S_{\Gamma}^{\text{NNLL}}(\nu,\mu)$
$\displaystyle=\frac{\Gamma_{0}\alpha_{s}(\nu)}{32\pi\beta_{0}^{2}}\bigg{[}\left(\frac{\beta_{1}^{2}}{\beta_{0}^{2}}-\frac{\beta_{2}}{\beta_{0}}\right)\left(1-r^{2}+2\ln
r\right)+\left(\frac{\beta_{1}\Gamma_{1}}{\beta_{0}\Gamma_{0}}-\frac{\Gamma_{2}}{\Gamma_{0}}\right)(1-r)^{2}$
$\displaystyle+2\left(\frac{\beta_{1}\Gamma_{1}}{\beta_{0}\Gamma_{0}}-\frac{\beta_{1}^{2}}{\beta_{0}^{2}}\right)(1-r+r\ln
r)\bigg{]}\,.$ (C.7)

Here we define $r\equiv\frac{\alpha_{s}(\mu)}{\alpha_{s}(\nu)}$. Similarly, we
have

$\displaystyle
A_{\Gamma}(\nu,\mu)=A_{\Gamma}^{\text{LL}}(\nu,\mu)+A_{\Gamma}^{\text{NLL}}(\nu,\mu)+A_{\Gamma}^{\text{NNLL}}(\nu,\mu)+\cdots\,,$
$\displaystyle
A_{\gamma_{j}}(\nu,\mu)=A_{\gamma_{j}}^{\text{LL}}(\nu,\mu)+A_{\gamma_{j}}^{\text{NLL}}(\nu,\mu)+A_{\gamma_{j}}^{\text{NNLL}}(\nu,\mu)+\cdots\,.$
(C.8)

The corresponding term at each order is

$\displaystyle A_{\Gamma}^{\text{LL}}(\nu,\mu)$
$\displaystyle=\frac{\Gamma_{0}}{2\beta_{0}}\ln r\,,$ $\displaystyle
A_{\Gamma}^{\text{NLL}}(\nu,\mu)$
$\displaystyle=\frac{\Gamma_{0}\alpha_{s}(\nu)}{8\pi\beta_{0}}\left(\frac{\beta_{1}}{\beta_{0}}-\frac{\Gamma_{1}}{\Gamma_{0}}\right)(1-r)\,,$
$\displaystyle A_{\Gamma}^{\text{NNLL}}(\nu,\mu)$
$\displaystyle=\frac{\Gamma_{0}\alpha(\nu)^{2}}{64\pi^{2}\beta_{0}}\left(\frac{\beta_{1}^{2}}{\beta_{0}^{2}}-\frac{\beta_{2}}{\beta_{0}}+\frac{\Gamma_{2}}{\Gamma_{0}}-\frac{\beta_{1}\Gamma_{1}}{\beta_{0}\Gamma_{0}}\right)\left(r^{2}-1\right)\,.$
(C.9)

##  Appendix D Dictionary between position and momentum space logs

In this section, we provide the analytic Fourier transformation rules that can
be used to relate singular cross sections and soft and jet functions in
position ($z$) and momentum ($r$) space, which is required e.g. to extract the
correct $\mathcal{O}(\alpha_{s})$ soft constant term in $z$ space. These rules
are distinct from the commonly considered case of plus distributions with
support only on the positive reals, and further require care due to the
additional double integral over the second derivative
$\frac{d^{3}\sigma}{d\rho^{3}}$ of the spectrum. For example, the singular
cross section in position space at a fixed scale $\mu=Q$ is given by

$\displaystyle\tilde{\sigma}(z,\mu_{h,j,s}=Q)$
$\displaystyle=3+\Bigl{(}\frac{\alpha_{s}}{4\pi}\Bigr{)}\biggl{\\{}-6(C_{A}+2C_{F})\ln^{2}\left(|z|e^{\gamma_{E}}\right)$
(D.1) $\displaystyle\quad+\Bigl{[}-4n_{f}T_{F}+C_{A}(11+12\ln 3)+6C_{F}(3+4\ln
3)$
$\displaystyle\quad\quad-2(C_{A}+2C_{F})\,i\pi\,\text{sgn}(z)\Bigl{]}\ln\left(|z|e^{\gamma_{E}}\right)$
$\displaystyle\quad+\frac{1}{6}\Bigl{(}-4n_{f}T_{F}+C_{A}(11+12\ln
48)+6C_{F}(3+4\ln 48)\Bigl{)}\,i\pi\,\text{sgn}(z)+\text{const}\biggr{\\}}.$

To obtain the singular expression in momentum ($r$) space, we need to take the
inverse Fourier transform Eq. (D) and integrate over $r$ twice.

To achieve this, we derive the map between logs in the momentum $z$ space and
position $r$ space. Specifically, we consider the dictionary between plus
distributions of the form

$\displaystyle\delta(r),\quad\left(\frac{\theta(r)\ln^{m}r}{r}\right)_{+},\quad\left(\frac{\theta(-r)\ln^{m}(-r)}{-r}\right)_{+}\,,\quad
m=0,1,2,\dots\,,$ (D.2)

which appear in $\frac{d^{3}\sigma}{d\rho^{3}}$, and Fourier-space logarithms

$\displaystyle\ln^{m}(|z|e^{\gamma_{E}}),\quad\ln^{m}(|z|e^{\gamma_{E}})\,\text{sgn}(z)\,,\quad
m=0,1,2,\dots$ (D.3)

appearing in $\tilde{\sigma}(z)$. When integrating over $r$ twice in addition,
which is accounted for by the kernel $K(z,r)$ in Eq.(4.13), we also obtain the
transformation rules between position logs and the spectrum logs:

$\displaystyle\theta(r)r\ln^{m}r,\quad\theta(-r)(-r)\ln^{m}(-r)\,,\quad
m=0,1,2,\dots\,.$ (D.4)

To derive the dictionary for the first part of the basis in Eq. (D.3), we
consider the generating functional

$f_{1}(z,\epsilon)=(-ize^{\gamma_{E}})^{\epsilon}(ize^{\gamma_{E}})^{\epsilon}\,,$
(D.5)

and its inverse Fourier transformation

$\mathcal{F}^{-1}[f_{1}(z,\epsilon)]=-\frac{e^{2\epsilon\gamma_{E}}\Gamma(1+2\epsilon)\sin(\pi\epsilon)}{\pi}|r|^{-1-2\epsilon}\,,$
(D.6)

where $|r|^{-1-2\epsilon}$ is expanded in terms of plus distributions as

$|r|^{-1-2\epsilon}=-\frac{\delta\left(|r|\right)}{2\epsilon}+\left(\frac{1}{|r|}\right)_{+}-2\epsilon\left(\frac{\ln|r|}{|r|}\right)_{+}+\cdots\,.$
(D.7)

We remind the reader that eqn (D.6) is simply a special case of eqn (3.2) with
$a=b=-\epsilon$. The distribution $\delta\left(|r|\right)$ needs to be defined
carefully, and we choose the following definition

$\delta(|r|)\equiv 2\frac{d^{2}}{dr^{2}}\left[\xi
r\theta(r)+r(\xi-1)\theta(-r)\right]\overset{\xi=1}{=}2\delta(r)\,.$ (D.8)

As explained in Sec. 4.1, we choose the LO condition $\xi=1$. Then if we
integrate $f_{2}(r,\epsilon)$ over $r$ twice subject to the $\xi=1$ condition,
which is equivalent to using the integral kernel Eq. (4.13) (let us call this
combined operation $\mathcal{I}$), we find

$\mathcal{I}[f_{1}(z,\epsilon)]=-\frac{e^{2\epsilon\gamma_{E}}\Gamma(-1+2\epsilon)\sin(\pi\epsilon)}{\pi}\left(r+|r|^{1-2\epsilon}\right)\,.$
(D.9)

Expanding $f_{1}(z,\epsilon)$ and $\mathcal{F}^{-1}[f_{1}(z,\epsilon)]$ in
$\epsilon$, we find the inverse Fourier transformation rules to any power of
$\ln^{m}\left(|z|e^{\gamma_{E}}\right)$. The first few orders read,

$\displaystyle\mathcal{F}^{-1}\left[1\right]$
$\displaystyle=\frac{1}{2}\delta\left(|r|\right)\,,$
$\displaystyle\mathcal{F}^{-1}\left[\ln(|z|e^{\gamma_{E}})\right]$
$\displaystyle=-\frac{1}{2}\left(\frac{1}{|r|}\right)_{+}=-\frac{1}{2}\left[\left(\frac{\theta(r)}{r}\right)_{+}+\left(\frac{\theta(-r)}{-r}\right)_{+}\right]\,,$
$\displaystyle\mathcal{F}^{-1}\left[\ln^{2}(|z|e^{\gamma_{E}})\right]$
$\displaystyle=\left(\frac{\ln|r|}{|r|}\right)_{+}+\frac{\pi^{2}}{24}\delta\left(|r|\right)=\left[\left(\frac{\theta(r)\ln
r}{r}\right)_{+}+\left(\frac{\theta(-r)\ln(-r)}{-r}\right)_{+}\right]+\frac{\pi^{2}}{24}\delta\left(|r|\right)\,,$
$\displaystyle\mathcal{F}^{-1}\left[\ln^{3}(|z|e^{\gamma_{E}})\right]$
$\displaystyle=-\frac{\pi^{2}}{8}\left(\frac{1}{|r|}\right)_{+}-\frac{3}{2}\left(\frac{\ln^{2}|r|}{|r|}\right)_{+}+\frac{1}{2}\psi^{(2)}(1)\delta\left(|r|\right)\,,$
$\displaystyle=-\frac{\pi^{2}}{8}\left[\left(\frac{\theta(r)}{r}\right)_{+}+\left(\frac{\theta(-r)}{-r}\right)_{+}\right]-\frac{3}{2}\left[\left(\frac{\theta(r)\ln^{2}r}{r}\right)_{+}+\left(\frac{\theta(-r)\ln^{2}(-r)}{-r}\right)_{+}\right]$
$\displaystyle\quad+\frac{1}{2}\psi^{(2)}(1)\delta\left(|r|\right)\,.$ (D.10)

where the digamma function is defined as
$\psi(z)\equiv\frac{\Gamma^{\prime}(z)}{\Gamma(z)}$ and
$\psi^{(n)}(z)\equiv\frac{d^{n}\psi(z)}{dz^{n}}$. From the RHS, we can see
that the first log basis $\ln^{m}(|z|e^{\gamma_{E}})$ is mapped to the sum of
both left and right shoulder r space logs, which is symmetric. Furthermore, we
can also expand $\mathcal{I}\left[f_{1}(z,\epsilon)\right]$ in $\epsilon$ and
read off the rules between $\ln^{m}(|z|e^{\gamma_{E}})$ and spectrum logs

$\displaystyle\mathcal{I}[1]$
$\displaystyle=\frac{1}{2}|r|\left(1+\text{sgn}(r)\right)=r\theta(r)\,,$
(D.11) $\displaystyle\mathcal{I}\left[\ln(|z|e^{\gamma_{E}})\right]$
$\displaystyle=\frac{1}{2}|r|\left(1+\text{sgn}(r)\right)-\frac{1}{2}|r|\ln|r|=\theta(r)\left(r-\frac{1}{2}r\ln
r\right)+\theta(-r)\frac{1}{2}r\ln(-r)\,,$
$\displaystyle\mathcal{I}\left[\ln^{2}(|z|e^{\gamma_{E}})\right]$
$\displaystyle=\left(1+\frac{\pi^{2}}{24}\right)|r|\left(1+\text{sgn}(r)\right)+\frac{1}{2}|r|\ln^{2}|r|-|r|\ln|r|$
$\displaystyle=\theta(r)\left(2r+\frac{\pi^{2}}{12}r-r\ln
r+\frac{1}{2}r\ln^{2}r\right)+\theta(-r)\left(r\ln(-r)-\frac{1}{2}r\ln^{2}(-r)\right)\,.$

For the second group of terms in the basis,
$\ln^{m}(|z|e^{\gamma_{E}})\,\text{sgn}(z)$, we consider the generating
functional

$f_{2}(z,\epsilon)=i\pi
e^{2\epsilon\gamma_{E}}|z|^{2\epsilon}\text{sgn}(z)\,,$ (D.12)

and its inverse Fourier transformation