# Joint Two-Dimensional Resummation in $q_{T}$ and $0$-Jettiness at NNLL

Gillian Lustermans  Johannes K. L. Michel  Frank J. Tackmann
and Wouter J. Waalewijn

(January 10, 2019)

###### Abstract

We consider Drell-Yan production $pp\to Z/\gamma^{*}\to\ell^{+}\ell^{-}$ with
the simultaneous measurement of the $Z$-boson transverse momentum $q_{T}$ and
$0$-jettiness $\mathcal{T}_{0}$. Since both observables resolve the initial-
state QCD radiation, the double-differential cross section in $q_{T}$ and
$\mathcal{T}_{0}$ contains Sudakov double logarithms of both $q_{T}/Q$ and
$\mathcal{T}_{0}/Q$, where $Q\sim m_{Z}$ is the dilepton invariant mass. We
simultaneously resum the logarithms in $q_{T}$ and $\mathcal{T}_{0}$ to next-
to-next-to-leading logarithmic order (NNLL) matched to next-to-leading fixed
order (NLO). Our results provide the first genuinely two-dimensional analytic
Sudakov resummation for initial-state radiation. Integrating the resummed
double-differential spectrum with an appropriate scale choice over either
$\mathcal{T}_{0}$ or $q_{T}$ recovers the corresponding single-differential
resummation for the remaining variable. We discuss in detail the required
effective field theory setups and their combination using two-dimensional
resummation profile scales. We also introduce a new method to perform the
$q_{T}$ resummation where the underlying resummation is carried out in impact-
parameter space, but is consistently turned off depending on the momentum-
space target value for $q_{T}$. Our methods apply at any order and for any
color-singlet production process, such that our results can be systematically
extended when the relevant perturbative ingredients become available.

##  1 Introduction

The increasing accuracy of measurements at the LHC places high demands on the
precision and versatility of theoretical predictions. Fixed-order perturbation
theory has proven to be a powerful tool to describe a large number of LHC
processes, provided the measurement is sufficiently inclusive. With increasing
data sets, however, more fine-grained measurements become possible and
increasingly differential quantities come into focus. These more exclusive
cross sections often involve several physical scales set by the hard
interaction and the differential measurements or cuts applied on the final
state. When these scales are widely separated, the perturbative series at each
order is dominated by logarithms of their ratios. The resummation of these
logarithms to all orders is crucial to arrive at the best possible
predictions.

The resummation for measurements sensitive to infrared (soft and/or collinear)
physics can, in part, be achieved through the use of parton-shower Monte Carlo
event generators; popular examples include Pythia [1, 2], Herwig [3, 4], or
Sherpa [5]. Parton showers provide fully exclusive final states so that in
principle, any desired measurements or cuts can be imposed on the generated
events. Existing implementations of parton showers are only formally accurate
at about leading-logarithmic (LL) level, depending on the shower’s evolution
variable (and other implementation details) and the observable in question. (A
recent detailed analysis can be found in ref. [6].) Furthermore, estimating
the perturbative uncertainties of parton showers is challenging, which is in
part due to their limited perturbative accuracy.

Analytic methods for the higher-order resummation of infrared-sensitive
observables are available. These include the CSS formalism [7, 8, 9],
seminumerical methods based on the coherent-branching formalism [10, 11, 12,
13], and methods using renormalization group evolution (RGE) in effective
field theories (EFTs) of QCD such as soft-collinear effective theory (SCET)
[14, 15, 16, 17, 18, 19]. The common drawback of analytic resummation methods
is that they only apply after a sufficient amount of emissions have been
integrated over, which is why they have been primarily used for the
resummation of single-differential observables. Their crucial advantage is
that they can be systematically extended to higher orders, and theoretical
uncertainties can be addressed in a more reliable way.

There has been much progress in extending analytic resummation methods to
cases involving multiple resummation variables. Examples include the joint
resummation of transverse momentum $q_{T}$ and threshold (large $x$)
logarithms [20, 21, 22, 23, 24, 25, 26], $q_{T}$ and small $x$ [27],
$N$-jettiness (or jet mass) together with dijet invariant masses [28, 29], two
angularities [30, 31], jet mass and jet radius [32], jet vetoes and jet
rapidity [33, 34], or threshold and jet radius in inclusive jet production
[35, 36]. Most of these examples either involve different variables that
effectively resolve different subsequent emissions, or involve a primary
resummation variable that is modified by an auxiliary measurement or
constraint. Another well-understood case is when an infrared-sensitive
measurement is separated into its contributions from mutually exclusive
regions of phase space [37, 38, 39].111Yet another case, which will not be
relevant here, arises when different infrared-sensitive measurements are
performed in different regions of phase space, which may require the
resummation of nonglobal logarithms [40, 41, 42, 43, 44, 45].

In contrast, here we are interested in resolving emissions at the same level
by simultaneously measuring two independent infrared-sensitive observables.
Extending analytic resummation to such genuinely multi-dimensional resolution
variables is of key theoretical concern, as it allows for a more complete
description of the emission pattern beyond LL, effectively filling a gap
between analytic resummations and parton showers. So far, this has been
achieved at NNLL for the case of simultaneously measuring two angularities in
$e^{+}e^{-}$ collisions [31].

In this paper, we consider Drell-Yan, $pp\to
Z/\gamma^{\ast}\to\ell^{+}\ell^{-}$, with a simultaneous measurement of (1)
the transverse momentum $q_{T}$ of the Drell-Yan lepton pair and (2) the
hadronic resolution variable 0-jettiness
${\mathcal{T}}\equiv{\mathcal{T}}_{0}$ [46, 47]. Achieving their combined
resummation is important conceptually because $q_{T}$ and ${\mathcal{T}}$ are
prototypes for two large classes of infrared-sensitive observables: $q_{T}$
constrains the transverse momentum of initial-state radiation, while
${\mathcal{T}}$ constrains its virtuality. These different behaviors lead to
very different logarithmic structures already at LL, which in SCET is
reflected in the RGE structure of two distinct effective theories, SCETI and
SCETII. (For parton showers, these correspond to evolution variables based on
either transverse momentum or virtuality, respectively.)

Beyond providing a prototype for combining SCETI and SCETII resummations, the
joint resummation of $q_{T}$ and ${\mathcal{T}}$ is also of direct
phenomenological interest. First, they are important variables individually.
The measurement of ${\mathcal{T}}$ in bins of $q_{T}$ [48] can probe the so-
called underlying event in hadronic collisions. Furthermore, the Geneva Monte
Carlo event generator [49, 50] uses ${\mathcal{T}}$ as the underlying
resolution variable for the event generation, achieving NNLL${}^{\prime}+$NNLO
accuracy in ${\mathcal{T}}$ in conjunction with fully showered and hadronized
events. While other observables, such as $q_{T}$, benefit from the underlying
high resummation order, they do not enjoy the same level of formal accuracy in
Geneva as ${\mathcal{T}}$ itself. The joint resummation of ${\mathcal{T}}$ and
$q_{T}$ to a given order enables extending the event generation in Geneva to
also be accurate in $q_{T}$ to the same order.

Figure -1019: The Drell-Yan cross section double-differential in the
transverse momentum $q_{T}$ of the $Z$ boson and the 0-jettiness event shape
${\mathcal{T}}$ at NNLL$+$NLO. For better visibility, the spectrum is plotted
with respect to $\log_{10}q_{T}$ and $\log_{10}{\mathcal{T}}$. On the two side
walls we show the corresponding single-differential spectra in $q_{T}$ and
${\mathcal{T}}$ obtained by integrating the double-differential spectrum up to
${\mathcal{T}_{\mathrm{cut}}}=100\,\mathrm{GeV}$ and
${q_{T}^{\mathrm{cut}}}=100\,\mathrm{GeV}$, respectively.

The double-differential factorization for $q_{T}$ and ${\mathcal{T}}$ was
first considered in ref. [51]. There, the regions of phase space where $q_{T}$
(SCETII) and ${\mathcal{T}}$ (SCETI) determine the resummation structure were
identified, together with the appropriate intermediate effective theory SCET+
[28, 51] that connects them. Here, we develop an explicit matching procedure
that combines the three different theories, SCETI, SCET+, and SCETII, such
that the resummation structure of each is recovered in its respective region
of phase space. In particular, our method ensures that the single-differential
resummation in one variable is recovered upon integration over the other. We
discuss in detail the technical challenges involved. These include the
construction of appropriate two-dimensional profile scales to combine the
SCETII resummation for $q_{T}$, which is performed in position (impact-
parameter) space, with the SCETI resummation for ${\mathcal{T}}$, which is
performed in momentum space, the estimation of perturbative uncertainties, and
the matching to full QCD at large $q_{T}$ and/or ${\mathcal{T}}$ in a flexible
way and consistent with the corresponding single-differential cases. We obtain
explicit numerical predictions for the double-differential
$(q_{T},{\mathcal{T}})$ spectrum, achieving its complete and fully two-
dimensional Sudakov resummation at NNLL$+$NLO. Our main result is shown in
figure -1019, featuring a nice two-dimensional Sudakov peak structure.

We like to stress that our methods are completely general and can be applied
to any color-singlet production process and at any order for which the
relevant perturbative ingredients are available. (Some of the double-
differential ingredients required at NNLL′ and N3LL are already known [52].)
Furthermore, our matching procedure is generic and can be applied to any type
of two-dimensional resummation for which the relevant EFTs on the boundaries
and in the bulk are known.

The remainder of the paper is organized as follows. In section 2, we discuss
the three different parametric regimes and the factorization and resummation
for each individually. In section 3, we then discuss in detail our method for
consistently combining them to obtain a complete description of the two-
dimensional $(q_{T},{\mathcal{T}})$ plane. Our numerical results for the
double-differential spectrum at NNLL$+$NLO are presented in section 4. We
conclude in section 5. In appendix A we summarize our conventions for plus
distributions and Fourier transforms. All required perturbative ingredients
are collected in appendix B.

##  2 Resummation framework

###  2.1 Overview of parametric regimes

Figure -1018: Parametric regimes in the $(q_{T},{\mathcal{T}})$ plane and
their SCET description. The solid lines correspond to the phase-space
boundaries $q_{T}={\mathcal{T}}$ (green) and $q_{T}=\sqrt{Q{\mathcal{T}}}$
(blue).

We consider color-singlet production at hadron colliders. Although the process
dependence is not important for our discussion, we consider the example of
Drell-Yan production, $pp\to Z/\gamma^{\ast}(\to\ell^{+}\ell^{-})$, for
concreteness. We measure the total invariant mass $Q$ and rapidity $Y$ of the
color-singlet final state (the lepton pair). The two resolution variables we
measure are the transverse momentum $q_{T}$ of the color-singlet final state
and the 0-jettiness ${\mathcal{T}}$ (aka beam thrust) [46, 47, 53, 38],
defined as

$\displaystyle{\mathcal{T}}\equiv{\mathcal{T}}_{0}=\sum_{i}\,\min\Bigl{\\{}\frac{2q_{a}\cdot
k_{i}}{Q_{a}},\frac{2q_{b}\cdot k_{i}}{Q_{b}}\Bigr{\\}}\,.$ (2.1)

The sum runs over all particles $i$ with momentum $k_{i}$ in the final state,
excluding the color-singlet final state. We choose the massless reference
momenta $q_{a}$ and $q_{b}$ as

$\displaystyle q_{a}^{\mu}=\frac{Qe^{+Y}}{2}\,n_{a}^{\mu}\,,\qquad
q_{b}^{\mu}=\frac{Qe^{-Y}}{2}\,n_{b}^{\mu}\,,$ (2.2)

where $n_{a}^{\mu}=(1,+\hat{z})$ and $n_{b}^{\mu}=(1,-\hat{z})$ are lightlike
vectors along the beam axis $\hat{z}$. For definiteness we use the leptonic
definition of 0-jettiness for all numerical results in this paper, for which
the measure factors $Q_{a,b}$ are simply given by

$Q_{a}=Q_{b}=Q\,.$ (2.3)

Our setup applies equally well to other definitions of ${\mathcal{T}}$, so we
keep $Q_{a}$ and $Q_{b}$ (with $Q_{a}Q_{b}=Q^{2}$) generic for the rest of
this section.

We are interested in the contribution of initial-state radiation (ISR) to the
simultaneous measurement of $q_{T},{\mathcal{T}}\ll Q$, where
$Q\gg\Lambda_{\mathrm{QCD}}$ sets the scale of the hard interaction. The
dynamics of perturbative ISR is then governed by three distinct momentum
scales set by the measurement of $q_{T}$ and ${\mathcal{T}}$. First, the
typical transverse momentum of emissions that recoil against the lepton pair
is set by $q_{T}$. Second, isotropic (soft) emissions at central rapidities
can contribute to ${\mathcal{T}}$ via either of the projections onto
$q_{a}^{\mu}$ and $q_{b}^{\mu}$ in eq. (2.1). This implies that their
characteristic transverse momentum is $\sim{\mathcal{T}}$. Third, ISR with
typical energy $\sim Q$ can contribute to ${\mathcal{T}}$ as long as it is
collinear to either of the incoming beams, such that its contribution to
${\mathcal{T}}$ in eq. (2.1) is small. These collinear emissions then have a
typical transverse momentum $\sim\sqrt{Q{\mathcal{T}}}$. The factorization and
resummation structure of the cross section for $q_{T},{\mathcal{T}}\ll Q$
depends on the parametric hierarchy between these scales. There are three
relevant parametric regimes [51], which are illustrated in figure -1018 and
are discussed in the following.

In the first (blue) regime, ${\mathcal{T}}\ll q_{T}\sim\sqrt{Q{\mathcal{T}}}$,
soft emissions with transverse momentum $\sim{\mathcal{T}}$ and collinear
emissions with transverse momentum $\sim\sqrt{Q{\mathcal{T}}}$ both contribute
to the ${\mathcal{T}}$ measurement. Due to the separation in transverse
momentum, the $q_{T}$ measurement is determined by collinear emissions, while
soft emissions do not contribute to it. The appropriate EFT description for
this regime is SCETI. It has the same RG structure as the single-differential
${\mathcal{T}}$ spectrum, with $q_{T}$ acting as an auxiliary variable. The
SCETI regime is discussed in more detail in section 2.2.

In the opposite (green) regime, ${\mathcal{T}}\sim
q_{T}\ll\sqrt{Q{\mathcal{T}}}$, both soft and collinear emissions have
transverse momentum $\sim q_{T}$ and thus contribute to $q_{T}$. On the other
hand, only soft radiation at central rapidities contributes to
${\mathcal{T}}$, while the contribution from collinear radiation is
suppressed. This regime is described by SCETII, whose RG structure is
analogous to that of the single-differential $q_{T}$ spectrum, with
${\mathcal{T}}$ as the auxiliary variable. The SCETII regime is discussed in
more detail in section 2.3.

Third, the intermediate (orange) regime in the bulk, ${\mathcal{T}}\ll
q_{T}\ll\sqrt{Q{\mathcal{T}}}$, shares features with both boundary cases. As
in the SCETI regime, central soft radiation contributes to ${\mathcal{T}}$,
while as in the SCETII regime, collinear radiation contributes to $q_{T}$. In
addition, this regime requires a distinct collinear-soft mode at an
intermediate rapidity scale that can contribute to both measurements [51]. The
relevant EFT description is provided by SCET+, which in this case shares
elements of both SCETI and SCETII. The SCET+ regime, as well as its relation
to the regimes on the two boundaries, is discussed in section 2.4. We briefly
comment on the regions beyond the phase-space boundaries (left blank in figure
-1018) in section 2.5.

All numerical results for the SCET predictions in the following are obtained
from our implementation in SCETlib [54]. All fixed NLO results in full QCD are
obtained from MCFM 8.0 [55, 56, 57]. Throughout this paper we use
MMHT2014nnlo68cl [58] NNLO PDFs with $\alpha_{s}(m_{Z})=0.118$ and five active
quark flavors.

###  2.2 SCETI: ${\mathcal{T}}\ll q_{T}\sim\sqrt{Q{\mathcal{T}}}$

In this regime, both soft and collinear modes are constrained by
${\mathcal{T}}$, while only collinear modes can contribute to $q_{T}$, whose
characteristic transverse momentum $\sqrt{Q{\mathcal{T}}}$ coincides
parametrically with $q_{T}$. The scaling of the relevant EFT modes reads

$\displaystyle n_{a}\text{-collinear:}$ $\displaystyle\quad
p^{\mu}\sim\Bigl{(}{\mathcal{T}},Q,\sqrt{Q{\mathcal{T}}}\Bigr{)}\sim\Bigl{(}\tfrac{q_{T}^{2}}{Q},Q,q_{T}\Bigr{)}\,,$
$\displaystyle n_{b}\text{-collinear:}$ $\displaystyle\quad
p^{\mu}\sim\Bigl{(}Q,{\mathcal{T}},\sqrt{Q{\mathcal{T}}}\Bigr{)}\sim\Bigl{(}Q,\tfrac{q_{T}^{2}}{Q},q_{T}\Bigr{)}\,,$
soft: $\displaystyle\quad
p^{\mu}\sim\Bigl{(}{\mathcal{T}},{\mathcal{T}},{\mathcal{T}}\Bigr{)}\,,$ (2.4)

in terms of lightcone coordinates defined by (with $n\equiv n_{a}$,
${\bar{n}}\equiv n_{b}$)

$p^{\mu}=n\\!\cdot\\!p\,\frac{{\bar{n}}^{\mu}}{2}+{\bar{n}}\\!\cdot\\!p\,\frac{n^{\mu}}{2}+p_{\perp}^{\mu}\equiv(n\\!\cdot\\!p,{\bar{n}}\\!\cdot\\!p,p_{\perp})\equiv(p^{+},p^{-},p_{\perp})\,.$
(2.5)

This leads to the following factorization formula for the cross section [46,
59],

$\displaystyle\frac{\mathrm{d}\sigma_{\mathrm{I}}}{\mathrm{d}Q\,\mathrm{d}Y\,\mathrm{d}q_{T}\,\mathrm{d}{\mathcal{T}}}$
$\displaystyle=H_{\kappa}(Q,\mu)\int\\!\mathrm{d}t_{a}\int\\!\mathrm{d}^{2}\vec{\kern
0.0ptk}_{a}\,B_{a}(t_{a},x_{a},\vec{\kern
0.0ptk}_{a},\mu)\int\\!\mathrm{d}t_{b}\int\\!\mathrm{d}^{2}\vec{\kern
0.0ptk}_{b}\,B_{b}(t_{b},x_{b},\vec{\kern 0.0ptk}_{b},\mu)$
$\displaystyle\quad\times\int\\!\mathrm{d}k\,S_{\kappa}(k,\mu)\,\delta\bigl{(}q_{T}-\lvert\vec{\kern
0.0ptk}_{a}+\vec{\kern
0.0ptk}_{b}\rvert\bigr{)}\,\delta\Bigl{(}{\mathcal{T}}-\frac{t_{a}}{Q_{a}}-\frac{t_{b}}{Q_{b}}-k\Bigr{)}\,,$
(2.6)

which holds up to power corrections of the form222Lorentz invariance suggests
that power corrections in $q_{T}$ always appear in terms of $q_{T}^{2}$. This
distinction is irrelevant for our discussion.

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}Q\,\mathrm{d}Y\,\mathrm{d}q_{T}\,\mathrm{d}{\mathcal{T}}}$
$\displaystyle=\frac{\mathrm{d}\sigma_{\mathrm{I}}}{\mathrm{d}Q\,\mathrm{d}Y\,\mathrm{d}q_{T}\,\mathrm{d}{\mathcal{T}}}\,\Big{[}1+\mathcal{O}\Bigl{(}\frac{{\mathcal{T}}}{Q},\frac{q_{T}^{2}}{Q^{2}},\frac{{\mathcal{T}}^{2}}{q_{T}^{2}}\Bigr{)}\Big{]}\,.$
(2.7)

The hard function $H_{\kappa}(Q,\mu)$ describes the short-distance scattering
that produces the lepton pair through the off-shell $\gamma^{\ast}$ or $Z$. In
addition to $Q$, it depends on the partonic channel $\kappa\equiv\\{a,b\\}$,
which is implicitly summed over all relevant combinations of quark and
antiquark flavors $a,b$ on the right-hand side of eq. (2.2). The beam
functions $B_{q}(t,x,\vec{\kern 0.0ptk}_{T},\mu)$ describe extracting a quark
(or antiquark) from the proton with momentum fraction $x$, virtuality $t$, and
transverse momentum $\vec{\kern 0.0ptk}_{T}$. The momentum fractions are
directly related to $Q$ and $Y$,

$x_{a}=\frac{Q}{E_{\mathrm{cm}}}\,e^{+Y}\,,\qquad
x_{b}=\frac{Q}{E_{\mathrm{cm}}}\,e^{-Y}\,.$ (2.8)

The $t$ and $\vec{\kern 0.0ptk}_{T}$ encode the contribution of the collinear
radiation to the ${\mathcal{T}}$ and $q_{T}$ measurement, as captured by the
measurement $\delta$ functions on the last line of eq. (2.2). For $t\sim
k_{T}^{2}\gg\Lambda_{\mathrm{QCD}}^{2}$, these beam functions can be matched
onto PDFs [46, 60, 59],

$\displaystyle B_{q}(t,x,\vec{\kern 0.0ptk}_{T},\mu)$
$\displaystyle=\sum_{j}\int_{x}^{1}\frac{\mathrm{d}z}{z}\,\mathcal{I}_{qj}(t,z,\vec{\kern
0.0ptk}_{T},\mu)\,f_{j}\Bigl{(}\frac{x}{z},\mu\Bigr{)}\bigg{[}1+\mathcal{O}\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}^{2}}{t},\frac{\Lambda_{\mathrm{QCD}}^{2}}{k_{T}^{2}}\Bigr{)}\bigg{]}\,.$
(2.9)

The soft function $S_{\kappa}(k,\mu)$ encodes the contribution from soft
radiation to the 0-jettiness measurement, and depends on the color charge of
the colliding partons.

The factorization in eq. (2.2) separates the physics at the canonical SCETI
scales

$\mu_{H}^{\mathrm{I}}\sim
Q\,,\qquad\mu_{B}^{\mathrm{I}}\sim\sqrt{Q{\mathcal{T}}}\,,\qquad\mu_{S}^{\mathrm{I}}\sim{\mathcal{T}}\,.$
(2.10)

By evaluating the ingredients at their natural scale and evolving them to a
common scale, all logarithms of
${\mathcal{T}}/Q\sim\mu^{\mathrm{I}}_{S}/\mu^{\mathrm{I}}_{H}\sim(\mu^{\mathrm{I}}_{B}/\mu^{\mathrm{I}}_{H})^{2}\sim(\mu^{\mathrm{I}}_{S}/\mu^{\mathrm{I}}_{B})^{2}$
are resummed.

The hard and soft function in eq. (2.2) are the same as in the single-
differential ${\mathcal{T}}$ spectrum and do not depend on $q_{T}$. The RG
consistency of the cross section then implies that the RGE of the double-
differential beam functions cannot depend on $q_{T}$, such that the overall RG
structure of the cross section is equivalent to the single-differential case,
i.e., $q_{T}$ takes the role of an auxiliary measurement in the SCETI
resummation, with no large logarithms of $q_{T}$ appearing in the cross
section as long as $q_{T}\sim\sqrt{Q{\mathcal{T}}}$ is satisfied. We stress
that eq. (2.2) nevertheless provides a nontrivial and genuinely double-
differential extension of the single-differential case. This is already
visible from the structure of power corrections in eq. (2.7). Furthermore, the
$q_{T}$ dependence does affect and is affected by the ${\mathcal{T}}$
resummation because the double-differential beam functions enter in a
convolution with the beam and soft renormalization group kernels. Physically,
they account for the total $q_{T}$ recoil from all collinear emissions that
are being resummed in ${\mathcal{T}}$.

The factorization of the double-differential spectrum in eq. (2.2) (and in the
following sections) does not account for effects from Glauber gluon exchange.
For active-parton scattering, they are expected to enter at
$\mathcal{O}(\alpha_{s}^{4})$ (N4LL′) [61, 62], which is well beyond the order
we are interested in. They can be included using the Glauber operator
framework of ref. [63]. For proton initial states the factorization formula
also does not account for spectator forward scattering effects. Their complete
treatment for the single-differential ${\mathcal{T}}$ spectrum is not yet
available, but we expect that their treatment for the double-differential case
would follow in a similar way.

#### Scale setting and fixed-order matching.

Figure -1017: Left: Comparison of singular and nonsingular contributions to
the fixed $\mathcal{O}(\alpha_{s})$ double spectrum as a function of
${\mathcal{T}}$, with $q_{T}=\sqrt{Q{\mathcal{T}}}$ kept fixed. The orange
solid line shows the full QCD result and the dashed blue line the singular
contributions contained in the SCETI result eq. (2.2). The dotted green line
shows their difference, which corresponds to the power corrections indicated
in eq. (2.7). Right: SCETI profile scales and their associated variations. The
dotted lines (and the yellow band) indicate common up/down variations of
$\mu_{S}^{\mathrm{I}}$ and $\mu_{B}^{\mathrm{I}}$ from varying $\alpha$. The
dashed lines (and the green band) are variations of $\beta$ that only act on
$\mu_{B}^{\mathrm{I}}$. In both plots, the thin vertical lines correspond to
the transition points $(x_{0},x_{1},x_{2},x_{3})$ given in the text.

To extend the description of the cross section to large ${\mathcal{T}}\sim
q_{T}^{2}/Q\lesssim Q$, we have to reinstate the power corrections dropped in
eq. (2.7). This is achieved by matching to the full fixed-order result, for
which we use the standard additive matching,

$\mathrm{d}\sigma^{\mathrm{match}}_{\mathrm{I}}=\mathrm{d}\sigma_{\mathrm{I}}\bigr{|}_{\mu^{\mathrm{I}}}+\bigl{[}\mathrm{d}\sigma_{\mathrm{FO}}-\mathrm{d}\sigma_{\mathrm{I}}\bigr{]}_{\mu_{\mathrm{FO}}}\,.$
(2.11)

Here we abbreviated
$\mathrm{d}\sigma\equiv\mathrm{d}\sigma/(\mathrm{d}Q\,\mathrm{d}Y\,\mathrm{d}q_{T}\,\mathrm{d}{\mathcal{T}})$,
and $\mathrm{d}\sigma_{\mathrm{FO}}$ denotes the fixed-order cross section in
full QCD. The scale subscripts on the right-hand side indicate whether
$\mathrm{d}\sigma_{\mathrm{I}}$ is RG evolved using the SCETI resummation
scales $\mu^{\mathrm{I}}$, with their precise choices given below, or whether
it is evaluated with all scales set to a common fixed-order scale
$\mu_{\mathrm{FO}}$.

By construction, $\mathrm{d}\sigma_{\mathrm{I}}$ evaluated at common scales
$\mu_{\mathrm{FO}}$ exactly reproduces the singular limit of
$\mathrm{d}\sigma_{\mathrm{FO}}$, such that the term in square brackets in eq.
(2.11) is a pure nonsingular power correction at small ${\mathcal{T}}$, which
we can simply add to the resummed cross section. In the left panel of figure
-1017, we explicitly check that this is satisfied at fixed
$\mathcal{O}(\alpha_{s})$, and numerically assess the size of the power
corrections. We compare the full QCD result (solid orange) to the SCETI
singular limit (dashed blue) as a function of ${\mathcal{T}}$, while keeping
$q_{T}=\sqrt{Q{\mathcal{T}}}$ fixed to ensure that all classes of power
corrections in eq. (2.7) uniformly vanish as ${\mathcal{T}}\to 0$. This is
indeed satisfied, as the difference (dotted green) vanishes like a power.

For ${\mathcal{T}}\sim Q$, the SCETI singular contribution and the power
corrections are of the same size, implying that the resummation must be turned
off to not upset the $\mathcal{O}(1)$ cancellation between them and correctly
reproduce the fixed-order result. This is commonly achieved by using profile
scales [64, 65], i.e., by having
$\mu_{B}^{\mathrm{I}}\equiv\mu_{B}^{\mathrm{I}}({\mathcal{T}})$ and
$\mu_{S}^{\mathrm{I}}\equiv\mu_{S}^{\mathrm{I}}({\mathcal{T}})$ transition
from their canonical values eq. (2.10) at small ${\mathcal{T}}$ to a common
high scale for large ${\mathcal{T}}$, schematically,

$\mu^{\mathrm{I}}_{B}({\mathcal{T}})\,,\mu^{\mathrm{I}}_{S}({\mathcal{T}})\to\mu_{H}^{\mathrm{I}}=\mu_{\mathrm{FO}}\quad\text{for}\quad{\mathcal{T}}\to
Q\,.$ (2.12)

As a result, the first and third term in eq. (2.11) exactly cancel in this
limit, so the matched result reproduces $\mathrm{d}\sigma_{\mathrm{FO}}$ as
desired.

For the concrete choices of $\mu^{\mathrm{I}}_{B}$, $\mu^{\mathrm{I}}_{S}$ we
can rely on those used for the single-differential spectrum due to the
equivalent RG structure. We use the profile scale setup developed for the
closely related case of SCETI-like jet vetoes in ref. [66] and used for the
${\mathcal{T}}$ resummation in Geneva [50]. The profile scales are chosen as

$\displaystyle\mu_{S}^{\mathrm{I}}={\mu_{\mathrm{FO}}}\,f_{\mathrm{run}}^{\mathrm{I}}\Bigl{(}\frac{{\mathcal{T}}}{Q}\Bigr{)}\,,\qquad\mu_{B}^{\mathrm{I}}={\mu_{\mathrm{FO}}}\,\Bigl{[}f_{\mathrm{run}}^{\mathrm{I}}\Bigl{(}\frac{{\mathcal{T}}}{Q}\Bigr{)}\Bigr{]}^{1/2}\,,\qquad\mu_{H}^{\mathrm{I}}=\mu_{\mathrm{FO}}\,,$
(2.13)

with the profile function $f_{\mathrm{run}}^{\mathrm{I}}$ given by [67]

$\displaystyle f_{\mathrm{run}}^{\mathrm{I}}(x)$
$\displaystyle=\begin{cases}x_{0}\Bigl{(}1+\frac{x^{2}}{4x_{0}^{2}}\Bigr{)}&x\leq
2x_{0}\,,\\\ x&2x_{0}<x\leq x_{1}\,,\\\
x+\frac{(2-x_{2}-x_{3})(x-x_{1})^{2}}{2(x_{2}-x_{1})(x_{3}-x_{1})}&x_{1}<x\leq
x_{2}\,,\\\
1-\frac{(2-x_{1}-x_{2})(x-x_{3})^{2}}{2(x_{3}-x_{1})(x_{3}-x_{2})}&x_{2}<x\leq
x_{3}\,,\\\ 1&x_{3}<x\,.\end{cases}$ (2.14)

Based on figure -1017, we take $(x_{1},x_{2},x_{3})=(0.2,0.5,0.8)$ for the
transition points towards the fixed-order region $x\sim 1$. In addition, eq.
(2.14) turns off the resummation in the nonperturbative region $x\lesssim
2x_{0}$, where we set $x_{0}=1\,\mathrm{GeV}/Q$. This cuts off the
nonperturbative region and ensures that RG running induced by perturbative
anomalous dimensions always starts from a perturbative boundary condition. For
$\mu_{\mathrm{FO}}$ itself we use $\mu_{\mathrm{FO}}=Q$ as the central scale.
Our central scale choices are illustrated as solid lines in the right panel of
figure -1017.

#### Perturbative uncertainties.

We estimate perturbative uncertainties in
$\mathrm{d}\sigma_{\mathrm{I}}^{\mathrm{match}}$ by considering two different
sources. The first uncertainty contribution $\Delta_{\mathrm{I}}$ is inherent
to the SCETI resummation. It is estimated by varying the individual SCETI
scales while keeping ${\mu_{\mathrm{FO}}}$ fixed, effectively probing the
tower of higher-order logarithms that are being resummed. For this we use the
profile scale variations [66]

$\displaystyle\mu_{S}^{\mathrm{I}}$
$\displaystyle={\mu_{\mathrm{FO}}}\,\Bigl{[}f_{\mathrm{vary}}\Bigl{(}\frac{{\mathcal{T}}}{Q}\Bigr{)}\Bigr{]}^{\alpha}\,f_{\mathrm{run}}^{\mathrm{I}}\Bigl{(}\frac{{\mathcal{T}}}{Q}\Bigr{)}\,,$
$\displaystyle\mu_{B}^{\mathrm{I}}$
$\displaystyle={\mu_{\mathrm{FO}}}\,\Bigl{\\{}\Bigl{[}f_{\mathrm{vary}}\Bigl{(}\frac{{\mathcal{T}}}{Q}\Bigr{)}\Bigr{]}^{\alpha}\,f_{\mathrm{run}}^{\mathrm{I}}\Bigl{(}\frac{{\mathcal{T}}}{Q}\Bigr{)}\Bigr{\\}}^{1/2-\beta}\,,$
(2.15)

where $\alpha=\beta=0$ corresponds to the central scale choice in eq. (2.13),
and the variation factor is defined as

$\displaystyle f_{\mathrm{vary}}(x)=\begin{cases}2(1-x^{2}/x_{3}^{2})&0\leq
x<x_{3}/2\,,\\\ 1+2(1-x/x_{3})^{2}&x_{3}/2\leq x<x_{3}\,,\\\ 1&x_{3}\leq
x\,.\end{cases}$ (2.16)

It approaches a factor of two in the resummation region at small $x$ and
reduces to unity toward the fixed-order regime at $x=x_{3}$, where the
resummation is turned off. The estimate for $\Delta_{\mathrm{I}}$ is obtained
by computing $\mathrm{d}\sigma^{\mathrm{I}}_{\mathrm{match}}$ for each of the
four profile scale variations

$(\alpha,\beta)=\\{(+1,0),(-1,0),(0,+1/6),(0,-1/6)\\}\,,$ (2.17)

and taking the maximum absolute deviation from the central result. These
variations are also indicated in the right panel of figure -1017. Note that
for simplicity we do not perform explicit variations of the transition points
since they are known to have a subdominant effect, and the uncertainty in the
fixed-order matching is not essential to this paper.

For the second uncertainty contribution, $\Delta_{\mathrm{FO}}$, we consider
common variations of $\mu_{\mathrm{FO}}$ up and down by a factor of two in all
pieces of eq. (2.11). Since $\mu_{\mathrm{FO}}$ enters all $\mu^{\mathrm{I}}$
scales as a common overall factor, they inherit the same variation, which
keeps all resummed logarithms invariant. Hence, the $\mu_{\mathrm{FO}}$
variation effectively probes the effect of missing higher-order corrections in
the fixed-order contributions. The final uncertainty estimate for
$\mathrm{d}\sigma_{\mathrm{I}}^{\mathrm{match}}$ is obtained by adding both
contributions in quadrature,

$\Delta^{\mathrm{I}}_{\mathrm{total}}=\Delta_{\mathrm{I}}\oplus\Delta_{\mathrm{FO}}\equiv\bigl{(}\Delta_{\mathrm{I}}^{2}+\Delta_{\mathrm{FO}}^{2}\bigr{)}^{1/2}\,.$
(2.18)

The matched result $\mathrm{d}\sigma_{\mathrm{I}}^{\mathrm{match}}$ in eq.
(2.11) on its own constitutes a prediction for the double-differential
spectrum that covers the part of phase space where
$q_{T}\sim\sqrt{Q{\mathcal{T}}}$.

###  2.3 SCETII: ${\mathcal{T}}\sim q_{T}\ll\sqrt{Q{\mathcal{T}}}$

In this regime, both soft and collinear emissions are constrained by $q_{T}$.
Only soft radiation is constrained by the ${\mathcal{T}}$ measurement, while
collinear radiation at transverse momenta $\sim q_{T}\ll\sqrt{Q{\mathcal{T}}}$
is not affected by it. The relevant EFT modes scale as

$\displaystyle n_{a}\text{-collinear:}$ $\displaystyle\quad
p^{\mu}\sim\Bigl{(}\tfrac{q_{T}^{2}}{Q},Q,q_{T}\Bigr{)}\,,$ $\displaystyle
n_{b}\text{-collinear:}$ $\displaystyle\quad
p^{\mu}\sim\Bigl{(}Q,\tfrac{q_{T}^{2}}{Q},q_{T}\Bigr{)}\,,$ soft:
$\displaystyle\quad
p^{\mu}\sim\Bigl{(}q_{T},q_{T},q_{T}\Bigr{)}\sim\Bigl{(}{\mathcal{T}},{\mathcal{T}},{\mathcal{T}}\Bigr{)}\,.$
(2.19)

In this case, the cross section factorizes as [51]

$\displaystyle\frac{\mathrm{d}\sigma_{\mathrm{II}}}{\mathrm{d}Q\,\mathrm{d}Y\,\mathrm{d}q_{T}\,\mathrm{d}{\mathcal{T}}}$
$\displaystyle=H_{\kappa}(Q,\mu)\int\\!\mathrm{d}^{2}\vec{\kern
0.0ptk}_{a}\,B_{a}(\omega_{a},\vec{\kern
0.0ptk}_{a},\mu,\nu)\int\\!\mathrm{d}^{2}\vec{\kern
0.0ptk}_{b}\,B_{b}(\omega_{b},\vec{\kern 0.0ptk}_{b},\mu,\nu)$ (2.20)
$\displaystyle\quad\times\int\\!\mathrm{d}k\int\\!\mathrm{d}^{2}\vec{\kern
0.0ptk}_{s}\,S_{\kappa}(k,\vec{\kern
0.0ptk}_{s},\mu,\nu)\,\delta\bigl{(}q_{T}-\lvert\vec{\kern
0.0ptk}_{a}+\vec{\kern 0.0ptk}_{b}+\vec{\kern
0.0ptk}_{s}\rvert\bigr{)}\,\delta\bigl{(}{\mathcal{T}}-k\bigr{)}\,.$

The factorization receives power corrections of the form

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}Q\,\mathrm{d}Y\,\mathrm{d}q_{T}\,\mathrm{d}{\mathcal{T}}}$
$\displaystyle=\frac{\mathrm{d}\sigma_{\mathrm{II}}}{\mathrm{d}Q\,\mathrm{d}Y\,\mathrm{d}q_{T}\,\mathrm{d}{\mathcal{T}}}\,\Big{[}1+\mathcal{O}\Bigl{(}\frac{{\mathcal{T}}}{Q},\frac{q_{T}^{2}}{{\mathcal{T}}Q}\Bigr{)}\Big{]}\,.$
(2.21)

The hard function $H_{\kappa}(Q,\mu)$ is the same as in eq. (2.2). In SCETII
an additional regulator is required to handle rapidity divergences, for which
we use the $\eta$ regulator of refs. [68, 69] as implemented at two loops in
ref. [70], with $\nu$ the corresponding rapidity renormalization scale. The
$B_{q}(\omega,\vec{\kern 0.0ptk}_{T},\mu,\nu)$ are the same transverse
momentum-dependent beam functions as in the single-differential $q_{T}$
spectrum. The large momentum components $\omega$ in eq. (2.20) are given by

$\omega_{a}=Qe^{+Y}=x_{a}E_{\mathrm{cm}}\,,\qquad\omega_{b}=Qe^{-Y}=x_{b}E_{\mathrm{cm}}\,,$
(2.22)

and we suppress the trivial dependence of the beam function on
$E_{\mathrm{cm}}$. For $k_{T}^{2}\gg\Lambda_{\mathrm{QCD}}^{2}$, the beam
functions satisfy a matching relation similar to eq. (2.9) [71, 72, 73, 74,
69],

$\displaystyle B_{q}(\omega,\vec{\kern 0.0ptk}_{T},\mu,\nu)$
$\displaystyle=\sum_{j}\int_{x}^{1}\frac{\mathrm{d}z}{z}\,\mathcal{I}_{qj}(\omega,\vec{\kern
0.0ptk}_{T},z,\mu,\nu)\,f_{j}\Bigl{(}\frac{\omega}{zE_{\mathrm{cm}}},\mu\Bigr{)}\bigg{[}1+\mathcal{O}\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}^{2}}{k_{T}^{2}}\Bigr{)}\bigg{]}\,.$
(2.23)

The double-differential soft function $S_{\kappa}(k,\vec{\kern
0.0ptk}_{s},\mu,\nu)$ encodes the contribution of soft radiation to both
${\mathcal{T}}$ and $q_{T}$. The RG consistency of the cross section implies
that its $\mu$ and $\nu$ RGEs do not depend on ${\mathcal{T}}$. Hence, the
overall RG structure of the double-differential cross section is equivalent to
the single-differential $q_{T}$ spectrum, with ${\mathcal{T}}$ acting as an
auxiliary measurement.

The factorization in eq. (2.20) separates the physics at the canonical SCETII
invariant-mass and rapidity scales

$\displaystyle\mu_{H}^{\mathrm{II}}$ $\displaystyle\sim Q\,,\qquad$
$\displaystyle\mu_{B}^{\mathrm{II}}$ $\displaystyle\sim q_{T}\,,\qquad$
$\displaystyle\mu_{S}^{\mathrm{II}}$ $\displaystyle\sim q_{T}\,,$
$\displaystyle\nu_{B}^{\mathrm{II}}$ $\displaystyle\sim Q\,,\qquad$
$\displaystyle\nu_{S}^{\mathrm{II}}$ $\displaystyle\sim q_{T}\,.$ (2.24)

It has been known for a long time [75] that directly resumming the logarithms
of $q_{T}/Q$ in momentum space is challenging due to the vectorial nature of
$\vec{\kern 0.0ptq}_{T}$, though by now approaches for doing so exist [76,
77]. The same complications arise here for the double-differential spectrum.
We bypass this issue, as is commonly done, by carrying out the resummation in
conjugate ($b_{T}$) space [78, 79, 80, 72]. The Fourier transform from
$\vec{\kern 0.0ptq}_{T}$ to $\vec{\kern 0.0ptb}_{T}$ turns the vectorial
convolutions in eq. (2.20) into simple products at
$b_{T}\equiv\lvert\vec{\kern 0.0ptb}_{T}\rvert$. The canonical SCETII scales
in $b_{T}$-space are then given by

$\displaystyle\mu_{H}^{\mathrm{II}}$ $\displaystyle\sim Q\,,\qquad$
$\displaystyle\mu_{B}^{\mathrm{II}}$ $\displaystyle\sim b_{0}/b_{T}\,,\qquad$
$\displaystyle\mu_{S}^{\mathrm{II}}$ $\displaystyle\sim b_{0}/b_{T}\,,$
$\displaystyle\nu^{\mathrm{II}}_{B}$ $\displaystyle\sim Q\,,\qquad$
$\displaystyle\nu^{\mathrm{II}}_{S}$ $\displaystyle\sim b_{0}/b_{T}\,,$ (2.25)

where $b_{0}\equiv 2e^{-\gamma_{E}}\approx 1.12291$ is conventional. By
evaluating the functions in the factorization theorem at their canonical
scales and evolving them to a common scale in both $\mu$ and $\nu$, all
logarithms of
$\mu_{B}/\mu_{H}\sim\mu_{S}/\mu_{H}\sim\nu_{S}/\nu_{B}\sim(b_{0}/b_{T})/Q$ are
resummed. In ref. [77] it was shown that the canonical resummation in $b_{T}$
space is in fact equivalent to the exact solution of the RGE in momentum
space, except for the fact that one effectively uses a shifted set of finite
terms in the boundary conditions (similar to the difference between
renormalization schemes). We exploit this and require that for $q_{T}\ll Q$,
eq. (2.3) is exactly satisfied, such that the resummed $q_{T}$ spectrum in
this region is obtained from the inverse Fourier transform of the canonical
$b_{T}$-space result.

A key feature of the resummed $q_{T}$ spectrum is that the anomalous dimension
$\gamma^{i}_{\nu}$, driving the $\nu$ running of the soft (or beam) function
at fixed $\mu$, is itself perturbatively renormalized at its intrinsic scale
$\mu_{0}$ and requires resummation when $\mu_{0}\neq\mu$. Specifically, in the
exponent of the $b_{T}$-space rapidity evolution factor we have

$\displaystyle\tilde{\gamma}^{i}_{\nu}(b_{T},\mu)$
$\displaystyle=-4\eta^{i}_{\Gamma}(\mu_{0},\mu)+\tilde{\gamma}^{i}_{\nu,\mathrm{FO}}(b_{T},\mu_{0})+\tilde{\gamma}^{i}_{\nu,\mathrm{np}}(b_{T})\,,$
(2.26)

where all logarithms of $\mu/\mu_{0}$ are resummed inside

$\eta_{\Gamma}^{i}(\mu_{0},\mu)=\int_{\mu_{0}}^{\mu}\\!\frac{\mathrm{d}\mu^{\prime}}{\mu^{\prime}}\,\Gamma^{i}_{\rm
cusp}[\alpha_{s}(\mu^{\prime})]\,.$ (2.27)

[See eq. (4.26) in ref. [77] for the analogous expressions in momentum space.]
The canonical choice of $\mu_{0}$ that eliminates all large logarithms in the
fixed-order boundary condition
$\tilde{\gamma}_{\nu}^{\mathrm{FO}}(b_{T},\mu_{0})$ is

$\mu_{0}\sim b_{0}/b_{T}\,.$ (2.28)

By choosing $\mu_{0}$ as a function of $b_{T}$ such that it freezes out to a
perturbative value at large $b_{T}$, we avoid the Landau pole at $b_{T}\sim
1/\Lambda_{\mathrm{QCD}}$.333In addition, this leaves fixed-order logarithms
of $\mu_{0}b_{T}$ in $\tilde{\gamma}^{i}_{\nu,\mathrm{FO}}$ that lead to an
exponential suppression of the $b_{T}$ space cross section as
$b_{T}\to\infty$. This increases the numerical stability of the inverse
Fourier transform. The mismatch to the full result can in principle be
captured by a nonperturbative model
$\tilde{\gamma}^{i}_{\nu,\mathrm{np}}(b_{T})$, which can be extracted from
experimental measurements at small $q_{T}$. Recently, it was shown that it
could also be determined from lattice calculations [81]. For our purposes we
set $\tilde{\gamma}^{i}_{\nu,\mathrm{np}}=0$ for simplicity. (We similarly
ignore nonperturbative effects in the SCETII beam and soft function boundary
conditions.) Our concrete choice of $\mu_{0}$ is given below.

#### Scale setting and fixed-order matching.

Figure -1016: Left: Comparison of singular and nonsingular contributions to
the fixed $\mathcal{O}(\alpha_{s})$ double spectrum as a function of $q_{T}$,
with ${\mathcal{T}}=q_{T}$ kept fixed. The orange solid line shows the full
QCD result and the dashed blue line the singular contributions contained in
the SCETII result eq. (2.20). The dotted green line shows their difference,
which corresponds to the power corrections indicated in eq. (2.21). The thin
vertical lines correspond to the transition points $(x_{1},x_{2},x_{3})$ given
in the text. Right: SCETII hybrid profile scales as a function of
$b_{0}/b_{T}$ for representative values of $q_{T}$. The thin vertical line in
the main plot corresponds to $Q$. The inset shows the behavior of the profile
in the nonperturbative region $b_{0}/b_{T}\sim\Lambda_{\mathrm{QCD}}$, where
the gray horizontal line indicates the scale $b_{0}/b_{\mathrm{max}}$ at which
we freeze out the resummation. The dashed orange line in the inset indicates
the canonical value of
$\mu_{B}^{\mathrm{II}},\mu_{S}^{\mathrm{II}},\nu_{S}^{\mathrm{II}}$.

We again extend the description of the cross section to the fixed-order region
$q_{T}\sim{\mathcal{T}}\lesssim Q$ by an additive matching,

$\mathrm{d}\sigma^{\mathrm{match}}_{\mathrm{II}}=\mathrm{d}\sigma_{\mathrm{II}}\bigr{|}_{\mu^{\mathrm{II}}}+\bigl{[}\mathrm{d}\sigma_{\mathrm{FO}}-\mathrm{d}\sigma_{\mathrm{II}}\bigr{]}_{\mu_{\mathrm{FO}}}\,.$
(2.29)

Here the subscript $\mu^{\mathrm{II}}$ indicates that we evaluate
$\mathrm{d}\sigma_{\mathrm{II}}$ at the SCETII resummation scales
$\mu^{\mathrm{II}}$ (given below) in $b_{T}$ space, and take a numerical
inverse Fourier transform in the end. The subscript $\mu_{\mathrm{FO}}$
indicates that it is instead evaluated at common fixed-order scales
$\mu_{\mathrm{FO}}$, which can be done directly in momentum space.

Analogous to the discussion for SCETI, the term in square brackets in eq.
(2.29) is by construction a pure nonsingular power correction at small
$q_{T}$. This is illustrated in the left panel of figure -1016, which shows
that the difference (green dotted) between the full QCD result (solid orange)
and the SCETII singular result (dashed blue) indeed vanishes like a power as
$q_{T}\to 0$ along the line of fixed ${\mathcal{T}}=q_{T}$.

Approaching $q_{T}\sim{\mathcal{T}}\sim Q$, the $q_{T}$ resummation must again
be turned off to ensure the delicate cancellations between singular and
nonsingular contributions and to properly recover the correct fixed-order
result for the spectrum. We achieve this by constructing hybrid profile scales
that depend on both $b_{T}$ and $q_{T}$, and undergo a continuous deformation
away from the canonical $b_{T}$ scales in eq. (2.3) as a function of the
target $q_{T}$ value, schematically,

$\mu_{B,S}^{\mathrm{II}}(q_{T},b_{T})\,,\nu_{B,S}^{\mathrm{II}}(q_{T},b_{T})\to\mu^{\mathrm{II}}_{H}=\mu_{\mathrm{FO}}\quad\text{for}\quad
q_{T}\to Q\,.$ (2.30)

We note that $\mu_{0}$ does not need to asymptote to $\mu_{\mathrm{FO}}$
towards large $q_{T}$ because its effect on the matched result is already
turned off as $\nu_{S}^{\mathrm{II}}\to\nu_{B}^{\mathrm{II}}$. In this limit,
the first and last term in eq. (2.29) exactly cancel, leaving the fixed-order
result $\mathrm{d}\sigma_{\mathrm{FO}}$.

Since the single-differential $q_{T}$ resummation is not the main focus of
this paper, we strive to achieve eq. (2.30) in the simplest possible way.
Specifically, we choose central scales as

$\displaystyle\mu^{\mathrm{II}}_{H}=\nu^{\mathrm{II}}_{B}={\mu_{\mathrm{FO}}}\,,\qquad\mu^{\mathrm{II}}_{B}=\mu^{\mathrm{II}}_{S}=\nu^{\mathrm{II}}_{S}={\mu_{\mathrm{FO}}}\,f_{\mathrm{run}}^{\rm
II}\Bigl{(}\frac{q_{T}}{Q},\frac{b_{0}}{b^{\ast}(b_{T})\,Q}\Bigr{)}\,,\qquad\mu_{0}=\frac{b_{0}}{b^{\ast}(b_{T})}\,,$
(2.31)

where $f_{\mathrm{run}}^{\mathrm{II}}$ is a hybrid profile function given by

$\displaystyle f_{\mathrm{run}}^{\mathrm{II}}(x,y)$
$\displaystyle=1+g_{\mathrm{run}}(x)(y-1)\,.$ (2.32)

It controls the amount of resummation by adjusting the slope of the scales in
$b_{T}$ space as a function of $q_{T}/Q$ via the function

$\displaystyle g_{\mathrm{run}}(x)$ $\displaystyle=\begin{cases}1&0<x\leq
x_{1}\,,\\\ 1-\frac{(x-x_{1})^{2}}{(x_{2}-x_{1})(x_{3}-x_{1})}&x_{1}<x\leq
x_{2}\,,\\\ \frac{(x-x_{3})^{2}}{(x_{3}-x_{1})(x_{3}-x_{2})}&x_{2}<x\leq
x_{3}\,,\\\ 0&x_{3}\leq x\,.\end{cases}$ (2.33)

As a result, for $q_{T}\leq x_{1}Q$, the slope is unity yielding the canonical
resummation, while for $q_{T}\geq x_{3}Q$, the slope vanishes so the
resummation is fully turned off. In between, the slope smoothly transitions
from one to zero, which transitions the resummation from being canonical to
being turned off. This is illustrated in the right panel of figure -1016. We
use the same transition points $(x_{1},x_{2},x_{3})=(0.2,0.5,0.8)$ as for
SCETI, which is supported by figure -1016.

We note that our approach differs from the hybrid profile scales introduced in
ref. [82]. While the latter also satisfy the requirement in eq. (2.30), they
do not reproduce the exact canonical $b_{T}$-space scales for $q_{T}\ll Q$
because they introduce a profile shape directly in $b_{T}$ space.

As discussed below eq. (2.28), we require a nonperturbative prescription when
the canonical value of $\mu_{0}$ (or $\mu_{S}^{\mathrm{II}}$, or
$\mu_{B}^{\mathrm{II}}$) approaches the Landau pole
$b_{0}/b_{T}\sim\Lambda_{\mathrm{QCD}}$. This is encoded in evaluating the
hybrid scales at $b^{\ast}(b_{T})$ rather than $b_{T}$ itself,

$\displaystyle
b^{\ast}(b_{T})=\frac{b_{T}}{\sqrt{1+b_{T}^{2}/b_{\mathrm{max}}^{2}}}\,,$
(2.34)

where $b_{0}/b_{\mathrm{max}}\gtrsim\Lambda_{\mathrm{QCD}}$ ensures that all
scales are canonical for small $b_{T}\approx b^{\ast}$, but remain
perturbative for large $b_{T}$ where $b^{\ast}\to b_{\mathrm{max}}$, as shown
in the inset in the right panel of figure -1016. In practice we pick

$b_{0}/b_{\mathrm{max}}=1\,\mathrm{GeV}\,,$ (2.35)

in keeping with our choice of nonperturbative turn-off parameter in the SCETI
case. The functional form of eq. (2.34) is the same as in the standard
$b^{\ast}$ prescription [79, 80], although any other functional form with the
same asymptotic behavior is also viable. We stress, however, that a key
difference in our case is that $b^{\ast}$ only affects the scales, so it
essentially serves the same purpose as the $x_{0}$ nonperturbative cutoff in
the SCETI scales in eq. (2.14). By contrast, the standard $b^{\ast}$
prescription corresponds to a _global_ replacement of $b_{T}$ by $b^{\ast}$,
including the measurement itself. For the single-differential $q_{T}$
spectrum, this global replacement induces power corrections
$\mathcal{O}(b_{T}^{2}/b_{\mathrm{max}}^{2})$ that scale like a generic
nonperturbative contribution. While they might complicate the extraction of
nonperturbative model parameters from data [83], they are not a critical
issue.

For the double-differential case, we find that a standard $b^{\ast}$
prescription does in fact not work. This is because substituting $b^{\ast}$
for $b_{T}$ in the physical measurement renders Fourier integrals of the
double-differential SCETII soft function divergent, at least at fixed order
(i.e., without Sudakov suppression). This can be seen from eqs. (B.40) and
(B.41), which only depend on $x=b_{T}{\mathcal{T}}$. Substituting $b^{\ast}$
for $b_{T}$ makes them asymptote to a constant for any given ${\mathcal{T}}$,
which upsets their required asymptotic behavior $\sim 1/x^{2}$. Physically
this means that the deformation of the measurement at large $b_{T}$ also
deforms the observable of interest, i.e., the dependence on ${\mathcal{T}}$.

#### Perturbative uncertainties.

To estimate the resummation uncertainty for
$\mathrm{d}\sigma_{\mathrm{II}}^{\mathrm{match}}$, we adopt the set of profile
scale variations introduced for the SCETII-like jet veto in ref. [67]. They
are given by

$\displaystyle\mu_{S}^{\mathrm{II}}$
$\displaystyle={\mu_{\mathrm{FO}}}\,\Bigl{[}f_{\mathrm{vary}}\Bigl{(}\frac{q_{T}}{Q}\Bigr{)}\Bigr{]}^{v_{\mu_{S}}}\,f_{\mathrm{run}}^{\mathrm{II}}\Bigl{(}\frac{q_{T}}{Q},\frac{b_{0}}{b^{\ast}Q}\Bigr{)}\,,$
$\displaystyle\nu_{S}^{\mathrm{II}}$
$\displaystyle={\mu_{\mathrm{FO}}}\,\Bigl{[}f_{\mathrm{vary}}\Bigl{(}\frac{q_{T}}{Q}\Bigr{)}\Bigr{]}^{v_{\nu_{S}}}\,f_{\mathrm{run}}^{\mathrm{II}}\Bigl{(}\frac{q_{T}}{Q},\frac{b_{0}}{b^{\ast}Q}\Bigr{)}\,,$
$\displaystyle\mu_{B}^{\mathrm{II}}$
$\displaystyle={\mu_{\mathrm{FO}}}\,\Bigl{[}f_{\mathrm{vary}}\Bigl{(}\frac{q_{T}}{Q}\Bigr{)}\Bigr{]}^{v_{\mu_{B}}}\,f_{\mathrm{run}}^{\mathrm{II}}\Bigl{(}\frac{q_{T}}{Q},\frac{b_{0}}{b^{\ast}Q}\Bigr{)}\,,$
$\displaystyle\nu_{B}^{\mathrm{II}}$
$\displaystyle={\mu_{\mathrm{FO}}}\,\Bigl{[}f_{\mathrm{vary}}\Bigl{(}\frac{q_{T}}{Q}\Bigr{)}\Bigr{]}^{v_{\nu_{B}}}\,,$
(2.36)

where each of the four variation exponents can be $v_{i}=\\{+1,0,-1\\}$, and
$f_{\mathrm{vary}}$ was given in eq. (2.16). The central scale choice
corresponds to $(v_{\mu_{S}},v_{\nu_{S}},v_{\mu_{B}},v_{\nu_{B}})=(0,0,0,0)$,
and a priori there are 80 possible different combinations of the $v_{i}$.
Since the arguments of the resummed logarithms are ratios of scales, some
combinations of scale variations will lead to variations of these arguments
that are larger than a factor of two, and therefore should be excluded [67].
After dropping these combinations we are left with 36 different scale
variations for the SCETII regime. We add two independent variations of
$b_{0}/b_{\mathrm{max}}=\\{0.5\,\mathrm{GeV},2\,\mathrm{GeV}\\}$ to probe the
uncertainty in our nonperturbative prescription. The SCETII resummation
uncertainty $\Delta_{\mathrm{II}}$ is then determined as the maximum absolute
deviation from the central result among all 38 variations. For simplicity we
again refrain from variations of the transition points. As for SCETI,
$\Delta_{\mathrm{FO}}$ is estimated by overall variations of
$\mu_{\mathrm{FO}}$ by a factor of two, which is inherited by all SCETII
scales, so it probes the fixed-order uncertainties while leaving the resummed
logarithms invariant. The total uncertainty estimate for
$\mathrm{d}\sigma_{\mathrm{II}}^{\mathrm{match}}$ is then obtained as

$\Delta^{\mathrm{II}}_{\mathrm{total}}=\Delta_{\mathrm{II}}\oplus\Delta_{\mathrm{FO}}\,.$
(2.37)

The matched result $\mathrm{d}\sigma_{\mathrm{II}}^{\mathrm{match}}$ in eq.
(2.29) provides a prediction for the double-differential spectrum that covers
the part of phase space where ${\mathcal{T}}\sim q_{T}$.

#### Results for the single-differential spectrum.

Figure -1015: The single-differential $q_{T}$ spectrum at NLL (blue) and
NNLL$+$NLO (orange), using the $q_{T}$ resummation method described in the
text. The bands indicate $\Delta_{\mathrm{II}}\oplus\Delta_{\mathrm{FO}}$. In
the right panel, the uncertainties are shown as percent differences relative
to the central result at each order.

Since we are using a new method to perform the $q_{T}$ resummation, we also
briefly consider the single-differential $q_{T}$ spectrum as a sanity check of
our setup. The setup described in this section immediately carries over to the
single-differential spectrum. In figure -1015 we show the $q_{T}$ spectrum at
the NNLL$+$NLO order we are aiming for for the double-differential spectrum,
as well as one order lower at NLL, and with the uncertainties estimated as
described above. The results look very reasonable, providing us with
confidence in our $q_{T}$ resummation procedure. Note that there is a slight
pinch in the uncertainty bands around $q_{T}=15\,\mathrm{GeV}$, indicating
that the uncertainties there are likely a bit underestimated. This is an
artifact of scale variations that is not unusual to be seen in resummed
spectrum predictions.

###  2.4 SCET+: ${\mathcal{T}}\ll q_{T}\ll\sqrt{Q{\mathcal{T}}}$

This regime is characterized by the presence of intermediate collinear-soft
modes that contribute both to the $q_{T}$ and the ${\mathcal{T}}$ measurement,
which uniquely fixes their scaling. Central soft modes only contribute to
${\mathcal{T}}$ as in SCETI, while the energetic collinear modes only
contribute to $q_{T}$ as in SCETII,

$\displaystyle n_{a}\text{-collinear:}$ $\displaystyle\quad
p^{\mu}\sim\Bigl{(}\tfrac{q_{T}^{2}}{Q},Q,q_{T}\Bigr{)}\,,\qquad$
$\displaystyle n_{a}\text{-collinear-soft:}\quad
p^{\mu}\sim\Bigl{(}{\mathcal{T}},\frac{q_{T}^{2}}{{\mathcal{T}}},q_{T}\Bigr{)}\,,$
$\displaystyle n_{b}\text{-collinear:}$ $\displaystyle\quad
p^{\mu}\sim\Bigl{(}Q,\tfrac{q_{T}^{2}}{Q},q_{T}\Bigr{)}\,,\qquad$
$\displaystyle n_{b}\text{-collinear-soft:}\quad
p^{\mu}\sim\Bigl{(}\frac{q_{T}^{2}}{{\mathcal{T}}},{\mathcal{T}},q_{T}\Bigr{)}\,,$
soft: $\displaystyle\quad
p^{\mu}\sim\Bigl{(}{\mathcal{T}},{\mathcal{T}},{\mathcal{T}}\Bigr{)}\,.$
(2.38)

The collinear-soft modes have the same virtuality as the collinear modes,
$p^{2}\sim q_{T}^{2}$, but live at more central rapidity $e^{\lvert
y\rvert}\sim q_{T}/{\mathcal{T}}$, which is small compared to the rapidity
$e^{\lvert y\rvert}\sim Q/q_{T}$ of the collinear modes. Hence, the two have a
SCETII-like relation and become a single collinear mode in the SCETI limit
$q_{T}\sim\sqrt{Q{\mathcal{T}}}$. At the same time, the collinear-soft and
soft modes have a SCETI-like relation, being separated in virtuality, and
become a single soft mode in the SCETII limit ${\mathcal{T}}\sim q_{T}$. In
this way, SCET+ is able to connect the SCETI and SCETII regimes. This is
similar to the collinear-soft mode originally introduced in ref. [28], which
instead connected two SCETI theories.

The cross section in SCET+ factorizes as [51]

$\displaystyle\frac{\mathrm{d}\sigma_{+}}{\mathrm{d}Q\,\mathrm{d}Y\,\mathrm{d}q_{T}\,\mathrm{d}{\mathcal{T}}}$
$\displaystyle=H_{\kappa}(Q,\mu)\int\\!\mathrm{d}^{2}\vec{\kern
0.0ptk}_{a}\,B_{a}(\omega_{a},\vec{\kern
0.0ptk}_{a},\mu,\nu)\int\\!\mathrm{d}^{2}\vec{\kern
0.0ptk}_{b}\,B_{b}(\omega_{b},\vec{\kern 0.0ptk}_{b},\mu,\nu)$ (2.39)
$\displaystyle\quad\times\int\\!\mathrm{d}\ell_{a}^{+}\int\\!\mathrm{d}^{2}\vec{\kern
0.0pt\ell}_{a}\,\mathcal{S}_{\kappa}(\ell_{a}^{+},\vec{\kern
0.0pt\ell}_{a},\mu,\nu)\int\\!\mathrm{d}\ell_{b}^{-}\int\\!\mathrm{d}^{2}\vec{\kern
0.0pt\ell}_{b}\,\mathcal{S}_{\kappa}(\ell_{b}^{-},\vec{\kern
0.0pt\ell}_{b},\mu,\nu)$
$\displaystyle\quad\times\int\\!\mathrm{d}k\,S_{\kappa}(k,\mu)\,\delta\bigl{(}q_{T}-\lvert\vec{\kern
0.0ptk}_{a}+\vec{\kern 0.0ptk}_{b}+\vec{\kern 0.0pt\ell}_{a}+\vec{\kern
0.0pt\ell}_{b}\rvert\bigr{)}\,\delta\Bigl{(}{\mathcal{T}}-\frac{\omega_{a}\ell_{a}^{+}}{Q_{a}}-\frac{\omega_{b}\ell_{b}^{-}}{Q_{b}}-k\Bigr{)}\,,$

which holds up to power corrections

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}Q\,\mathrm{d}Y\,\mathrm{d}q_{T}\,\mathrm{d}{\mathcal{T}}}$
$\displaystyle=\frac{\mathrm{d}\sigma_{+}}{\mathrm{d}Q\,\mathrm{d}Y\,\mathrm{d}q_{T}\,\mathrm{d}{\mathcal{T}}}\,\Big{[}1+\mathcal{O}\Bigl{(}\frac{q_{T}^{2}}{{\mathcal{T}}Q},\frac{{\mathcal{T}}^{2}}{q_{T}^{2}}\Bigr{)}\Big{]}\,.$
(2.40)

The hard function is the same as before. The beam functions are the
$q_{T}$-dependent ones from SCETII, while the soft function is the
${\mathcal{T}}$-dependent one from SCETI. The new ingredient is the double-
differential collinear-soft function $\mathcal{S}_{\kappa}(k,\vec{\kern
0.0ptk}_{T},\mu,\nu)$, which encodes the contributions of the collinear-soft
modes to both $q_{T}$ and ${\mathcal{T}}$. Like the soft function it is
defined as a matrix element of eikonal Wilson lines, but like the beam
functions it describes radiation that goes into a definite hemisphere.

Equation (2.39) can be interpreted as a refactorization of the double-
differential SCETI and SCETII cross sections [51], which precisely reflects
the relation between the EFT modes described above. Expanding the SCETI
double-differential beam function in the limit
$q_{T}\ll\sqrt{Q{\mathcal{T}}}$, it factorizes into the SCETII beam function
and the collinear-soft function,

$\displaystyle B_{q}(\omega k,\omega/E_{\mathrm{cm}},\vec{\kern
0.0ptk}_{T},\mu)$ $\displaystyle=\int\\!\mathrm{d}^{2}\vec{\kern
0.0pt\ell}_{T}\,B_{q}(\omega,\vec{\kern 0.0ptk}_{T}-\vec{\kern
0.0pt\ell}_{T},\mu,\nu)\,\mathcal{S}_{\kappa}(k,\vec{\kern
0.0pt\ell}_{T},\mu,\nu)\,\Bigl{[}1+\mathcal{O}\Bigl{(}\frac{k_{T}^{2}}{\omega
k}\Bigr{)}\Bigr{]}\,.$ (2.41)

The $\nu$ dependence of the two terms on the right-hand side must cancel,
while their $\mu$ dependence must combine into that of the left-hand side.
This allows us to derive the RGE for the collinear-soft function given in eq.
(B.1).

Similarly, expanding the SCETII double-differential soft function in the limit
${\mathcal{T}}\ll q_{T}$, it factorizes into the SCETI soft function and the
two $n_{a}$-collinear-soft and $n_{b}$-collinear-soft functions,

$\displaystyle S_{\kappa}(k,\vec{\kern 0.0ptk}_{T},\mu,\nu)$
$\displaystyle=\int\\!\mathrm{d}^{2}\vec{\kern
0.0pt\ell}_{T}\,\int\\!\mathrm{d}\ell_{a}^{+}\,\mathcal{S}_{\kappa}(\ell_{a}^{+},\vec{\kern
0.0pt\ell}_{T},\mu,\nu)\,\int\\!\mathrm{d}\ell_{b}^{-}\,\mathcal{S}_{\kappa}(\ell_{b}^{-},\vec{\kern
0.0ptk}_{T}-\vec{\kern 0.0pt\ell}_{T},\mu,\nu)\,$ $\displaystyle\quad\times
S_{\kappa}\Bigl{(}k-\frac{\omega_{a}\ell_{a}^{+}}{Q_{a}}-\frac{\omega_{b}\ell_{b}^{-}}{Q_{b}}\Bigr{)}\,\Bigl{[}1+\mathcal{O}\Bigl{(}\frac{k^{2}}{k_{T}^{2}}\Bigr{)}\Bigr{]}\,.$
(2.42)

Since the left-hand side does not depend on $\omega_{a,b}$ and $Q_{a,b}$, this
dependence must also drop out on the right-hand side, and therefore in the
whole SCET+ cross section in eq. (2.39). To see this explicitly, first recall
that $\omega_{a}\omega_{b}=Q_{a}Q_{b}=Q^{2}$. In addition, boost invariance at
the level of the collinear-soft matrix element implies that
$\mathrm{d}\ell_{a}^{+}\,\mathcal{S}_{\kappa}(\ell_{a}^{+},\vec{\kern
0.0ptk}_{T},\mu,\nu)$ can only depend on the product $\ell_{a}^{+}\nu$ (and
analogously for $\ell_{b}^{-}$).444More explicitly, the rapidity regulator
breaks the RPI-III invariance of SCET [84, 85], which is equivalent to boost
invariance that must hold separately in each collinear sector. To restore it,
$\nu$ must transform under RPI-III like ${\bar{n}}\cdot p$ in each
$n$-collinear-soft sector. This is most straightforward to see when strictly
expanding the rapidity regulator to leading power in
${\mathcal{T}}^{2}/q_{T}^{2}$ using the soft-collinear mode scaling in eq.
(2.4). The RPI-III transformation of the explicit measurement $\delta$
function in the matrix element is canceled by the corresponding integration
measure in eqs. (2.39) and (2.4). Therefore, RPI-III invariance implies that
each collinear-soft function can only depend on the RPI-III invariant
combination $\nu\,n\cdot k$. Hence, we can rewrite

$\displaystyle\mathrm{d}\ell_{a}^{+}\,\mathcal{S}_{\kappa}(\ell_{a}^{+},\vec{\kern
0.0pt\ell}_{a},\mu,\nu)\,\mathrm{d}\ell_{b}^{-}\,\mathcal{S}_{\kappa}(\ell_{b}^{-},\vec{\kern
0.0pt\ell}_{b},\mu,\nu)\,\delta\Bigl{(}{\mathcal{T}}-\frac{\omega_{a}\ell_{a}^{+}}{Q_{a}}-\frac{\omega_{b}\ell_{b}^{-}}{Q_{b}}-k\Bigr{)}$
$\displaystyle\qquad=\mathrm{d}k_{a}^{+}\,\mathcal{S}_{\kappa}\Bigl{(}k_{a}^{+},\vec{\kern
0.0pt\ell}_{a},\mu,\frac{Q_{a}\nu}{\omega_{a}}\Bigr{)}\,\mathrm{d}k_{b}^{-}\,\mathcal{S}_{\kappa}\Bigl{(}k_{b}^{-},\vec{\kern
0.0pt\ell}_{b},\mu,\frac{Q_{b}\nu}{\omega_{b}}\Bigr{)}\,\delta({\mathcal{T}}-k_{a}^{+}-k_{b}^{-}-k)$
$\displaystyle\qquad=\mathrm{d}k_{a}^{+}\,\mathcal{S}_{\kappa}(k_{a}^{+},\vec{\kern
0.0pt\ell}_{a},\mu,\nu)\,\mathrm{d}k_{b}^{-}\,\mathcal{S}_{\kappa}(k_{b}^{-},\vec{\kern
0.0pt\ell}_{b},\mu,\nu)\,\delta({\mathcal{T}}-k_{a}^{+}-k_{b}^{-}-k)\,,$
(2.43)

where in the first step we changed variables from $\ell_{a,b}^{\pm}$ to
$k_{a}^{+}=\omega_{a}\ell_{a}^{+}/Q_{a}$ and
$k_{b}^{-}=\omega_{b}\ell_{b}^{-}/Q_{b}$. In the second step we performed the
rapidity evolution from $\nu_{a,b}\equiv Q_{a,b}\,\nu/\omega_{a,b}$ back to a
common $\nu$ at fixed $\mu$ [see eq. (B.3)], for which the rapidity evolution
factors exactly cancel because

$\ln\frac{\nu_{a}}{\nu}+\ln\frac{\nu_{b}}{\nu}=\ln\frac{Q_{a}Q_{b}}{\omega_{a}\omega_{b}}=0\,.$
(2.44)

The SCET+ factorization in eq. (2.39) fully disentangles the physics at the
following SCET+ canonical energy and rapidity scales:

$\displaystyle\mu^{+}_{H}$ $\displaystyle\sim Q\,,\qquad$
$\displaystyle\mu^{+}_{B}$ $\displaystyle\sim q_{T}\,,\qquad$
$\displaystyle\mu^{+}_{\mathcal{S}}$ $\displaystyle\sim q_{T}\,,\qquad$
$\displaystyle\mu^{+}_{S}$ $\displaystyle\sim{\mathcal{T}}\,,$
$\displaystyle\nu^{+}_{B}$ $\displaystyle\sim Q\,,\qquad$
$\displaystyle\nu^{+}_{\mathcal{S}}$ $\displaystyle\sim
q_{T}^{2}/{\mathcal{T}}\,.\quad$ (2.45)

As for SCETII, we perform the $q_{T}$ resummation in $b_{T}$ space,
transforming the vectorial convolutions in eq. (2.39) into simple products. In
$b_{T}$ space, the canonical SCET+ scales are

$\displaystyle\mu^{+}_{H}$ $\displaystyle\sim Q\,,\qquad$
$\displaystyle\mu^{+}_{B}$ $\displaystyle\sim b_{0}/b_{T}\,,\qquad$
$\displaystyle\mu^{+}_{\mathcal{S}}$ $\displaystyle\sim b_{0}/b_{T}\,,\qquad$
$\displaystyle\mu^{+}_{S}$ $\displaystyle\sim{\mathcal{T}}\,,$
$\displaystyle\nu^{+}_{B}$ $\displaystyle\sim Q\,,\qquad$
$\displaystyle\nu^{+}_{\mathcal{S}}$
$\displaystyle\sim(b_{0}/b_{T})^{2}/{\mathcal{T}}\,.\quad$ (2.46)

By evaluating all functions at their natural scales and evolving them to
common scales, all logarithms of large scale ratios in the problem are
resummed, e.g.,

$\frac{(b_{0}/b_{T})^{2}}{Q{\mathcal{T}}}\sim\frac{\nu^{+}_{\mathcal{S}}}{\nu^{+}_{B}}\,,\qquad\frac{{\mathcal{T}}}{b_{0}/b_{T}}\sim\frac{\mu^{+}_{S}}{\mu^{+}_{\mathcal{S}}}\,,\qquad\frac{b_{0}/b_{T}}{Q}\sim\frac{\mu^{+}_{B}}{\mu^{+}_{H}}\sim\frac{\mu^{+}_{\mathcal{S}}}{\mu^{+}_{H}}\,,\qquad\frac{{\mathcal{T}}}{Q}\sim\frac{\mu^{+}_{S}}{\mu^{+}_{H}}\,.$
(2.47)

The logarithms of the first ratio appear in the double-differential SCETI beam
function in the limit $q_{T}\ll\sqrt{Q{\mathcal{T}}}$, and are resummed in
SCET+ by the additional $\nu$ evolution in the refactorization in eq. (2.41).
Similarly, logarithms of the second ratio appear in the double-differential
SCETII soft function in the limit ${\mathcal{T}}\ll q_{T}$, and are resummed
in SCET+ by the additional $\mu$ evolution in eq. (2.4). Our framework to
match between the rich logarithmic structure predicted by eq. (2.39) and the
two boundary regimes is the subject of section 3.

###  2.5 Outer space

We now briefly discuss the outer phase-space regions left blank in figure
-1018. The region above the SCETII regime is characterized by the hierarchy
$q_{T}\ll{\mathcal{T}}\ll\sqrt{Q{\mathcal{T}}}$, while the region to the right
of the SCETI regime corresponds to ${\mathcal{T}}\ll\sqrt{Q{\mathcal{T}}}\ll
q_{T}$. Both regions are power suppressed.

As we have discussed in section 2.3, only the soft function contributes to
${\mathcal{T}}$ in SCETII, as the collinear contribution is power suppressed.
However, for $q_{T}\ll{\mathcal{T}}$, even the soft contribution to
${\mathcal{T}}$ becomes power suppressed. In particular, for a single real
emission at fixed $\mathcal{O}(\alpha_{s})$, the region ${\mathcal{T}}>q_{T}$
is kinematically forbidden both in SCETII as well as in full QCD. At higher
orders only (soft) emissions that are mostly back-to-back such that their
transverse momenta largely cancel can fill out this region. The cross section
in this region is power suppressed by
$\mathcal{O}(q_{T}^{2}/{\mathcal{T}}^{2})$. Equivalently, expanding the SCETII
factorization of the double-differential cross section in the limit
$q_{T}\ll{\mathcal{T}}$ reduces it to the single-differential $q_{T}$ spectrum
with an overall $\delta({\mathcal{T}})$, which we exploit in our numerical
implementation, cf. eq. (B.33). Physically this means that by integrating the
double spectrum in SCETII up to some ${\mathcal{T}}_{\mathrm{cut}}\gg q_{T}$,
we recover the single-differential $q_{T}$ spectrum, while the effect of the
cut is power suppressed in this limit. Note that there is also a contribution
from double-parton scattering [86, 87, 88, 89] in this region, where the two
jets produced in the second interaction alongside the $Z$ boson are naturally
back to back and not power suppressed. This contribution is still not expected
to much exceed the single-parton scattering contribution because double-parton
scattering itself is power suppressed by $\mathcal{O}(\Lambda_{\rm
QCD}^{2}/{\mathcal{T}}^{2})$, with ${\mathcal{T}}$ setting the scale of the
second hard scatter producing the back-to-back jets.

Similarly, in the limit $\sqrt{Q{\mathcal{T}}}\ll q_{T}$, even the
contribution from collinear radiation to $q_{T}$ becomes power suppressed in
SCETI [cf. eq. (B.24)], and at leading power we recover the single-
differential ${\mathcal{T}}$ spectrum with an overall $\delta(q_{T})$. This is
analogous to the relation between the regimes 1 and 2 for a jet veto with a
jet rapidity cut in ref. [34], where the effect of a very forward jet rapidity
cut (the auxiliary measurement) on collinear radiation becomes power
suppressed. An additional subtlety for $\sqrt{Q{\mathcal{T}}}\ll q_{T}$ is
that very energetic forward radiation with energy $\sim
q_{T}^{2}/{\mathcal{T}}$ can theoretically contribute [51], pushing the hard
scale up to $q_{T}^{2}/{\mathcal{T}}\gg Q$. However, the cross section in this
kinematic configuration is very strongly suppressed by the PDFs, so we choose
to describe it at fixed order in this paper.

The above analysis justifies focusing on the shaded regions of phase space in
figure -1018, corresponding to the main SCETI, SCETII, and SCET+ regimes.

##  3 Matching effective theories

###  3.1 Structure of power corrections

An important feature of our EFT setup is that the factorized cross section in
SCET+ differs from the ones in SCETI and SCETII only by a subset of the power
corrections it receives relative to the full QCD result,

$\displaystyle\frac{\mathrm{d}\sigma_{\mathrm{I}}}{\mathrm{d}Q\,\mathrm{d}Y\,\mathrm{d}q_{T}\,\mathrm{d}{\mathcal{T}}}$
$\displaystyle=\frac{\mathrm{d}\sigma_{+}}{\mathrm{d}Q\,\mathrm{d}Y\,\mathrm{d}q_{T}\,\mathrm{d}{\mathcal{T}}}\,\Bigl{[}1+\mathcal{O}\Bigl{(}\frac{q_{T}^{2}}{{\mathcal{T}}Q}\Bigr{)}\Bigr{]}\,,$
$\displaystyle\frac{\mathrm{d}\sigma_{\mathrm{II}}}{\mathrm{d}Q\,\mathrm{d}Y\,\mathrm{d}q_{T}\,\mathrm{d}{\mathcal{T}}}$
$\displaystyle=\frac{\mathrm{d}\sigma_{+}}{\mathrm{d}Q\,\mathrm{d}Y\,\mathrm{d}q_{T}\,\mathrm{d}{\mathcal{T}}}\,\Bigl{[}1+\mathcal{O}\Bigl{(}\frac{{\mathcal{T}}^{2}}{q_{T}^{2}}\Bigr{)}\Bigr{]}\,.$
(3.1)

This is illustrated in figure -1014, and follows from comparing eq. (2.40) to
eq. (2.7) and eq. (2.21), respectively. Crucially, eq. (3.1) also holds when
the cross sections are evaluated at common (but not necessarily fixed-order)
scales.

Figure -1014: Venn diagram of power corrections to the factorized double-
differential spectrum. SCETI (blue) and SCETII (green) each capture a set of
power corrections that is expanded away in the SCET+ factorization (red) and
the opposite boundary regime. A third class of power corrections to the
overall soft-collinear limit is captured by the fixed-order calculation in
full QCD (gray).

For example, both $\sigma_{\mathrm{I}}$ and $\sigma_{+}$ share a logarithmic
singularity with respect to ${\mathcal{T}}/Q$, which can be resummed by
running between the scales of the hard, soft, and (refactorized) beam
functions. In SCET+ this amounts to setting the $\mu^{+}$ scales to be equal
to their $\mu^{\mathrm{I}}$ counterparts,

$\mu_{B}^{+}=\mu_{\mathcal{S}}^{+}=\mu_{B}^{\mathrm{I}}\,,\qquad\nu_{B}^{+}=\nu_{\mathcal{S}}^{+}={\mu_{\mathrm{FO}}}\,,\qquad\mu_{S}^{+}=\mu_{S}^{\mathrm{I}}\,,$
(3.2)

such that any large logarithms inside the refactorized beam function in eq.
(2.41) are treated at fixed order. We write
$\mathrm{d}\sigma_{+}\big{|}_{\mu^{\mathrm{I}}}$ to indicate that
$\mathrm{d}\sigma_{+}$ is evaluated at scales that satisfy eq. (3.2). A
natural way to judge the size of the power corrections in eq. (3.1) then is to
compare $\mathrm{d}\sigma_{+}\big{|}_{\mu^{\mathrm{I}}}$ to
$\mathrm{d}\sigma_{\mathrm{I}}\big{|}_{\mu^{\mathrm{I}}}$, with our choices
for $\mu^{\mathrm{I}}$ as given in section 2.2, i.e., including the whole set
of all-order terms from the ${\mathcal{T}}$ resummation in both of them. This
comparison is shown in figure -1013 for representative choices of fixed
${\mathcal{T}}$ and $q_{T}$ at NNLL. We can clearly read off a power-like
behavior of the difference
$\bigl{[}\mathrm{d}\sigma_{\mathrm{I}}-\mathrm{d}\sigma_{+}\bigr{]}_{\mu^{\mathrm{I}}}$
(dotted green) as either $q_{T}\to 0$ for fixed ${\mathcal{T}}$ (left panel)
or ${\mathcal{T}}\to\infty$ for fixed $q_{T}$ (right panel). This also
provides a nontrivial check on our implementation of $\sigma_{\mathrm{I}}$ and
$\sigma_{+}$. This comparison in figure -1013 is analogous to the usual
procedure of comparing the full-theory result for a cross section with its
singular EFT limit at a common scale $\mu_{\mathrm{FO}}$. Here, SCETI takes on
the role of the full theory, while SCET+ provides the singular limit, and the
comparison is performed at common scales $\mu^{\mathrm{I}}$.

Figure -1013: Singular/nonsingular comparison between SCETI and SCET+ at NNLL
as a function of $q_{T}$ for fixed ${\mathcal{T}}=5\,\mathrm{GeV}$ (left) and
as a function of ${\mathcal{T}}$ for fixed $q_{T}=15\,\mathrm{GeV}$ (right).
The orange solid lines show the full SCETI result including resummation. The
dashed blue lines show the corresponding SCET+ singular limit with _only_
SCETI resummation. The dotted green lines show their difference, corresponding
to the power corrections indicated in eq. (3.1). The thin vertical lines
indicate our choice of transition points $(a_{1},a_{2},a_{3})$ with respect to
the regime parameter $a$ (upper horizontal axis), as discussed in section 3.3.

Similarly, both $\sigma_{\mathrm{II}}$ and $\sigma_{+}$ have a common singular
structure as $q_{T}/Q\to 0$. In this case, resumming the shared logarithmic
terms requires running between the hard, beam, and (refactorized) soft
function. In SCET+ this amounts to setting the $\mu^{+}$ scales to be equal to
their $\mu^{\mathrm{II}}$ counterparts,

$\mu_{\mathcal{S}}^{+}=\mu_{S}^{+}=\mu_{S}^{\mathrm{II}}\,,\qquad\nu_{\mathcal{S}}^{+}=\nu_{S}^{\mathrm{II}}\,,$
(3.3)

which treats the large logarithms in the refactorized double-differential soft
function in eq. (2.4) at fixed order. We denote this choice of scales by
$\mathrm{d}\sigma_{+}\bigr{|}_{\mu^{\mathrm{II}}}$, with scale setting in
$b_{T}$ space and the inverse Fourier transform understood as in section 2.3.
In figure -1012 we compare $\mathrm{d}\sigma_{+}\bigr{|}_{\mu^{\mathrm{II}}}$
to $\mathrm{d}\sigma^{\mathrm{II}}\bigr{|}_{\mu^{\mathrm{II}}}$ at NNLL as a
function of ${\mathcal{T}}$ at fixed $q_{T}$ (left) and vice versa (right). It
is clear that even when evaluated at its intrinsic scales,
$\mathrm{d}\sigma_{\mathrm{II}}\bigr{|}_{\mu^{\mathrm{II}}}$ (solid orange)
exhibits an unresummed singularity as ${\mathcal{T}}/q_{T}\ll 1$, which, as
expected, is captured by $\mathrm{d}\sigma_{+}\bigr{|}_{\mu^{\mathrm{II}}}$
(dashed blue) up to power corrections (dotted green). This check is highly
nontrivial as it involves an additional Fourier transform on both sides of the
comparison. We note that the strong kinematic suppression of the double
spectrum for ${\mathcal{T}}\gtrsim q_{T}$ is correctly captured by SCETII,
where central soft modes resolve the phase-space boundary. In SCET+, soft
modes have too little energy and collinear-soft modes are too forward to
resolve it, leading to large power corrections in this region.

Figure -1012: Singular/nonsingular comparison between SCETII and SCET+ at NNLL
as a function of ${\mathcal{T}}$ for fixed $q_{T}=15\,\mathrm{GeV}$ (left) and
as a function of $q_{T}$ for fixed ${\mathcal{T}}=5\,\mathrm{GeV}$ (right).
The orange solid lines show the full SCETII result including resummation. The
dashed blue lines show the corresponding SCET+ singular limit with _only_
SCETII resummation. The dotted green lines show their difference,
corresponding to the power corrections indicated in eq. (3.1). The thin
vertical lines indicate our choice of transition points $(a_{4},a_{5},a_{6})$
with respect to the regime parameter $a$ (upper horizontal axis), as discussed
in section 3.3.

Figure -1011: Singular/nonsingular comparison between the matched SCET
descriptions and QCD at fixed $\mathcal{O}(\alpha_{s})$ as a function of
$q_{T}$ for ${\mathcal{T}}=q_{T}/2$ (left) and ${\mathcal{T}}=q_{T}$ (right).
The orange solid line shows the fixed-order QCD double spectrum, the dashed
blue the matched SCET result in eq. (3.1), and the dotted green the
difference.

As a final important consequence of figure -1014, the complete infrared
structure of the double-differential spectrum for $q_{T}\ll Q$ and
${\mathcal{T}}\ll Q$, i.e., for any hierarchy between $q_{T}$ and
${\mathcal{T}}$, is described by adding the SCETI and SCETII cross sections
and removing the overlap between the two by subtracting the SCET+ cross
section,

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}Q\,\mathrm{d}Y\,\mathrm{d}q_{T}\,\mathrm{d}{\mathcal{T}}}$
$\displaystyle=\Bigl{[}\frac{\mathrm{d}\sigma_{\mathrm{I}}}{\mathrm{d}Q\,\mathrm{d}Y\,\mathrm{d}q_{T}\,\mathrm{d}{\mathcal{T}}}+\frac{\mathrm{d}\sigma_{\mathrm{II}}}{\mathrm{d}Q\,\mathrm{d}Y\,\mathrm{d}q_{T}\,\mathrm{d}{\mathcal{T}}}-\frac{\mathrm{d}\sigma_{+}}{\mathrm{d}Q\,\mathrm{d}Y\,\mathrm{d}q_{T}\,\mathrm{d}{\mathcal{T}}}\Bigr{]}$
$\displaystyle\quad\times\Bigl{[}1+\mathcal{O}\Bigl{(}\frac{q_{T}^{2}}{Q^{2}},\frac{{\mathcal{T}}}{Q}\Bigr{)}\Bigr{]}\,.$
(3.4)

In figure -1011 we numerically check this relation at fixed
$\mathcal{O}(\alpha_{s})$, which requires setting all scales equal to a common
$\mu_{\mathrm{FO}}$. We plot the comparison as a function of $q_{T}$ along
lines of fixed ${\mathcal{T}}/q_{T}=1/2$ (left) and ${\mathcal{T}}/q_{T}=1$
(right), finding excellent agreement between the full result (solid orange)
and the first line on the right-hand side of eq. (3.1) (dashed blue), as
evident from the power-like behavior of their difference (dotted green) as
$q_{T},{\mathcal{T}}\to 0$.

This singular/nonsingular comparison is qualitatively different from the
structure of power corrections in either SCETI or SCETII alone, which we
already verified in figure -1017 and figure -1016. Because SCETI and SCETII
both involve an additional expansion about a specific hierarchy between
$q_{T}$ and ${\mathcal{T}}$, they incur power corrections
$\mathcal{O}({\mathcal{T}}^{2}/q_{T}^{2})$ or
$\mathcal{O}(q_{T}^{2}/(Q{\mathcal{T}}))$, respectively. Accordingly, they
only recover the singular limit of full QCD when approaching it along specific
lines in the $(q_{T},{\mathcal{T}})$ plane. This is different from figure
-1011, where the combined expression in eq. (3.1) (dashed blue) describes the
singular limit $q_{T},{\mathcal{T}}\to 0$ along an arbitrary line of approach,
with the ratio $q_{T}/{\mathcal{T}}$ effectively controlling the “admixture”
of power corrections $\mathcal{O}(q_{T}^{2}/Q^{2})$ and
$\mathcal{O}({\mathcal{T}}/Q)$, respectively. We have verified that also for
other fixed ratios of $q_{T}$ and ${\mathcal{T}}$, the singular behavior of
full QCD is correctly described.

As a final remark, as noted in ref. [90], this property actually qualifies the
expression
$\mathrm{d}\sigma_{\mathrm{I}}+\mathrm{d}\sigma_{\mathrm{II}}-\mathrm{d}\sigma_{+}$
for use as a double-differential subtraction term to treat infrared
divergences in fixed-order calculations.

###  3.2 Matching formula

The structure of power corrections discussed in the previous section, together
with the all-order resummation shared between SCET+ and SCETI or SCETII,
suggests the following matching formula to describe all regions of the double-
differential spectrum:

$\displaystyle\mathrm{d}\sigma^{\mathrm{match}}=\mathrm{d}\sigma_{+}\bigr{|}_{\mu^{+}}$
$\displaystyle+\bigl{[}\mathrm{d}\sigma_{\mathrm{I}}-\mathrm{d}\sigma_{+}\bigr{]}_{\mu^{\mathrm{I}}}+\bigl{[}\mathrm{d}\sigma_{\mathrm{II}}-\mathrm{d}\sigma_{+}\bigr{]}_{\mu^{\mathrm{II}}}$
$\displaystyle+\bigl{[}\mathrm{d}\sigma_{\mathrm{FO}}-\mathrm{d}\sigma_{\mathrm{I}}-\mathrm{d}\sigma_{\mathrm{II}}+\mathrm{d}\sigma_{+}\bigr{]}_{\mu_{\mathrm{FO}}}\,.$
(3.5)

The only ingredient in this matching formula we have not yet discussed is
$\mathrm{d}\sigma_{+}\bigr{|}_{\mu^{+}}$, for which all ingredients in the
SCET+ factorization are evaluated at the SCET+ scales $\mu^{+}$, such that the
full RGE of SCET+ is in effect. In the following we describe the requirements
on $\mu^{+}$ to ensure the best possible prediction across phase space. Our
precise construction of $\mu^{+}$ to satisfy all requirements is the subject
of section 3.3.

In the simplest case, i.e., when the power corrections in eq. (2.40) are
small, and thus the SCET+ parametric assumptions are satisfied, $\mu^{+}$ is
given by the canonical SCET+ scales in eq. (2.4). As for $\mu^{\mathrm{II}}$,
these scales are set in $b_{T}$ space, followed by an inverse Fourier
transform.

As we approach the SCETI region, the resummation inside the refactorization of
the beam function in eq. (2.41) must be turned off,

$\left.\begin{aligned}
\mu_{B}^{+}(q_{T},{\mathcal{T}},b_{T})&\to\mu_{B}^{\mathrm{I}}({\mathcal{T}})\\\
\mu_{\mathcal{S}}^{+}(q_{T},{\mathcal{T}},b_{T})&\to\mu_{B}^{\mathrm{I}}({\mathcal{T}})\\\
\nu_{\mathcal{S}}^{+}(q_{T},{\mathcal{T}},b_{T})&\to\nu_{B}^{+}(q_{T},{\mathcal{T}},b_{T})\end{aligned}\quad\right\\}\quad\text{for}\quad
q_{T}\to\sqrt{Q{\mathcal{T}}}\,.$ (3.6)

In addition we can identify the soft scales in SCETI and SCET+ because the
soft functions are identical,

$\mu_{S}^{+}(q_{T},{\mathcal{T}},b_{T})\to\mu_{S}^{\mathrm{I}}({\mathcal{T}})\quad\text{for}\quad
q_{T}\to\sqrt{Q{\mathcal{T}}}\,.$ (3.7)

These relations must hold for every value of the $b_{T}$ argument of the
scale.

Similarly, as we approach the SCETII region, the scales inside the
refactorized soft function eq. (2.4) must become equal

$\left.\begin{aligned}
\mu_{S}^{+}(q_{T},{\mathcal{T}},b_{T})&\to\mu_{S}^{\mathrm{II}}(q_{T},b_{T})\\\
\mu_{\mathcal{S}}^{+}(q_{T},{\mathcal{T}},b_{T})&\to\mu_{S}^{\mathrm{II}}(q_{T},b_{T})\\\
\nu_{\mathcal{S}}^{+}(q_{T},{\mathcal{T}},b_{T})&\to\nu_{S}^{\mathrm{II}}(q_{T},b_{T})\end{aligned}\quad\right\\}\quad\text{for}\quad
q_{T}\to{\mathcal{T}}\,,$ (3.8)

and we can identify the scales of the common beam function in SCETII and
SCET+,

$\left.\begin{aligned}
\mu_{B}^{+}(q_{T},{\mathcal{T}},b_{T})&\to\mu_{B}^{\mathrm{II}}(q_{T},b_{T})\\\
\nu_{B}^{+}(q_{T},{\mathcal{T}},b_{T})&\to\nu_{B}^{\mathrm{II}}(q_{T},b_{T})\end{aligned}\quad\right\\}\quad\text{for}\quad
q_{T}\to{\mathcal{T}}\,.$ (3.9)

Some of the above requirements for the behavior at the boundary are already
satisfied by the canonical SCET+ scales, e.g., the canonical soft scales in
SCET+ and SCETI are simply equal. The challenge in these cases is to extend
the scale choice onto the opposite boundary, where they are constrained in a
nontrivial way. The nontrivial all-order information in SCET+ is mostly
encoded in the canonical choice of

$\nu_{\mathcal{S}}^{+}(q_{T},{\mathcal{T}},b_{T})=\frac{(b_{0}/b_{T})^{2}}{{\mathcal{T}}}\quad\text{for}\quad{\mathcal{T}}\ll
q_{T}\ll\sqrt{Q{\mathcal{T}}}\,,$ (3.10)

which does not coincide with any scale on either boundary.

It is instructive to explicitly consider the behavior of eq. (3.2) on the
SCETI and SCETII phase-space boundaries, as well as in the fixed-order region.
By construction, for any choice of $\mu^{+}$ scales satisfying eqs. (3.6) and
(3.7) we have

$\mathrm{d}\sigma_{+}\bigr{|}_{\mu^{+}}\to\mathrm{d}\sigma_{+}\bigr{|}_{\mu^{\mathrm{I}}}\quad\text{for}\quad
q_{T}\to\sqrt{Q{\mathcal{T}}}\,.$ (3.11)

It follows that

$\displaystyle\mathrm{d}\sigma^{\mathrm{match}}\to\mathrm{d}\sigma_{\mathrm{I}}\bigr{|}_{\mu^{\mathrm{I}}}$
$\displaystyle+\bigl{[}\mathrm{d}\sigma_{\mathrm{FO}}-\mathrm{d}\sigma_{\mathrm{I}}\bigr{]}_{\mu_{\mathrm{FO}}}$
$\displaystyle+\bigl{[}\mathrm{d}\sigma_{\mathrm{II}}-\mathrm{d}\sigma_{+}\bigr{]}_{\mu^{\mathrm{II}}}-\bigl{[}\mathrm{d}\sigma_{\mathrm{II}}-\mathrm{d}\sigma_{+}\bigr{]}_{\mu_{\mathrm{FO}}}\quad\text{for}\quad
q_{T}\to\sqrt{Q{\mathcal{T}}}\,.$ (3.12)

This mostly coincides with the result in eq. (2.11) of matching
$\mathrm{d}\sigma_{\mathrm{I}}$ to the fixed-order result
$\mathrm{d}\sigma_{\mathrm{FO}}$, and is guaranteed to capture all large
logarithms of ${\mathcal{T}}/Q$ captured by the SCETI RGE. It improves over
eq. (2.11) by also resumming logarithms of $q_{T}/Q$ in the power corrections
$\mathcal{O}({\mathcal{T}}^{2}/q_{T}^{2})$, encoded in
$\bigl{[}\sigma_{\mathrm{II}}-\sigma_{+}\bigr{]}_{\mu^{\mathrm{II}}}$. This is
not a numerically large effect and cannot be exploited to achieve the
resummation of ${\mathcal{T}}$ at next-to-leading power, as it is only a
subset of all power corrections.

Similarly, eqs. (3.8) and (3.9) imply that

$\mathrm{d}\sigma_{+}\bigr{|}_{\mu^{+}}\to\mathrm{d}\sigma_{+}\bigr{|}_{\mu^{\mathrm{II}}}\quad\text{for}\quad
q_{T}\to{\mathcal{T}}\,,$ (3.13)

and consequently

$\displaystyle\mathrm{d}\sigma^{\mathrm{match}}\to\mathrm{d}\sigma_{\mathrm{II}}\bigr{|}_{\mu^{\mathrm{II}}}$
$\displaystyle+\bigl{[}\mathrm{d}\sigma_{\mathrm{FO}}-\mathrm{d}\sigma_{\mathrm{II}}\bigr{]}_{\mu_{\mathrm{FO}}}$
$\displaystyle+\bigl{[}\mathrm{d}\sigma_{\mathrm{I}}-\mathrm{d}\sigma_{+}\bigr{]}_{\mu^{\mathrm{I}}}-\bigl{[}\mathrm{d}\sigma_{\mathrm{I}}-\mathrm{d}\sigma_{+}\bigr{]}_{\mu_{\mathrm{FO}}}\quad\text{for}\quad
q_{T}\to{\mathcal{T}}\,.$ (3.14)

This mostly coincides with the result in eq. (2.29) of matching
$\mathrm{d}\sigma_{\mathrm{II}}$ to the fixed-order result
$\mathrm{d}\sigma_{\mathrm{FO}}$, and thus captures all large logarithms of
$q_{T}/Q$ captured by the SCETII RGE. In addition, it resums logarithms of
${\mathcal{T}}/Q$ in the $\mathcal{O}(q_{T}^{2}/({\mathcal{T}}Q))$ power
corrections encoded in
$[\mathrm{d}\sigma_{\mathrm{I}}-\mathrm{d}\sigma_{+}]_{\mu^{\mathrm{I}}}$.

Finally, in the fixed-order region, all $\mu^{+}$, $\mu^{\mathrm{I}}$, and
$\mu^{\mathrm{II}}$ scales become equal to $\mu_{\mathrm{FO}}$. Thus as
desired, the matched prediction reduces to the fixed-order result,

$\mathrm{d}\sigma^{\mathrm{match}}\to\mathrm{d}\sigma_{\mathrm{FO}}\bigr{|}_{\mu_{\mathrm{FO}}}\quad\text{for}\quad
q_{T},{\mathcal{T}}\to Q\,.$ (3.15)

###  3.3 Profile scales

In this section, we describe our choice of the central $\mu^{+}$ scales for
the various ingredients in the SCET+ factorized cross section, taking into
account the transition to the SCETI and SCETII boundary theories as well as
the transition to the fixed-order region. The SCET+ scales are obtained using
a regime parameter that selects the appropriate combination of scales from the
boundary theories in each region of phase space, and selects a third,
independent choice in the SCET+ “bulk” when necessary. The profile functions
that handle the transition to fixed order can conveniently be reused from
SCETI and SCETII.

Scale | SCETI | SCET+ | SCETII
---|---|---|---
$\mu_{H}$ | $Q$ | $Q$ | $Q$
$\mu_{B}$ | $\sqrt{{\mathcal{T}}Q}$ | $b_{0}/b_{T}$ | $b_{0}/b_{T}$
$\nu_{B}$ |  | $Q$ | $Q$
$\mu_{\mathcal{S}}$ |  | $b_{0}/b_{T}$ |
$\nu_{\mathcal{S}}$ |  | $(b_{0}/b_{T})^{2}/{\mathcal{T}}$ |
$\mu_{S}$ | ${\mathcal{T}}$ | ${\mathcal{T}}$ | $b_{0}/b_{T}$
$\nu_{S}$ |  |  | $b_{0}/b_{T}$
Table 1: Summary of canonical scales in SCETI, SCET+, and SCETII [see eq.
(2.10), (2.3) and (2.4)]. For SCET+ and SCETII we give the canonical scales in
$b_{T}$ space.

We start by summarizing the canonical scales for SCETI, SCETII, SCET+ in table
1. At these scales, the arguments of logarithms in the ingredients of the
factorized cross section are order one, i.e., all large logarithms are
resummed by RG evolution. To interpolate between the canonical scales in
different regimes, we find it convenient to introduce the regime parameter

$\displaystyle
a=3-\frac{\lvert\ln({\mathcal{T}}/Q)\rvert}{\lvert\ln(q_{T}/Q)\rvert}\,.$
(3.16)

Its definition is carefully chosen such that $a=1$ when the SCETI parametric
relation is exactly satisfied, $q_{T}=\sqrt{{\mathcal{T}}Q}$, and $a=2$ on the
SCETII boundary of phase space, $q_{T}={\mathcal{T}}$. As illustrated in the
left panel of figure -1010, the canonical SCET+ region lies at intermediate
$a\sim 1.5$. The requirements on the SCET+ scales were given in eqs. (3.6) and
(3.7) for the transition to SCETI, and in eqs. (3.8) and (3.9) for SCETII. To
satisfy these requirements, we take weighted products of scales on the
boundary and in the bulk, schematically,

$\displaystyle\mu^{+}=\bigl{[}\mu^{\mathrm{I}}\bigr{]}^{h_{\mathrm{I}}(a)}\,\bigl{[}\mu^{+}_{\mathrm{bulk}}\bigr{]}^{h_{+}(a)}\bigl{[}\mu^{\mathrm{II}}\bigr{]}^{h_{\mathrm{II}}(a)}\,.$
(3.17)

The weights in the exponent are given by helper functions that depend on $a$,
as illustrated in the right panel of figure -1010. They satisfy

$h_{\mathrm{I}}(a)+h_{+}(a)+h_{\mathrm{II}}(a)=1\,,$ (3.18)

for any $a$ and are given explicitly in eq. (3.3) below. The helper functions
ensure that the appropriate scales are used in each region, e.g.,
$h_{\mathrm{II}}(a)$ is one in the vicinity of $a=2$ and vanishes for $a<1.5$,
with a smooth transition between regions.

Figure -1010: Left: Illustration of the regime parameter $a$ that governs the
matching between EFTs. We show lines of constant $a$ in the
$(q_{T},{\mathcal{T}})$ plane. For $a=1$ the SCETI parametric relation is
exactly satisfied, $q_{T}=\sqrt{Q{\mathcal{T}}}$, whereas for $a=2$, the
SCETII parametric relation is exactly satisfied, $q_{T}={\mathcal{T}}$. Right:
Helper functions used to interpolate between scales on the boundaries (SCETI,
SCETII) and in the bulk (SCET+). The helper functions have continuous
derivatives and always sum to one. The individual helper functions are exactly
one in their respective canonical regions. See the text for a detailed
discussion.

For the soft and collinear-soft scales, eq. (3.17) takes the following
concrete form:

$\displaystyle\mu_{\mathcal{S}}^{+}$
$\displaystyle=\bigl{[}\mu_{B}^{\mathrm{I}}\bigr{]}^{h_{\mathrm{I}}(a)}\,$
$\displaystyle\bigl{[}\mu_{\mathcal{S},\mathrm{bulk}}^{+}\bigr{]}^{h_{+}(a)}\,$
$\displaystyle\bigl{[}\mu_{S}^{\mathrm{II}}\bigr{]}^{h_{\mathrm{II}}(a)}\,,$
$\displaystyle\nu_{\mathcal{S}}^{+}$
$\displaystyle=\bigl{[}\nu\bigr{]}^{h_{\mathrm{I}}(a)}\,$
$\displaystyle\bigl{[}\nu_{\mathcal{S},\mathrm{bulk}}^{+}\bigr{]}^{h_{+}(a)}\,$
$\displaystyle\bigl{[}\nu_{S}^{\mathrm{II}}\bigr{]}^{h_{\mathrm{II}}(a)}\,\,,$
$\displaystyle\mu_{S}^{+}$
$\displaystyle=\bigl{[}\mu_{S}^{\mathrm{I}}\bigr{]}^{h_{\mathrm{I}}(a)}\,$
$\displaystyle\bigl{[}\mu_{S,\mathrm{bulk}}^{+}\bigr{]}^{h_{+}(a)}\,$
$\displaystyle\bigl{[}\mu_{S}^{\mathrm{II}}\bigr{]}^{h_{\mathrm{II}}(a)}\,\,.$
(3.19)

The most nontrivial of these cases is $\nu_{\mathcal{S}}$, which must be equal
to the overall $\nu$ in the SCETI region to turn off the rapidity resummation
there, has a distinct canonical value in the SCET+ bulk, and must asymptote to
yet another value on the SCETII boundary. We note that $\mu_{S}^{+}$ also
requires a distinct treatment in the bulk to ensure that the hierarchy
$\mu_{S}^{+}<\mu_{\mathcal{S}}^{+}$ inside the refactorized soft function, as
implied by the SCET+ power counting, is not upset by variations (see next
subsection). Our central choices for the above scales in the bulk are

$\displaystyle\mu_{\mathcal{S},\mathrm{bulk}}^{+}$
$\displaystyle=\mu_{\mathrm{FO}}\,f_{\mathrm{run}}^{\mathrm{II}}\Bigl{(}\frac{q_{T}}{Q},\frac{b_{0}}{b^{\ast}(b_{T})\,Q}\Bigr{)}\,,\quad\nu_{\mathcal{S},\mathrm{bulk}}^{+}=\mu_{\mathrm{FO}}\frac{\Bigl{[}f_{\mathrm{run}}^{\mathrm{II}}\Bigl{(}\frac{q_{T}}{Q},\frac{b_{0}}{b^{\ast}(b_{T})Q}\Bigr{)}\Bigr{]}^{2}}{f_{\mathrm{run}}^{\mathrm{I}}\Bigl{(}\frac{{\mathcal{T}}}{Q}\Bigr{)}}\,,$
$\displaystyle\mu_{S,\mathrm{bulk}}^{+}$
$\displaystyle=\mu_{\mathrm{FO}}\,f_{\mathrm{run}}^{\mathrm{I}}\Bigl{(}\frac{{\mathcal{T}}}{Q}\Bigr{)}\,.$
(3.20)

The profile function $f_{\mathrm{run}}^{\mathrm{I}}$ was introduced for the
transition between SCETI and fixed-order QCD in eq. (2.14), and similarly for
the hybrid profile $f_{\mathrm{run}}^{\mathrm{II}}$ in eq. (2.32) and the
nonperturbative $b^{\ast}(b_{T})$ prescription in eq. (2.34). These functions
turn off the resummation of logarithms involving $q_{T}$ ($b_{T}$) and
${\mathcal{T}}$, respectively, as the fixed-order regime is approached, and
also ensure that scales are frozen in the nonperturbative regime to avoid the
Landau pole. It is straightforward to check that away from the nonperturbative
region, the above bulk scales all assume their canonical values for
$q_{T},{\mathcal{T}}\ll Q$ as given in table 1, and asymptote to
$\mu_{\mathrm{FO}}$ when simultaneously taking $q_{T},{\mathcal{T}}\to Q$. The
beam function scales in the bulk can simply be associated with their SCETII
counterparts and only require a transition towards the SCETI boundary,

$\displaystyle\mu_{B}^{+}$
$\displaystyle=\bigl{[}\mu_{B}^{\mathrm{I}}\bigr{]}^{h_{\mathrm{I}}(a)}\,\bigl{[}\mu_{B}^{\mathrm{II}}\bigr{]}^{h_{+}(a)+h_{\mathrm{II}}(a)}\,,$
$\displaystyle\nu_{B}^{+}$
$\displaystyle=\bigl{[}\nu\bigr{]}^{h_{\mathrm{I}}(a)}\,\bigl{[}\nu_{B}^{\mathrm{II}}\bigr{]}^{h_{+}(a)+h_{\mathrm{II}}(a)}\,.$
(3.21)

In our numerical implementation, we choose the helper functions
$h_{\mathrm{I},\mathrm{II},+}$ as

$\displaystyle h_{\mathrm{I}}(a)$
$\displaystyle\equiv\begin{cases}1&a<a_{1}\,,\\\ 1-c_{123}(a)&a_{1}\leq
a<a_{2}\,,\\\ c_{312}(a)&a_{2}\leq a<a_{3}\,,\\\ 0&a_{3}\leq
a\,,\end{cases}\qquad h_{\mathrm{II}}(a)\equiv\begin{cases}0&a<a_{4}\,,\\\
c_{456}(a)&a_{4}\leq a<a_{5}\,,\\\ 1-c_{645}(a)&a_{5}\leq a<a_{6}\,,\\\
1&a_{6}\leq a\,,\end{cases}$ $\displaystyle h_{+}(a)$ $\displaystyle\equiv
1-h_{\mathrm{I}}(a)-h_{\mathrm{II}}(a)\,,$ (3.22)

where the polynomials governing the interpolation between zero and one are

$\displaystyle c_{ijk}(a)=\frac{(a-a_{i})^{2}}{(a_{i}-a_{j})(a_{i}-a_{k})}\,.$
(3.23)

The transition points $a_{1,\ldots,6}$ determine the transition between the
different regions, as can be seen from the helper functions in figure -1010:
For values $a_{3}\leq a<a_{4}$, the exact canonical SCET+ scales are selected,
implying that the resummation of logarithms of both $q_{T}$ and
${\mathcal{T}}$ is fully turned on. For lower values $a_{1}\leq a<a_{3}$, the
additional $q_{T}$ resummation is smoothly turned off and for $a<a_{1}$, SCETI
scales are used so that only logarithms of ${\mathcal{T}}$ are resummed.
Conversely, for higher values of the regime parameter $a_{4}\leq a<a_{6}$, the
resummation of ${\mathcal{T}}$ logarithms is smoothly turned off. At values
$a_{6}\leq a$, SCETII scales are selected by the helper functions, and the
additional resummation of logarithms of ${\mathcal{T}}$ is completely turned
off.

In practice we use
$(a_{1},a_{2},a_{3},a_{4},a_{5},a_{6})=(0.5,\,1.0,\,1.5,\,1.5,\,1.75,\,2.0)$.
This choice ensures that for $a\geq a_{6}=2$, we fully recover SCETII
resummation and faithfully describe the kinematic edge at
$q_{T}\sim{\mathcal{T}}$ by preserving the $\mathcal{O}(1)$ cancellation
between $\sigma_{+}\bigr{|}_{\mu^{\mathrm{II}}}$ and the SCETII nonsingular
contribution visible at $a\sim 2$ in the left panel of figure -1012. (In both
figures -1013 and -1012, corresponding values of $a$ are indicated on the
horizontal axis at the top of the panels.) On the other hand, from figure
-1013 we observe that power corrections from SCETI are smaller and tend to set
in at values of $a$ lower than the naively expected $a=1$. E.g., an
$\mathcal{O}(1)$ cancellation between $\sigma_{+}\bigr{|}_{\mu^{\mathrm{I}}}$
and the SCETI nonsingular only is in effect around $a\sim 0.5$ in the right
panel of figure -1013, leaving more room for slowly turning off the SCET+
resummation down towards $a_{1}=0.5$. This is expected because the SCETI
nonsingular encodes the suppression of collinear recoil beyond the naive
phase-space boundary at $a\sim 1$ ($q_{T}\sim\sqrt{Q{\mathcal{T}}}$) that is
washed out by the PDFs, unlike the sharp kinematic edge at
$q_{T}\sim{\mathcal{T}}$ encoded in the SCETII nonsingular. For simplicity we
set $a_{3}=a_{4}$ for our central prediction, i.e., we shrink the canonical
SCET+ region to a point at $a=1.5$, and fix $a_{2}$ ($a_{5}$) to be the
midpoint between $a_{1}$ and $a_{3}$ ($a_{4}$ and $a_{6}$). Variations of the
transition points, including independent variations of $a_{3}$ and $a_{4}$,
are considered as part of our uncertainty estimate described in the next
section.

###  3.4 Perturbative uncertainties

In this section we describe how we assess perturbative uncertainties by
varying the scales entering the matched prediction in eq. (3.2). Following the
same approach as for SCETI and SCETII on their own (see sections 2.2 and 2.3),
we distinguish between resummation uncertainties and a fixed-order
uncertainty. The fixed-order uncertainty $\Delta_{\mathrm{FO}}$ is estimated
by varying $\mu_{\mathrm{FO}}$ up and down by a factor of two, i.e., by
setting ${\mu_{\mathrm{FO}}}=\\{Q/2,2Q\\}$. Since all scales (in any piece of
the matching formula) include an overall factor of $\mu_{\mathrm{FO}}$, the
ratios between the various scales remain unchanged and the same logarithms are
resummed. The fixed-order uncertainty $\Delta_{\mathrm{FO}}$ is then taken to
be the maximum deviation from the central cross section.

We consider several sources of resummation uncertainty entering the matched
prediction in eq. (3.2). To probe the tower of logarithms of ${\mathcal{T}}/Q$
predicted by the SCETI RGE, we perform variations of $\mu_{B}^{\mathrm{I}}$
and $\mu_{S}^{\mathrm{I}}$ parametrized by $\alpha$ and $\beta$ as in eq.
(2.2). This directly affects the resummed power corrections
$\bigl{[}\mathrm{d}\sigma_{\mathrm{I}}-\mathrm{d}\sigma_{+}\bigr{]}_{\mu_{\mathrm{I}}}$
captured by SCETI. In addition, however,
$\mathrm{d}\sigma_{+}\bigr{|}_{\mu_{+}}$ near the SCETI boundary also
undergoes variations because for large $h_{\mathrm{I}}$, the SCET+ scales in
eqs. (3.3) and (3.3) strongly depend on their SCETI counterparts and inherit
their variations. Our setup thus ensures that in (or near) the SCETI region,
variations probing resummed logarithms of ${\mathcal{T}}/Q$ are properly
treated as correlated between the SCET+ cross section and the SCETI matching
correction. When referring to the matched prediction in eq. (3.2), we take
$\Delta_{\mathrm{I}}$ to be the maximum deviation of
$\mathrm{d}\sigma_{\mathrm{match}}$ from its central value under these
correlated variations of $\alpha$, $\beta$.

In complete analogy, we define $\Delta_{\mathrm{II}}$ as the maximum deviation
under correlated variations of $\mu_{\mathrm{II}}$ as described in section
2.3. These variations act on both
$\bigl{[}\mathrm{d}\sigma_{\mathrm{II}}-\mathrm{d}\sigma_{+}\bigr{]}_{\mu_{\mathrm{II}}}$
and $\mathrm{d}\sigma_{+}\bigr{|}_{\mu_{+}}$, where now the SCET+ scales
inherit variations from $\mu_{\mathrm{II}}$ near the SCETII boundary (where
$h_{\mathrm{II}}$ is large). As a result, $\Delta_{\mathrm{II}}$ probes an
all-order set of logarithms of $(b_{0}/b_{T})/Q$ predicted and resummed by the
SCETII RGE, and properly captures the correlated tower of logarithms in SCET+.
We like to stress that our setup is fully general with respect to the method
chosen to perform scale variations for the boundary theories, as any variation
will automatically be inherited by SCET+.

As a final source of uncertainty, we consider the uncertainty inherent in our
matching procedure and in our choice of SCET+ scales in the bulk. To estimate
this we perform the following 8 variations of the (in principle arbitrary)
transition points $(a_{1},a_{3},a_{4},a_{6})$,

$\displaystyle(\uparrow,-,-,-)\,,\qquad$
$\displaystyle(-,\downarrow,-,-)\,,\qquad$
$\displaystyle(-,-,-,\downarrow)\,,\qquad$
$\displaystyle(-,\uparrow,\uparrow,-)\,,$
$\displaystyle(\downarrow,-,-,-)\,,\qquad$
$\displaystyle(-,-,\uparrow,-)\,,\qquad$
$\displaystyle(-,-,-,\uparrow)\,,\qquad$
$\displaystyle(-,\downarrow,\downarrow,-)\,,$ (3.24)

where $\uparrow\,(\downarrow)$ indicates a variation by $+0.2\,(-0.2)$, a dash
indicates keeping the transition point fixed, and we always maintain
$a_{2}=(a_{1}+a_{3})/2$ and $a_{5}=(a_{4}+a_{6})/2$. In addition, we perform
the following two variations of the SCET+ bulk scales,

$\displaystyle\mu_{\mathcal{S},\mathrm{bulk}}^{+}$
$\displaystyle=\mu_{\mathrm{FO}}\,\Bigl{(}\frac{q_{T}}{{\mathcal{T}}}\Bigr{)}^{+\gamma/2}\,f_{\mathrm{run}}^{\mathrm{II}}\Bigl{(}\frac{q_{T}}{Q},\frac{b_{0}}{b^{\ast}Q}\Bigr{)}\,,$
$\displaystyle\mu_{S,\mathrm{bulk}}^{+}$
$\displaystyle=\mu_{\mathrm{FO}}\,\Bigl{(}\frac{q_{T}}{{\mathcal{T}}}\Bigr{)}^{-\gamma/2}\,f_{\mathrm{run}}^{\mathrm{I}}\Bigl{(}\frac{{\mathcal{T}}}{Q}\Bigr{)}\,,\qquad\qquad\gamma=\\{+1/6,-1/6\\}\,,$
(3.25)

where $\gamma=0$ corresponds to the central scales in eq. (3.3). Similarly to
the role of $\beta$ in the SCETI variations [see eq. (2.2)], making the
strength of the $\gamma$ variations depend on the ratio $q_{T}/{\mathcal{T}}$
ensures that the hierarchy $\mu_{S}<\mu_{\mathcal{S}}$ implied by the SCET+
power counting is not upset by variations, counting $b_{0}/b_{T}\sim q_{T}$.
We note that the third independent bulk scale
$\nu_{\mathcal{S},\mathrm{bulk}}^{+}$ does not require independent variation
because it only enters through rapidity logarithms of
$\nu_{B}^{+}/\nu_{\mathcal{S}}^{+}$, which are already being probed by
variations of $\nu_{B}^{+}$ inherited from SCETII. Taking the envelope of the
eight transition point variations and the two bulk scale variations, we obtain
a third contribution to the resummation uncertainty denoted by $\Delta_{+}$.
The total uncertainty assigned to the matched prediction is then given by
adding all contributions in quadrature,

$\displaystyle\Delta_{\mathrm{total}}=\Delta_{+}\oplus\Delta_{\mathrm{I}}\oplus\Delta_{\mathrm{II}}\oplus\Delta_{\mathrm{FO}}\,.$
(3.26)

###  3.5 Differential and cumulant scale setting

We will now discuss the issue of differential versus cumulant scale setting,
starting with the simpler case of a cross section differential in a single
observable and using 0-jettiness ${\mathcal{T}}$ as an example. There are two
equivalent quantities of interest in this case, namely the spectrum
$\mathrm{d}\sigma/\mathrm{d}{\mathcal{T}}$ with respect to ${\mathcal{T}}$ and
the cumulant cross section $\sigma({\mathcal{T}_{\mathrm{cut}}})$ with a cut
on ${\mathcal{T}}$. The two quantities are related by

$\sigma({\mathcal{T}_{\mathrm{cut}}})=\int_{0}^{{\mathcal{T}_{\mathrm{cut}}}}\\!\mathrm{d}{\mathcal{T}}\,\frac{\mathrm{d}\sigma}{\mathrm{d}{\mathcal{T}}}\,,$
(3.27)

where we suppress the dependence on $Q^{2}$ and $Y$ for the purposes of this
subsection. Accordingly, in a resummation analysis one can implement the
resummation scales either in terms of the differential variable
${\mathcal{T}}$ to directly predict the spectrum, or in terms of the cumulant
variable ${\mathcal{T}}_{\mathrm{cut}}$ to predict the cross section
integrated up to ${\mathcal{T}_{\mathrm{cut}}}$. The other observable then
follows from eq. (3.27).

Explicitly, with differential scale setting (indicated by the subscript), the
differential and cumulant cross section are given by

$\displaystyle\frac{\mathrm{d}\sigma_{\mathrm{diff}}}{\mathrm{d}{\mathcal{T}}}$
$\displaystyle=\frac{\mathrm{d}\sigma}{\mathrm{d}{\mathcal{T}}}\Big{|}_{\mu({\mathcal{T}})}\,,$
$\displaystyle\sigma_{\mathrm{diff}}({\mathcal{T}}_{\mathrm{cut}})$
$\displaystyle=\int^{{\mathcal{T}}_{\mathrm{cut}}}\\!\mathrm{d}{\mathcal{T}}\,\biggl{[}\theta({\mathcal{T}}>{\mathcal{T}}_{\mathrm{np}})\frac{\mathrm{d}\sigma}{\mathrm{d}{\mathcal{T}}}\Big{|}_{\mu({\mathcal{T}})}+\theta({\mathcal{T}}\leq{\mathcal{T}}_{\mathrm{np}})\frac{\mathrm{d}\sigma}{\mathrm{d}{\mathcal{T}}}\Big{|}_{\mu({\mathcal{T}}_{\mathrm{np}})}\biggr{]}\,.$
(3.28)

In the first term under the integral in the cumulant cross section, all scales
$\mu$ entering the resummed and matched prediction depend on the integration
variable ${\mathcal{T}}$. Because our setup only reliably predicts the
spectrum away from the nonperturbative region, we choose to integrate the
resummed spectrum with differential scale setting up from some small cutoff
${\mathcal{T}}_{\mathrm{np}}$, and include an “underflow” contribution given
by the second term under the integral. For the underflow contribution for
${\mathcal{T}}\leq{\mathcal{T}}_{\mathrm{np}}$, the spectrum is evaluated at
fixed scales corresponding to ${\mathcal{T}}_{\mathrm{np}}$, such that the
integral can be done analytically. The underflow contribution is Sudakov
suppressed and thus typically small.

Using cumulant scale setting, we instead use

$\displaystyle\sigma_{\mathrm{cumul}}({\mathcal{T}}_{\mathrm{cut}})$
$\displaystyle=\int^{{\mathcal{T}}_{\mathrm{cut}}}\\!\mathrm{d}{\mathcal{T}}\,\frac{\mathrm{d}\sigma}{\mathrm{d}{\mathcal{T}}}\Big{|}_{\mu({\mathcal{T}}_{\mathrm{cut}})}\,,$
$\displaystyle\frac{\mathrm{d}\sigma_{\mathrm{cumul}}}{\mathrm{d}{\mathcal{T}}}$
$\displaystyle=\frac{\mathrm{d}\sigma}{\mathrm{d}{\mathcal{T}}}\Big{|}_{\mu({\mathcal{T}})}+\sum_{i}\bigg{[}\frac{\mathrm{d}}{\mathrm{d}\mu_{i}}\int^{{\mathcal{T}}}\\!\mathrm{d}{\mathcal{T}}^{\prime}\,\frac{\mathrm{d}\sigma}{\mathrm{d}{\mathcal{T}}^{\prime}}\bigg{]}_{\mu({\mathcal{T}})}\,\frac{\mathrm{d}\mu_{i}({\mathcal{T}})}{\mathrm{d}{\mathcal{T}}}\,.$
(3.29)

In this case, the scales in the cumulant cross section depend on
${\mathcal{T}_{\mathrm{cut}}}$, and not the integration variable
${\mathcal{T}}$, so the integral up to ${\mathcal{T}_{\mathrm{cut}}}$ can
easily be performed analytically. The expression for the differential cross
section arises from taking the derivative of the cumulant cross section, where
the chain rule leads to the sum of derivatives of each of the scales $\mu_{i}$
in $\mu$ with respect to ${\mathcal{T}}$.

Cumulant scale setting ensures that for ${\mathcal{T}}_{\rm cut}\to Q$, the
resummed and matched cumulant cross section exactly reproduces the inclusive
fixed-order cross section. This follows from the generic requirement on
profile scales in the fixed-order region,

$\displaystyle\mu_{i}({\mathcal{T}_{\mathrm{cut}}})\to\mu_{\mathrm{FO}}\quad\text{for}\quad{\mathcal{T}_{\mathrm{cut}}}\to
Q\,.$ (3.30)

Thus for cumulant scale setting, the spectrum has the correct (fixed-order)
normalization. However, the additional derivatives of the scales in eq. (3.5)
tend to produce artifacts in the spectrum if the profile functions
$\mu_{i}({\mathcal{T}})$ used to interpolate between the resummation region
${\mathcal{T}}\ll Q$ to the fixed-order region ${\mathcal{T}}\sim Q$ undergo a
rapid transition. In particular, a smooth matching to the fixed-order
prediction at the level of the differential spectrum typically requires
differential scale setting. Moreover, the scale variations using cumulant
scale setting tend to produce unreliable uncertainties for the spectrum.

If instead differential scale setting is used, the spectrum is free from such
artifacts. However, the integral of the spectrum will not exactly recover the
inclusive fixed-order cross section, and the uncertainties obtained for the
cumulant by integrating the spectrum scale variations tend to accumulate and
end up being much larger than they should be for the total cross section. As
in the case of the spectrum with cumulant scale setting, this mismatch purely
arises from residual scale dependence, and therefore is formally beyond the
working order. It can however still be numerically significant.

Therefore, in general one should use the scale setting that is appropriate for
the quantity of interest, i.e., one should use cumulant scale setting when
making predictions for the cumulant, and differential scale setting when one
is interested in the spectrum. This issue of differential versus cumulant
scale setting is well appreciated in the literature for the single-
differential case, see e.g. refs. [65, 91, 50, 92]. It fundamentally results
from the fact that long-range correlations across the spectrum are not
accounted for by the profile scales used for the differential predictions.
Conversely, profile scales for the cumulant do not correctly capture the slope
of the cumulant and its uncertainty. An elaborate procedure for obtaining a
spectrum with differential scales that still produce the exact cross section
and uncertainties was developed in ref. [92]. In the Geneva Monte Carlo
generator, the mismatch between differential and cumulant scales is accounted
for by adding explicit higher-order terms [50].

For a simultaneous measurement of $q_{T}$ and ${\mathcal{T}}$, there are in
principle four quantities of interest, namely the double-differential spectrum
$\mathrm{d}\sigma/\mathrm{d}q_{T}\,\mathrm{d}{\mathcal{T}}$, the single-
differential spectra
$\mathrm{d}\sigma({q_{T}^{\mathrm{cut}}})/\mathrm{d}{\mathcal{T}}$ and
$\mathrm{d}\sigma({\mathcal{T}_{\mathrm{cut}}})/\mathrm{d}q_{T}$ with a cut on
the other variable, and the double cumulant
$\sigma({q_{T}^{\mathrm{cut}}},{\mathcal{T}_{\mathrm{cut}}})$. They are all
related by integration or differentiation, allowing for four different ways of
setting scales in each case. For our explicit numerical results in section 4,
we take a pragmatic approach and use the appropriate combination of
differential or cumulant scale setting with respect to either $q_{T}$ or
${\mathcal{T}}$ for each of these quantities. This is achieved by evaluating
the resummed prediction at profile scales given by the setup described in
sections 2.2 and 2.3 as well as section 3.3, but with $q_{T}$
(${\mathcal{T}}$) replaced by ${q_{T}^{\mathrm{cut}}}$
(${\mathcal{T}_{\mathrm{cut}}}$) as appropriate. In this way we are guaranteed
to avoid artifacts from profile functions in spectrum observables, and on the
other hand ensure that cumulant observables have the correct limiting
behavior; e.g., $\sigma({q_{T}^{\mathrm{cut}}},{\mathcal{T}_{\mathrm{cut}}})$
will by construction recover the inclusive fixed-order cross section when
lifting both cuts, while
$\mathrm{d}\sigma({q_{T}^{\mathrm{cut}}})/\mathrm{d}{\mathcal{T}}$ and
$\mathrm{d}\sigma({\mathcal{T}_{\mathrm{cut}}})/\mathrm{d}q_{T}$ exactly
recover the resummed and matched prediction for the respective inclusive
spectrum at large values of the cut.

Nevertheless, it is interesting to ask how well the different combinations of
differential and cumulant scale setting fare for observables other than the
one they are designed to describe. In particular we should ask how well the
$(q_{T},{\mathcal{T}})$ scale setting we described in earlier sections
performs at the level of cumulant observables and their inclusive limit. To do
so, we can always promote a spectrum using differential scale setting in
$q_{T}$ (${\mathcal{T}}$) to a prediction for the cumulant up to
${q_{T}^{\mathrm{cut}}}$ (${\mathcal{T}_{\mathrm{cut}}}$) using the analogue
of eq. (3.5). The only nontrivial new procedure is computing the double
cumulant directly from $(q_{T},{\mathcal{T}})$ scales, where we need to
account for an overlap in underflow contributions as

$\displaystyle\sigma_{\mathrm{diff},\mathrm{diff}}({q_{T}^{\mathrm{cut}}},{\mathcal{T}_{\mathrm{cut}}})$
$\displaystyle=\int^{{q_{T}^{\mathrm{cut}}}}\\!\mathrm{d}q_{T}\,\int^{{\mathcal{T}}_{\mathrm{cut}}}\\!\mathrm{d}{\mathcal{T}}\,\Bigl{[}\theta(q_{T}>q_{T}^{\mathrm{np}})\,\theta({\mathcal{T}}>{\mathcal{T}}_{\mathrm{np}})\frac{\mathrm{d}\sigma}{\mathrm{d}q_{T}\,\mathrm{d}{\mathcal{T}}}\Big{|}_{\mu(q_{T},{\mathcal{T}})}$
(3.31) $\displaystyle\qquad\qquad\qquad\qquad\quad+\theta(q_{T}\leq
q_{T}^{\mathrm{np}})\,\theta({\mathcal{T}}>{\mathcal{T}}_{\mathrm{np}})\frac{\mathrm{d}\sigma}{\mathrm{d}q_{T}\,\mathrm{d}{\mathcal{T}}}\Big{|}_{\mu(q_{T}^{\mathrm{np}},{\mathcal{T}})}$
$\displaystyle\qquad\qquad\qquad\qquad\quad+\theta(q_{T}>q_{T}^{\mathrm{np}})\,\theta({\mathcal{T}}\leq{\mathcal{T}}_{\mathrm{np}})\frac{\mathrm{d}\sigma}{\mathrm{d}q_{T}\,\mathrm{d}{\mathcal{T}}}\Big{|}_{\mu(q_{T},{\mathcal{T}}_{\mathrm{np}})}$
$\displaystyle\qquad\qquad\qquad\qquad\quad-\theta(q_{T}\leq
q_{T}^{\mathrm{np}})\,\theta({\mathcal{T}}\leq{\mathcal{T}}_{\mathrm{np}})\frac{\mathrm{d}\sigma}{\mathrm{d}q_{T}\,\mathrm{d}{\mathcal{T}}}\Big{|}_{\mu(q_{T}^{\mathrm{np}},{\mathcal{T}}_{\mathrm{np}})}\Bigr{]}\,.$

The distinction between differential or cumulant scale setting is only
relevant for $q_{T}$ versus ${q_{T}^{\mathrm{cut}}}$ but not for the
underlying resummation in $b_{T}$ space, so we suppress the dependence of the
hybrid scales on $b_{T}$. In practice we use
$q_{T}^{\mathrm{np}}={\mathcal{T}}_{\mathrm{np}}=1\,\mathrm{GeV}$, and
implement the integrals in eqs. (3.5) and (3.31) as sums over logarithmically
spaced bins with bin size
$\Delta(\log_{10}q_{T})=\Delta(\log_{10}{\mathcal{T}})=0.08$, where the
spectrum is evaluated at the logarithmic midpoint of the bin. Scale variations
in the integrated results are performed by integrating each instance of the
spectrum separately and computing maximum deviations from central in the end.
The final results are interpolated for clarity.

Figure -1009: The double cumulant cross section as a function of
${\mathcal{T}_{\mathrm{cut}}}$ for ${q_{T}^{\mathrm{cut}}}=100\,\mathrm{GeV}$
(left) and as a function of ${q_{T}^{\mathrm{cut}}}$ for
${\mathcal{T}_{\mathrm{cut}}}=100\,\mathrm{GeV}$ (right). The bands indicate
the total perturbative uncertainty $\Delta_{\mathrm{total}}$, see section 3.4.
The colors correspond to different scale setting prescriptions (default: solid
orange); see the text for details.

In figures -1009 to -1007, we compare our default scale setting for various
cumulant observables (solid orange) to more differential scale setting (dashed
blue and dotted green), i.e., choosing $\mu$ in terms of $q_{T}$ rather than
${q_{T}^{\mathrm{cut}}}$ and/or ${\mathcal{T}}$ rather than
${\mathcal{T}_{\mathrm{cut}}}$. In figure -1009, we show the double cumulant
cross section, for which our default is to use scales in terms of
$q_{T}^{\mathrm{cut}}$ and ${\mathcal{T}}_{\mathrm{cut}}$. The horizontal
reference line indicates the inclusive fixed-order cross section. In figure
-1008 we show the ${\mathcal{T}}$ spectrum with a cut on $q_{T}$, for which
our default scales are in terms of $q_{T}^{\mathrm{cut}}$ and ${\mathcal{T}}$,
and the converse for figure -1007. In figures -1008 and -1007 the left panel
shows the dependence on the cut at a representative point along the spectrum,
with the reference line indicating the resummed prediction for the inclusive
(strictly single-differential) spectrum. The right panel shows the spectrum at
a representative choice of the cut.

We start by observing that in all cases, the predictions obtained using the
default scale setting (solid orange) cleanly asymptote to the respective
target observable (the reference line) for large values of the cut. The
central double-differential prediction in the left panel of figure -1007
slightly overshoots the inclusive result beyond the phase-space boundary
${\mathcal{T}_{\mathrm{cut}}}\gtrsim q_{T}$ (where our calculation is
effectively a leading-order calculation), but is monotonic within
uncertainties. Furthermore, the uncertainty obtained using our default is
smaller than any of the ones obtained from more differential scale setting.
This is expected because differential scale setting cannot account for
correlations between different bins of the spectrum, giving rise to a larger
band in the cumulant cross sections.

We further note that predictions obtained using $q_{T}$ or
${q_{T}^{\mathrm{cut}}}$ scale setting are mutually compatible, i.e., their
uncertainty bands (very nearly) overlap, as long as the scale setting with
respect to ${\mathcal{T}}$ is done the same way in both cases. This can be
seen from the right panel of figure -1009 by contrasting the default
$({q_{T}^{\mathrm{cut}}},{\mathcal{T}_{\mathrm{cut}}})$ scales (solid orange)
and $(q_{T},{\mathcal{T}_{\mathrm{cut}}})$ scales (dotted green). Similarly,
in figure -1008 we find that the default
$({q_{T}^{\mathrm{cut}}},{\mathcal{T}})$ scales (solid orange) and
$(q_{T},{\mathcal{T}})$ scale setting (dashed blue) roughly differ by their
respective uncertainties. In principle these relations are expected since the
unphysical scale dependence is canceled by higher-order corrections, which our
scale variations are designed to probe. For the case of $q_{T}$ versus
${q_{T}^{\mathrm{cut}}}$ scales in particular, we note that due to our
specific choice of hybrid profile scales in eq. (2.32), differences between
the two prescriptions only start to appear when turning off the resummation,
such that $g_{\mathrm{run}}$ is nonzero. E.g. for a high
${\mathcal{T}_{\mathrm{cut}}}=100\,\mathrm{GeV}$, which is also a good proxy
for the inclusive $q_{T}$ spectrum, the two prescriptions fully agree in the
canonical region ${q_{T}^{\mathrm{cut}}}\leq 20\,\mathrm{GeV}$ (see the right
panel of figure -1009). This is responsible for the good overall agreement
because most of the cross section is concentrated in the canonical region.

Figure -1008: The ${\mathcal{T}}$ spectrum with a cut on $q_{T}$ as a function
of ${q_{T}^{\mathrm{cut}}}$ at fixed ${\mathcal{T}}=5\,\mathrm{GeV}$ (left)
and as a function of ${\mathcal{T}}$ at fixed
${q_{T}^{\mathrm{cut}}}=100\,\mathrm{GeV}$ (right). The bands indicate the
total perturbative uncertainty $\Delta_{\mathrm{total}}$, see section 3.4. The
colors correspond to different scale setting prescriptions (default: solid
orange); see the text for details.

Figure -1007: The $q_{T}$ spectrum with a cut on ${\mathcal{T}}$ as a function
of ${\mathcal{T}_{\mathrm{cut}}}$ for $q_{T}=15\,\mathrm{GeV}$ (left) and as a
function of $q_{T}$ for ${\mathcal{T}_{\mathrm{cut}}}=100\,\mathrm{GeV}$
(right). The bands indicate the total perturbative uncertainty
$\Delta_{\mathrm{total}}$, see section 3.4. The colors correspond to different
scale setting prescriptions (default: solid orange); see the text for detail.

The comparison of ${\mathcal{T}}$ versus ${\mathcal{T}_{\mathrm{cut}}}$ scales
is much less favorable, with the former failing to reproduce the latter’s
inclusive limit within uncertainties in all cases. This is in line with the
discrepancy reported in ref. [92] for a single-differential measurement of
thrust in $e^{+}e^{-}$ collisions and at a comparable working order
(NLL${}^{\prime}+$NLO). The mismatch is most striking between the default
scales (solid orange) and $(q_{T},{\mathcal{T}})$ scales (dashed blue) in
figures -1009 and -1007, implying that more effort is required to ensure both
a correct integral and the best possible prediction for the shape of the
double-differential spectrum.

From our previous discussion we conclude that the mismatch mostly reduces to
the question of differential versus cumulant scale setting in ${\mathcal{T}}$
alone, so that the methods developed for the single-differential case in refs.
[92, 50] can be brought to bear here as well if desired. However, since this
is a well-known issue that is merely inherited from the single-differential
case, we do not pursue this further here.

Instead, we consider a modification of our profile scales to illustrate that
the issue is indeed a correlated higher-order effect related to scale choices.
Specifically, we can consider lowering the canonical low scale
$\mu_{S}^{\mathrm{I}}\sim(\mu_{B}^{\mathrm{I}})^{2}/\mu_{H}^{\mathrm{I}}\sim{\mathcal{T}}$
in SCETI by a factor of $c=0.5$, which does not parametrically violate the
canonical scaling. Including a smooth interpolation to the fixed-order and
nonperturbative region, this can be achieved by replacing eq. (2.14) with

$\displaystyle f_{\mathrm{run}}^{\mathrm{I}}(c;x)$
$\displaystyle=\begin{cases}x_{0}\Bigl{(}1+\frac{c^{2}x^{2}}{4x_{0}^{2}}\Bigr{)}&x\leq
2x_{0}/c\,,\\\ cx&2x_{0}/c<x\leq x_{1}\,,\\\
cx+\frac{(2-cx_{2}-cx_{3})(x-x_{1})^{2}}{2(x_{2}-x_{1})(x_{3}-x_{1})}&x_{1}<x\leq
x_{2}\,,\\\
1-\frac{(2-cx_{1}-cx_{2})(x-x_{3})^{2}}{2(x_{3}-x_{1})(x_{3}-x_{2})}&x_{2}<x\leq
x_{3}\,,\\\ 1&x_{3}<x\,,\end{cases}$ (3.32)

and keeping the entire remaining profile setup unchanged; setting $c=1$
recovers eq. (2.14).

Figure -1006: Left: The double cumulant cross section as a function of
${\mathcal{T}_{\mathrm{cut}}}$ for ${q_{T}^{\mathrm{cut}}}=100\,\mathrm{GeV}$
for different scale setting prescriptions, with a modified slope $c=0.5$ of
the SCETI profile scales, see eq. (3.32). Right: The $q_{T}$ spectrum with a
cut on ${\mathcal{T}}$ as a function of ${\mathcal{T}_{\mathrm{cut}}}$ for
different scale setting prescriptions, also using modified SCETI profile
scales with $c=0.5$. The bands indicate the total perturbative uncertainty
$\Delta_{\mathrm{total}}$, see section 3.4.

Our results using eq. (3.32) are shown in figure -1006, where we repeat the
left panels of figures -1009 and -1007 using the modified setup. Note that for
simplicity, we use eq. (3.32) for all results in this figure, i.e., for both
differential and cumulant scale setting. We find that the simple modification
eq. (3.32) already substantially improves the agreement between differential
and cumulant scale setting, with the result from
$({q_{T}^{\mathrm{cut}}},{\mathcal{T}})$ scales (dotted green, left panel)
covering the inclusive fixed-order cross section and the result from
$(q_{T},{\mathcal{T}})$ scales (dashed blue, right panel) covering the result
from single-differential $q_{T}$ resummation, at the price of much larger
uncertainties.

We conclude that with additional effort, e.g. applying the methods used in
refs. [50, 92], it would be possible to fully reconcile the best possible
predictions for both the differential shape and the cumulant of the double-
differential spectrum. However, for our purposes we can simply use the
appropriate scale setting for the observable of interest. In particular, if
the experimental observable of interest has cumulant-like character in either
$q_{T}$ or ${\mathcal{T}}$, e.g. if large bins in either observable are
considered, the double-differential profile setup given in this paper, using
$({q_{T}^{\mathrm{cut}}},{\mathcal{T}})$ or
$(q_{T},{\mathcal{T}_{\mathrm{cut}}})$ scales as appropriate, will be
completely sufficient.

##  4 Results

In this section we present our results for Drell-Yan production $pp\to
Z/\gamma^{\ast}\to\ell^{+}\ell^{-}$ at the LHC, with a simultaneous
measurement of the transverse momentum $q_{T}$ of the lepton pair and the
0-jettiness event shape ${\mathcal{T}}$. The center-of-mass energy is taken to
be $E_{\mathrm{cm}}=13\,\mathrm{TeV}$. We assume that in addition, the
invariant mass $Q$ of the lepton pair is measured, and write $pp\to Z$ for
short if $Q=m_{Z}$, and $pp\to Z^{\ast}$ otherwise. The subsequent decay and
the contribution from the virtual photon are included in either case.

To obtain numerical results for the SCETI, SCETII, and SCET+ contributions, we
have implemented all pieces of the relevant double-differential factorized
cross sections to $\mathcal{O}(\alpha_{s})$ and their RGEs to NNLL in SCETlib
[54]. The fixed NLO contributions in full QCD are obtained from MCFM 8.0 [55,
56, 57]. We make use of the MMHT2014nnlo68cl [58] NNLO PDFs with five-flavor
running and $\alpha_{s}(m_{Z})=0.118$. Since we focus on the perturbative
calculation and do not include any nonperturbative effects, we provide the
results down to $1\,\mathrm{GeV}$ in $q_{T}$ and ${\mathcal{T}}$.

The outline of this section is as follows: In section 4.1 we present our fully
resummed prediction for the double-differential spectrum, both as surface
plots over the $(q_{T},{\mathcal{T}})$ plane and for selected slices along
lines of constant $q_{T}$ or ${\mathcal{T}}$. We demonstrate that our
prediction smoothly interpolates between the SCETI and SCETII boundary
theories, i.e., we show that our matching formula in eq. (3.2) recovers the
matched predictions on either boundary and improves over them by an additional
resummation of power-suppressed terms. Finally, in section 4.2 we present our
predictions for the single-differential spectra
$\mathrm{d}\sigma({q_{T}^{\mathrm{cut}}})/\mathrm{d}{\mathcal{T}}$ and
$\mathrm{d}\sigma({\mathcal{T}_{\mathrm{cut}}})/\mathrm{d}q_{T}$ with a cut on
the other variable, and show how they recover the inclusive single-
differential ${\mathcal{T}}$ and $q_{T}$ spectra for large values of
${q_{T}^{\mathrm{cut}}}$ and ${\mathcal{T}_{\mathrm{cut}}}$, respectively.

###  4.1 Double spectrum and comparison with boundary theories

Figure -1005: The double-differential Drell-Yan cross section at fixed NLO
(top), resummed NNLL${}_{\mathcal{T}}+$NLO (center), and NNLL${}_{q_{T}}+$NLO
(bottom). The resummed predictions are obtained by using only SCETI (SCETII)
renormalization group evolution to resum logarithms of ${\mathcal{T}}$
($q_{T}$), as outlined in section 2.2 (section 2.3), and matching the result
to the fixed-order cross section. For better visibility we show the spectrum
with respect to $\log_{10}q_{T}$ and $\log_{10}{\mathcal{T}}$. On the rear
walls we show the result of integrating the double spectrum over either
variable up to a cut at $100\,\mathrm{GeV}$.

To highlight the necessity of the simultaneous resummation of large logarithms
of both $q_{T}$ and ${\mathcal{T}}$, we start by showing results for the
double spectrum (the cross section double-differential in $q_{T}$ and
${\mathcal{T}}$) where only some of the logarithms are resummed. These results
are shown as surface plots in figure -1005, where we plot the double-
differential spectrum with respect to $\log_{10}q_{T}$ and
$\log_{10}{\mathcal{T}}$ for better visibility. In each case the left rear
wall of the surface plot shows the result of integrating the double-
differential spectrum up to ${\mathcal{T}_{\mathrm{cut}}}=100\,\mathrm{GeV}$,
but staying differential in $\log_{10}q_{T}$. Similarly, the right rear wall
shows the projection onto the single-differential spectrum in
$\log_{10}{\mathcal{T}}$, with a cut at
${q_{T}^{\mathrm{cut}}}=100\,\mathrm{GeV}$.555We refer the reader to section
3.5 for the precise way we perform these integrals.

The top left panel of figure -1005 shows the spectrum evaluated at fixed
$\mathcal{O}(\alpha_{s})$, without any resummation. The double-differential
fixed-order spectrum diverges logarithmically for small ${\mathcal{T}}$ at any
value of $q_{T}$, while its projections onto the single-differential spectra
in $q_{T}$ and ${\mathcal{T}}$ feature double-logarithmic singularities.
Notably, the double-differential spectrum has a sharp kinematic edge along
$q_{T}={\mathcal{T}}$. This sharp edge is unphysical because it only reflects
the kinematics of a single on-shell emission with transverse momentum $k_{T}$
at rapidity $\eta$, which contributes at most
${\mathcal{T}}=k_{T}\,e^{-\lvert\eta\rvert}\leq k_{T}=q_{T}$. Due to the
vectorial nature of $q_{T}$, however, back-to-back emissions can populate the
region ${\mathcal{T}}>q_{T}$ at higher orders, and the kinematic edge must be
smeared out.

Next, we consider the cases in which only logarithms of one variable are
resummed, while logarithms involving the auxiliary variable are treated at
fixed order. In the middle panel of figure -1005, we show the result of
resumming logarithms of ${\mathcal{T}}$ using the SCETI matched result in eq.
(2.11). The resummation is performed at NNLL and is matched to full NLO, which
we refer to as NNLL${}_{\mathcal{T}}+$NLO. As discussed in section 2.2, this
prediction is valid as long as the parametric relation ${\mathcal{T}}\ll
q_{T}\sim\sqrt{Q{\mathcal{T}}}$ is satisfied. This corresponds to the SCETI
phase-space boundary (blue) in figure -1018, running from the region of small
${\mathcal{T}}$ and intermediate $q_{T}$ towards the fixed-order region where
$q_{T}\sim{\mathcal{T}}\sim Q$. It is clear that away from its region of
validity, the NNLL${}_{\mathcal{T}}+$NLO result contains unresummed logarithms
of $q_{T}$ because at any point in ${\mathcal{T}}$, the prediction diverges
for very small $q_{T}$. In particular, power corrections of
$\mathcal{O}({\mathcal{T}}^{2}/q_{T}^{2})$ are only captured by the fixed-
order matching. They become $\mathcal{O}(1)$ as one approaches the diagonal
${\mathcal{T}}=q_{T}$ (the green line in figure -1018), and encode the phase-
space boundary at $q_{T}\sim{\mathcal{T}}$. As in the NLO case, treating this
phase-space boundary at fixed order leads to the sharp kinematic edge along
the diagonal; physically, the all-order tower of collinear emissions that
contribute to $q_{T}$ in SCETI cannot resolve the boundary because it arises
from the dynamics at central rapidities. The projections onto the rear walls
highlight that only ${\mathcal{T}}$ is resummed. The single-differential
$q_{T}$ spectrum still diverges as $q_{T}\to 0$, while the ${\mathcal{T}}$
spectrum features a physical Sudakov peak.

Figure -1004: The double-differential Drell-Yan cross section at NNLL$+$NLO,
at $Q=m_{Z}$ (top) and $Q=300\,\mathrm{GeV}$ (bottom), with respect to
$\log_{10}q_{T}$ and $\log_{10}{\mathcal{T}}$. On the rear walls we show the
result of integrating the double spectrum over either variable up to a cut at
$100\,\mathrm{GeV}$. The contour plots indicate total perturbative
uncertainties relative to the cross section,
$\Delta_{\mathrm{total}}=\Delta_{+}\oplus\Delta_{\mathrm{I}}\oplus\Delta_{\mathrm{II}}\oplus\Delta_{\mathrm{FO}}$.
The contour plots are left blank in the region where
$\mathrm{d}\sigma/(\mathrm{d}Q\,\mathrm{d}\log_{10}q_{T}\,\mathrm{d}\log_{10}{\mathcal{T}})$
is less than $3\%$ of its peak height.

In the bottom panel of figure -1005, we show the result of resumming
logarithms of (the $b_{T}$ variable conjugate to) $q_{T}$ to NNLL and matching
to fixed NLO, using the SCETII matched result in eq. (2.29). We denote this
order by NNLL${}_{q_{T}}$+NLO. This result is valid for ${\mathcal{T}}\sim
q_{T}\ll\sqrt{Q{\mathcal{T}}}$, i.e., around the SCETII phase-space boundary
(green) in figure -1018, where we find the onset of a Sudakov peak from the
$q_{T}$ resummation and a smooth kinematic suppression towards
${\mathcal{T}}\gg q_{T}$. However, the NNLL${}_{q_{T}}$+NLO result diverges
for smaller values of ${\mathcal{T}}$. This is due to unresummed logarithms of
${\mathcal{T}}$ in both the factorized cross section in SCETII and terms of
$\mathcal{O}(q_{T}^{2}/(Q{\mathcal{T}}))$ that are treated at fixed order as
part of the matching correction. In this case the single-differential
projections show a Sudakov peak in $q_{T}$, but a logarithmic divergence at
small ${\mathcal{T}}$.

Figure -1003: Breakdown of resummation uncertainties contributing to the
relative uncertainty in the top right panel of figure -1004, showing (from
left to right) $\Delta_{\mathrm{I}}$, $\Delta_{\mathrm{II}}$, and
$\Delta_{+}$. As in figure -1004 we leave regions blank where the cross
section is small.

Our final results for the Drell-Yan double spectrum are shown in figure -1004,
as given by the fully matched prediction in eq. (3.2). Here all resummed
contributions are evaluated at NNLL, and we match to fixed NLO. This achieves,
for the first time, the complete resummation of all large logarithms in the
double spectrum, so we simply refer to this order as NNLL$+$NLO. The top row
of plots is for $Q=m_{Z}$, i.e., for Drell-Yan production at the $Z$ pole. In
the bottom row we consider $Q=300\,\mathrm{GeV}$ as a representative phase-
space point at higher production energies. Our matched and fully resummed
double spectrum features a two-dimensional Sudakov peak that is situated
between the two parametric phase-space boundaries (cf. figure -1018), is
smoothly suppressed beyond, and shifts towards higher values of $q_{T}$ and
${\mathcal{T}}$ for $Q=300\,\mathrm{GeV}$, as expected. Integrating the double
spectrum over either variable also results in a physical Sudakov peak, as can
be seen from the projections onto the rear walls. Up to small differences in
scale setting discussed in section 3.5, the left and right rear walls agree
with the result of integrating the NNLL${}_{q_{T}}+$NLO and
NNLL${}_{\mathcal{T}}+$NLO results in figure -1005 over ${\mathcal{T}}$ and
$q_{T}$, respectively. The contour plots in figure -1004 show the total
perturbative uncertainties $\Delta_{\mathrm{total}}$ as percent deviations
from the central result for the double spectrum. As described in section 3.4,
$\Delta_{\mathrm{total}}$ combines estimates of all sources of resummation
uncertainty in the prediction.

In figure -1003, we break down the uncertainty for the Drell-Yan double-
differential spectrum at $Q=m_{Z}$ into its contributions from SCETI, SCETII