Joint Two-Dimensional Resummation in $q_T$ and $0$-Jettiness at NNLL

TL;DR

This paper develops a novel two-dimensional resummation technique for initial-state radiation in Drell-Yan processes, simultaneously resumming logarithms in both $q_T$ and $0$-jettiness at NNLL order, providing more precise theoretical predictions.

Contribution

It introduces the first analytic two-dimensional Sudakov resummation for initial-state radiation, combining $q_T$ and $0$-jettiness resummation at NNLL order with a new impact-parameter space method.

Findings

First two-dimensional analytic Sudakov resummation for initial-state radiation.

Resummation matches fixed order results when integrated over one variable.

Method applicable to any color-singlet production process.

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.

Tables

Table 1. Table 1 : Summary of canonical scales in SCET I , SCET + , and SCET II [see eq. ( 2.10 ), ( 2.3 ) and ( 2.4 )]. For SCET + and SCET II we give the canonical scales in b T subscript 𝑏 𝑇 b_{T} space.

Scale	SCETI	SCET+	SCETII
${\mu }_H$	$Q$	$Q$	$Q$
${\mu }_B$	$\sqrt{𝒯Q}$	$b_0/b_T$	$b_0/b_T$
${\nu }_B$		$Q$	$Q$
${\mu }_{𝒮}$		$b_0/b_T$
${\nu }_{𝒮}$		${\left(b_0/b_T\right)}^2/𝒯$
${\mu }_S$	$𝒯$	$𝒯$	$b_0/b_T$
${\nu }_S$			$b_0/b_T$

Equations

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

$$q_a^{\mu }=\frac{Q{e}^{+Y}}{2}{n}_a^{\mu },q_b^{\mu }=\frac{Q{e}^{-Y}}{2}{n}_b^{\mu },$$

$$Q_a=Q_b=Q.$$

$$n_a\text{-collinear:}$$

$$n_b\text{-collinear:}$$

$$\displaystyle\quad p^{\mu}\sim\Bigl{(}{\mathcal{T}},{\mathcal{T}},{\mathcal{T}}\Bigr{)}\,,$$

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

$$\frac{\mathrm{d}{\sigma }_{\mathrm{I}}}{\mathrm{d}Q\mathrm{d}Y\mathrm{d}{q}_T\mathrm{d}\mathcal{T}}$$

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

$$\frac{\mathrm{d}\sigma }{\mathrm{d}Q\mathrm{d}Y\mathrm{d}{q}_T\mathrm{d}\mathcal{T}}$$

$$x_a=\frac{Q}{E_{\mathrm{cm}}}{e}^{+Y},x_b=\frac{Q}{E_{\mathrm{cm}}}{e}^{-Y}.$$

$$B_q(t,x,{\vec{k}}_T,\mu )$$

$${\mu }_H^{\mathrm{I}}\sim Q,{\mu }_B^{\mathrm{I}}\sim \sqrt{Q\mathcal{T}},{\mu }_S^{\mathrm{I}}\sim \mathcal{T}.$$

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

$${\mu }_B^{\mathrm{I}}(\mathcal{T}),{\mu }_S^{\mathrm{I}}(\mathcal{T})\to {\mu }_H^{\mathrm{I}}={\mu }_{\mathrm{FO}}\text{for}\mathcal{T}\to Q.$$

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

$$f_{\mathrm{run}}^{\mathrm{I}}(x)$$

$${\mu }_S^{\mathrm{I}}$$

$${\mu }_B^{\mathrm{I}}$$

$$f_{\mathrm{vary}}(x)=\{\begin{array}{ll}2(1-x^2/x_3^2)& 0\le x<x_3/2,\\ 1+2(1-x/x_3{)}^2& x_3/2\le x<x_3,\\ 1& x_3\le x.\end{array}$$

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

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

$$n_a\text{-collinear:}$$

$$n_b\text{-collinear:}$$

$$\displaystyle\quad p^{\mu}\sim\Bigl{(}q_{T},q_{T},q_{T}\Bigr{)}\sim\Bigl{(}{\mathcal{T}},{\mathcal{T}},{\mathcal{T}}\Bigr{)}\,.$$

$$\frac{\mathrm{d}{\sigma }_{\mathrm{II}}}{\mathrm{d}Q\mathrm{d}Y\mathrm{d}{q}_T\mathrm{d}\mathcal{T}}$$

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

$$\frac{\mathrm{d}\sigma }{\mathrm{d}Q\mathrm{d}Y\mathrm{d}{q}_T\mathrm{d}\mathcal{T}}$$

$${\omega }_a=Q{e}^{+Y}=x_a{E}_{\mathrm{cm}},{\omega }_b=Q{e}^{-Y}=x_b{E}_{\mathrm{cm}},$$

$$B_q(\omega ,{\vec{k}}_T,\mu ,\nu )$$

$${\mu }_H^{\mathrm{II}}$$

$${\nu }_B^{\mathrm{II}}$$

$${\mu }_H^{\mathrm{II}}$$

$${\nu }_B^{\mathrm{II}}$$

$${\tilde{\gamma }}_{\nu }^i(b_T,\mu )$$

Full text

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.

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

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

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,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

Based on figure -1137, 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 -1137.

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]

[TABLE]

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

[TABLE]

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

[TABLE]

and taking the maximum absolute deviation from the central result. These variations are also indicated in the right panel of figure -1137. 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,

[TABLE]

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

[TABLE]

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

[TABLE]

The factorization receives power corrections of the form

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

[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

[TABLE]

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.

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,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

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,

[TABLE]

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 -1136. In practice we pick

[TABLE]

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

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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]

[TABLE]

which holds up to power corrections

[TABLE]

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,

[TABLE]

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,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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 -1138. 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 -1138, 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,

[TABLE]

This is illustrated in figure -1134, 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.

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,

[TABLE]

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

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,

[TABLE]

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

As a final important consequence of figure -1134, 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,

[TABLE]

In figure -1131 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 -1137 and figure -1136. 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 -1131, 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:

[TABLE]

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,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

It follows that

[TABLE]

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

[TABLE]

and consequently

[TABLE]

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,

[TABLE]

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.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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

[TABLE]

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,

[TABLE]

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

[TABLE]

where the polynomials governing the interpolation between zero and one are

[TABLE]

The transition points $a_{1,\ldots,6}$ determine the transition between the different regions, as can be seen from the helper functions in figure -1130: 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 -1132. (In both figures -1133 and -1132, corresponding values of $a$ are indicated on the horizontal axis at the top of the panels.) On the other hand, from figure -1133 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 -1133, 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})$,

[TABLE]

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,

[TABLE]

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,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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,

[TABLE]

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

[TABLE]

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.

In figures -1129 to -1127, 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 -1129, 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 -1128 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 -1127. In figures -1128 and -1127 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 -1127 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 -1129 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 -1128 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 -1129). This is responsible for the good overall agreement because most of the cross section is concentrated in the canonical region.

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 -1129 and -1127, 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

[TABLE]

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

Our results using eq. (3.32) are shown in figure -1126, where we repeat the left panels of figures -1129 and -1127 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

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 -1125, 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 -1125 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 -1125, 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 -1138, 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 -1138), 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.

In the bottom panel of figure -1125, 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 -1138, 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}}$.

Our final results for the Drell-Yan double spectrum are shown in figure -1124, 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 -1138), 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 -1125 over ${\mathcal{T}}$ and $q_{T}$, respectively. The contour plots in figure -1124 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 -1123, we break down the uncertainty for the Drell-Yan double-differential spectrum at $Q=m_{Z}$ into its contributions from SCETI, SCETII and SCET+ resummation uncertainties, respectively. As expected, the SCETI resummation uncertainty dominates in the SCETI region of phase space, and similarly for SCETII. The SCET+ resummation uncertainty is largest along the phase-space boundaries, indicating that it is mostly sensitive to variations of the transition points, i.e., the points where the intrinsic SCET+ resummation is turned off in our matched prediction.

To further highlight the differences between our fully double-differential resummation and the single-differential resummation at either NNLL${}_{q_{T}}$ or NNLLT, we take slices of the surface plots and overlay them in figure -1122, keeping $q_{T}$ (left) or ${\mathcal{T}}$ (right) fixed. The solid red curve corresponds to the matched and fully resummed cross section in eq. (3.2), with the uncertainty band given by the total perturbative uncertainty $\Delta_{\mathrm{total}}$, see eq. (3.26). The matched SCETI (dashed blue) and SCETII (dotted green) predictions correspond to the middle and bottom panel of figure -1125, respectively. Their uncertainty bands are given by $\Delta^{\mathrm{I}}_{\mathrm{total}}$ and $\Delta^{\mathrm{II}}_{\mathrm{total}}$, which only probe a subset of higher-order terms as predicted by the respective RGE, see eqs. (2.18) and (2.37). The SCETI prediction features an unphysical sharp edge at ${\mathcal{T}}=q_{T}$, cf. the middle panel of figure -1125, and for this reason is cut off at ${\mathcal{T}}=0.8\,q_{T}$.

All panels in figure -1122 show that our final prediction smoothly interpolates between the SCETI and SCETII boundary theories, both for the central values and for the uncertainties. Specifically, the matched prediction tends towards SCETI (SCETII) for small (large) values of ${\mathcal{T}}$ and large (small) values of $q_{T}$. In the left column one clearly sees that SCETII only captures logarithms of ${\mathcal{T}}$ at fixed order, leading to a diverging spectrum as ${\mathcal{T}}\to 0$, while the complete NNLL result features a physical Sudakov peak. Conversely, the SCETI result diverges as $q_{T}\to 0$, but is rendered physical by the additional $q_{T}$ resummation at NNLL.

We would like to stress that our fully resummed prediction does not exactly agree with either boundary theory, even beyond the final transition points $a_{1}$ and $a_{6}$ where the intrinsic SCET+ resummation is turned off. The reason for this is that even in these limits, the matched cross section in eq. (3.2) improves over the matched SCETI and SCETII cross sections in eq. (2.11) and eq. (2.29) by an additional resummation of power-suppressed terms, cf. eqs. (3.2) and (3.2). To assess the size of this effect, we again compare both single-differential resummations (dashed blue and dotted green) to our matched prediction (solid red) in figure -1121, but for reference include the case where $\sigma_{+}$ in the matched prediction is evaluated at $\mu^{\mathrm{I}}$ (solid blue) or $\mu^{\mathrm{II}}$ (solid green) directly. One can easily verify from e.g. the right panel that for $q_{T}$ above the right-most vertical line (where $a<a_{1}$), the difference between the solid blue and the dashed blue curves indeed amounts to a small power-suppressed set of higher-order terms, while our best prediction (solid red) recovers the solid blue curve as it must. Similarly, for $q_{T}$ below the left-most vertical line (where $a>a_{6}$), the difference between the solid green (and solid red) and dashed green curves can be seen to be a small correction, reflecting the size of power-suppressed higher-order terms predicted by the SCETI RGE in this region. The asymptotic limits are interchanged in the left panel, where $a<a_{1}$ towards the left and $a>a_{6}$ towards the very right of the plot.

4.2 Single-differential spectra with a cut on the other variable

So far we have turned our attention to the cross section differential in both $q_{T}$ and ${\mathcal{T}}$. In addition to this double spectrum, our setup also predicts the fully matched and resummed double cumulant cross section, and the single-differential $q_{T}$ (or ${\mathcal{T}}$) spectrum with a cut on the other variable; selected results for these observables where already discussed in section 3.5 from a more technical point of view. In figure -1120, we show some more detailed results for the single-differential spectra with an additional cut, where the left panel shows $\mathrm{d}\sigma({q_{T}^{\mathrm{cut}}})/\mathrm{d}{\mathcal{T}}$ as a function of ${\mathcal{T}}$ for various values of ${q_{T}^{\mathrm{cut}}}$, and the right panel shows $\mathrm{d}\sigma({\mathcal{T}_{\mathrm{cut}}})/\mathrm{d}q_{T}$ as a function of $q_{T}$ for various values of ${\mathcal{T}_{\mathrm{cut}}}$. By increasing the value of the cut, they can be seen to approach the inclusive single-differential spectra (orange solid), with which they must agree when sending ${q_{T}^{\mathrm{cut}}}\to\infty$ or ${\mathcal{T}_{\mathrm{cut}}}\to\infty$, respectively. This is by construction because we employ cumulant scale setting as appropriate for this prediction, cf. section 3.5. We observe that cuts on the other variable shape either spectrum in a very nontrivial way. Tight cuts $\lesssim 10\,\mathrm{GeV}$ push the peak to lower values and suppress the tail, where the $q_{T}$ spectrum is somewhat more resilient to cuts on ${\mathcal{T}}$ than vice versa. Intermediate cuts $\sim 10-15\,\mathrm{GeV}$ keep the peak and mostly lead to a suppression in the tail, while the effect of cuts $\gtrsim 40\,\mathrm{GeV}$ is almost negligible in the $q_{T}$ and ${\mathcal{T}}$ ranges of interest.

5 Conclusions

In this paper we calculated the Drell-Yan cross section double-differential in the transverse momentum $q_{T}$ of the lepton pair and the 0-jettiness ${\mathcal{T}}$. Both ${\mathcal{T}}$ and $q_{T}$ probe the initial state radiation, leading to Sudakov double logarithms of ${\mathcal{T}}/Q$ and $q_{T}/Q$ in the cross section. We performed, for the first time, the simultaneous resummation of both kinds of logarithms, achieving next-to-next-to-leading logarithmic accuracy and matching the result to next-to-leading fixed order. We accomplish this resummation by using SCETI and SCETII to describe the regions ${\mathcal{T}}\ll q_{T}\sim\sqrt{{\mathcal{T}}Q}$ and ${\mathcal{T}}\sim q_{T}\ll\sqrt{{\mathcal{T}}Q}$, respectively, and SCET+ to describe the bulk of phase space in between these boundaries [51].

Obtaining reliable numerical predictions required several nontrivial steps: (1) Matching several factorized cross sections for the different regions of phase space, for which we use a Venn-diagram method to avoid double counting. (2) Choosing appropriate profile scales for the various ingredients in the factorization formulas that respect all relevant canonical scaling relations and at the same time smoothly interpolate between the different regions of phase space, and varying these scales to estimate perturbative uncertainties. This is significantly more involved than in the usual single-differential case, and is further complicated by the requirement to choose scales in impact parameter ($b_{T}$) space for SCETII. For example, the rapidity scale for the collinear-soft function in SCET+ has a canonical scaling that does not coincide with any scale on the SCETI and SCETII boundaries. (3) Ensuring that scales and scale variations are still, to the extent possible, inherited from the single-differential resummation of ${\mathcal{T}}$ and $q_{T}$. This makes our setup flexible enough to incorporate other procedures for estimating the uncertainty in the individual resummations. (4) To handle the transition between SCETI, SCET+ and SCETII, we introduced profile scales in terms of a regime parameter $a$, designed such that $a=1$ for SCETI and $a=2$ for SCETII. The precise transition points in $a$ were chosen by comparing the various singular and nonsingular cross section, and are varied as part of the uncertainty estimate. (5) We also introduced a new hybrid (i.e., $q_{T}$ and $b_{T}$ dependent) scale choice for $q_{T}$ resummation, allowing the resummation to strictly take place in $b_{T}$ space, while turning the resummation on and off using $q_{T}$.

We demonstrated that our simultaneous resummation of ${\mathcal{T}}$ and $q_{T}$ yields the correct resummed single-differential cross sections after integrating over either ${\mathcal{T}}$ or $q_{T}$. This requires choosing scales at the level of the differential or integrated (cumulative) cross section as appropriate, which we discussed in detail.

While the predictions obtained here are of some direct phenomenological interest, as ${\mathcal{T}}$ has been measured in bins of $q_{T}$ [48], our analysis is also an important step towards precise and differential predictions for LHC cross sections in general. Specifically, the Monte Carlo event generator Geneva [49, 50] is based on a NNLL*′* resummed prediction for the cross section differential in ${\mathcal{T}}$, and would benefit from the simultaneous resummation of $q_{T}$. Indeed, our NNLL results clearly indicate that only resumming the logarithms of either ${\mathcal{T}}$ or $q_{T}$ gives a poor description of the double-differential cross section. Our methods apply at any order and for any color-singlet production process, allowing for a straightforward extension once the relevant perturbative ingredients become available. We hope that our analysis can pave the way for going beyond single-differential resummations in many other contexts as well.

Acknowledgments

We thank Goutam Das, Markus Diehl, and Markus Ebert for discussions. We would also like to thank Markus Ebert for his contributions to SCETlib. This work is supported by the ERC grant ERC-STG-2015-677323 and the D-ITP consortium, a program of the Netherlands Organization for Scientific Research (NWO) that is funded by the Dutch Ministry of Education, Culture and Science (OCW). FT and JM thank Nikhef for hospitality. GL thanks DESY for hospitality and the DESY theory group thanks GL for an ample supply of Dutch pepernoten.

Appendix A Plus distributions and Fourier transform

We use the following standard plus distributions with dimensionless arguments,

[TABLE]

In intermediate steps we need a two-dimensional plus distribution originally defined in ref. [51],

[TABLE]

Our shorthands for distributions with dimensionful arguments in one spatial dimension are

[TABLE]

In terms of the second class of (power-like) plus distributions, we further define

[TABLE]

which have a group property, assuming identical boundary condition $\mu$,

[TABLE]

Shifting the boundary condition of $\mathcal{V}_{a}(k,\mu)$ from $\mu$ to $\mu^{\prime}$ is also straightforward,

[TABLE]

We use the conventions from app. C of ref. [77] for logarithmic plus distributions in two integer spatial dimensions, with $k_{T}^{2}\equiv\vec{\kern 0.0ptk}_{T}^{2}\geq 0$,

[TABLE]

Our convention for the two-dimensional Fourier transform also follows ref. [77],

[TABLE]

Here we make the mass dimension of $\mathrm{d}f/\mathrm{d}\vec{\kern 0.0ptp}_{T}=\mathrm{d}f/(\mathrm{d}p_{x}\,\mathrm{d}p_{y})$ explicit. If $f$ is azimuthally symmetric, i.e., if for $p_{T}\equiv\lvert\vec{\kern 0.0ptp}_{T}\rvert,b_{T}\equiv\lvert\vec{\kern 0.0ptb}_{T}\rvert$,

[TABLE]

the azimuthal integral can be performed, leaving

[TABLE]

where $J_{0}(x)$ is the zeroth-order Bessel function of the first kind. Integrating the first expression in eq. (A.11) by parts, the cumulant in $p_{T}$ is given by

[TABLE]

where $J_{1}(x)$ is the first-order Bessel function of the first kind. Fourier transforms of $\mathcal{L}_{n}(\vec{\kern 0.0ptp}_{T},\mu)$ can be found in table 5 of ref. [77], and are most conveniently expressed in terms of

[TABLE]

Appendix B Perturbative ingredients

B.1 Anomalous dimensions

We expand the $\beta$ function of QCD as

[TABLE]

The coefficients in the $\overline{\text{MS}}$ scheme are, up to three loops [93, 94],

[TABLE]

We work with $n_{f}=5$ light flavors. The cusp and noncusp anomalous dimensions are expanded as

[TABLE]

The coefficients of the $\overline{\mathrm{MS}}$ cusp anomalous dimension to three loops are [95, 96, 97]

[TABLE]

where here and in the following $C_{i}=C_{F}\,(C_{A})$ for $i=q\,(g)$. The fixed-order boundary condition of the resummed rapidity anomalous dimension eq. (2.26) reads, through NNLL,

[TABLE]

where we have already used that $\gamma^{i}_{\nu\,0}=0$. For our choice of regulator, the two-loop boundary condition is given by [70]

[TABLE]

Hard function.

The hard function $H_{\kappa}(Q,\mu)$ is proportional to the square of the hard matching coefficient [see eq. (B.2)]. The hard matching coefficient for $q\bar{q}\to Z/\gamma^{\ast}$ is renormalized as

[TABLE]

The coefficients of the quark noncusp anomalous dimension up to two loops are

[TABLE]

Beam functions and PDFs.

In $b_{T}$ space, the TMD beam function is renormalized as

[TABLE]

We include a tilde to indicate that $\tilde{\gamma}_{B}^{q}$ is related to the SCET${}_{\text{II}}$ beam function, even though it does not depend on $b_{T}$. Its coefficients through two loops are [70]

[TABLE]

The resummed rapidity anomalous dimension $\tilde{\gamma}^{i}_{\nu}(b_{T},\mu)$ is given in eq. (2.26).

The double-differential beam function satisfies the same RGE as the inclusive SCET${}_{\text{I}}$ beam function,666We note that ref. [51] incorrectly did not distinguish between $\gamma^{i}_{B}(\alpha_{s})$ and $\tilde{\gamma}^{i}_{B}(\alpha_{s})$. This lead to the noncusp contribution to the collinear-soft anomalous dimension being missing in their eq. (3.26), cf. our corrected eq. (B.1) and the nonvanishing two-loop noncusp coefficient in our eq. (B.1).

[TABLE]

The coefficients of the SCET${}_{\text{I}}$ quark beam anomalous dimension are [60, 98]

[TABLE]

We also require the one-loop coefficients of the PDF anomalous dimension,

[TABLE]

Note that we expand $P_{ij}(\alpha_{s},z)$ in $\alpha_{s}/(4\pi)$. The one-loop coefficients are

[TABLE]

in terms of the standard color-stripped one-loop QCD splitting functions

[TABLE]

Soft and collinear-soft functions.

The RGE of the beam thrust soft function reads

[TABLE]

For the double-differential soft function in $b_{T}$ space we have

[TABLE]

where we again use a tilde on the $\mu$ anomalous dimension to indicate that it relates to the SCET${}_{\text{II}}$ soft function. The RGE of the collinear-soft function in $b_{T}$ space reads

[TABLE]

The soft and collinear-soft noncusp anomalous dimension coefficients are only nonzero starting at two loops and can be inferred from consistency,

[TABLE]

B.2 Fixed-order ingredients

Hard process.

The Born cross section for $q{\bar{q}}\to Z/\gamma^{\ast}\to\ell^{+}\ell^{-}$ is given by

[TABLE]

where $Q_{q}$ is the quark charge in units of $\lvert e\rvert$, $v_{\ell,q}$ and $a_{\ell,q}$ are the standard vector and axial couplings of the leptons and quarks, and $m_{Z}$ and $\Gamma_{Z}$ are the mass and width of the $Z$ boson. The one-loop Wilson coefficient $C_{q{\bar{q}}}^{V,A}(Q^{2},\mu)$ from matching the quark current in QCD onto SCET was computed in refs. [99, 100]. This leads to the following hard function [46],

[TABLE]

where $\operatorname{Re}$ denotes the real part.

Beam functions.

The one-loop matching coefficients for the single-differential SCET${}_{\text{II}}$ beam function in $b_{T}$ space are given by [101, 70]

[TABLE]

where $L_{b}$ was defined in eq. (A.13) and the boundary conditions at $L_{b}=0,\nu=\omega$ are

[TABLE]

As for its anomalous dimension, we use a tilde to indicate that these boundary conditions are part of the SCET${}_{\text{II}}$ beam function, even though they do not depend on $b_{T}$.

It is convenient to decompose the matching coefficients for the double-differential SCET${}_{\text{I}}$ quark beam function as

[TABLE]

where $\mathcal{I}_{qj}(t,z,\mu)$ is the matching coefficient for the inclusive quark beam function [46, 60],

[TABLE]

with the finite terms in this case given by

[TABLE]

and using the shorthand

[TABLE]

The $\Delta\mathcal{I}_{qj}$ piece can be interpreted as a correction over the limit $t\ll k_{T}^{2}$, where recoil from collinear radiation is power suppressed and the double-differential beam function becomes proportional to $\delta(\vec{\kern 0.0ptk}_{T})$. Specifically, it scales as

[TABLE]

and by construction satisfies

[TABLE]

At one loop it can be extracted from the full calculation of $\mathcal{I}_{qj}(t,z,\vec{\kern 0.0ptk}_{T},\mu)$ [102, 59] and has the compact form

[TABLE]

The second line is regular in $t$ because the term in square brackets vanishes as $t\to 0$. After accumulating over the transverse plane up to $q_{T}^{\mathrm{cut}}>0$, we have

[TABLE]

Soft and collinear-soft functions.

The (beam) thrust soft function is [103, 104, 46]

[TABLE]

The one-loop collinear-soft function in $b_{T}$ space is [51]

[TABLE]

It is again convenient to decompose the double-differential soft function computed in ref. [51] into separate pieces with distinct power counting,

[TABLE]

Here $S_{i}(\vec{\kern 0.0ptk}_{T},\mu,\nu)$ is the standard single-differential $q_{T}$ soft function, which in $b_{T}$ space at one loop is given by [69]

[TABLE]

The second term in eq. (B.33) can again be interpreted as a correction, in this case over the limit $k\gg\vec{\kern 0.0ptk}_{T}$ where the contribution of soft radiation to the ${\mathcal{T}}=k$ measurement becomes power suppressed. In momentum space this term satisfies

[TABLE]

Equivalently, in position space we have

[TABLE]

At one loop, $\Delta S_{i}$ is given by

[TABLE]

The second line is not yet manifestly independent of $\mu$, but can be simplified noting that

[TABLE]

It is straightforward to show this by writing all three distributions in terms of $\theta(x_{1}-\beta)$ and $\theta(x_{2}-\beta^{2})$ for infinitesimal $\beta$, collecting terms, and noting that the result is finite as $\beta\to 0$. From eq. (B.38) we can immediately read off the fixed-order double cumulant of $\Delta S_{i,1}$ for ${\mathcal{T}}_{\mathrm{cut}}>0$, $q_{T}^{\mathrm{cut}}>0$,

[TABLE]

where the dependence on $\mu$ drops out as expected. Inserting eq. (B.38) and integrating by parts also yields the cumulant up to ${\mathcal{T}}_{\mathrm{cut}}>0$ in position space,

[TABLE]

where $x\equiv b_{T}{\mathcal{T}}_{\mathrm{cut}}$ and $\phantom{}{}_{i}F_{j}(x_{1},\ldots,x_{i};y_{1},\ldots,y_{j};z)$ is the generalized hypergeometric function. The right hand side of eq. (B.40) asymptotes to $1/x^{2}$ as $x\to\infty$, as required by the scaling law in eq. (B.2). We also need the spectrum of $\Delta\tilde{S}_{i,1}$ at ${\mathcal{T}}>0$ in position space,

[TABLE]

where this time $x\equiv b_{T}{\mathcal{T}}$ and the term in square brackets again asymptotes to $1/x^{2}$ as $x\to\infty$.

B.3 Renormalization-group evolution

SCETI.

The closed-form solution of eq. (B.1) is [104, 64]

[TABLE]

where $\mathcal{V}_{\eta}$ was defined in eq. (A.5). Similarly, the solution of eq. (B.1) is

[TABLE]

The definitions of $\eta^{i}_{\Gamma}$, $K^{i}_{\Gamma}$, and $K^{i}_{\gamma}$ and their approximate analytical form at NNLL are given for example in app. A.4 of ref. [105]. The solution of the hard RGE in eq. (B.1) in our notation is also found there.

SCETII.

Evolving first in $\nu$ and then in $\mu$ (from right to left), the solution of eq. (B.1) is

[TABLE]

where the resummed rapidity anomalous dimension $\tilde{\gamma}^{q}_{\nu}(b_{T},\mu_{B})$ is given by eq. (2.26). For the double-differential soft function renormalized as in eq. (B.1), we have

[TABLE]

SCET**+**.

Again evolving first in $\nu$ and then in $\mu$, the solution of the collinear-soft RGE in eq. (B.1) in $b_{T}$ space is given by

[TABLE]

The rapidity evolution factor on the second line does not depend on $k^{\prime}$ and thus may be taken out of the convolution integral.

B.4 Beam function convolutions with RG kernels

In SCETI the combined beam and soft renormalization group running in momentum space has the functional form $\mathcal{V}_{\eta}(k,\mu)$ [see eq. (A.5)], where the convolution variable $k$ is ultimately fixed by the overall ${\mathcal{T}}$ measurement. To evaluate the resummed SCETI cross section as a function of ${\mathcal{T}}$ (${\mathcal{T}}_{\mathrm{cut}}$) we require convolutions of $\mathcal{V}_{\eta}$ with the double-differential beam function at finite ${\mathcal{T}}$ (integrated up to ${\mathcal{T}}_{\mathrm{cut}}$). Convolutions of $\mathcal{V}_{\eta}$ with the inclusive beam function reduce to a linear combination of

[TABLE]

with $Q_{i}$ as in eq. (2.1), and are straightforward to evaluate using results in app. B of ref. [64]. The same is true for convolutions of the evolution kernel with the single-differential ${\mathcal{T}}$ soft function and, in SCET+, with the collinear-soft function. (Cross terms between beam, csoft, or soft functions are absent at one loop.) The only nontrivial ingredient are convolutions of $\mathcal{V}_{\eta}(k,\mu)$ with the one-loop $\Delta I_{qj}^{(1)}$ piece defined in eq. (B.2). Depending on the measurement, we distinguish the following four cases for the resulting Mellin kernel:

cumulant up to ${\mathcal{T}}_{\mathrm{cut}}>0$, cumulant up to $q_{T}^{\mathrm{cut}}>0$, $r\equiv(q_{T}^{\mathrm{cut}})^{2}/(Q_{i}{\mathcal{T}}_{\mathrm{cut}})$:

[TABLE]

where $B_{x}(a,b)$ is the incomplete Beta function,

[TABLE] 2. 2.

spectrum at ${\mathcal{T}}>0$, cumulant up to $q_{T}^{\mathrm{cut}}>0$, $r\equiv(q_{T}^{\mathrm{cut}})^{2}/(Q_{i}{\mathcal{T}})$:

[TABLE] 3. 3.

cumulant up to ${\mathcal{T}}_{\mathrm{cut}}$, spectrum at $q_{T}>0$, $r\equiv(q_{T})^{2}/(Q_{i}{\mathcal{T}}_{\mathrm{cut}})$:

[TABLE] 4. 4.

spectrum at ${\mathcal{T}}>0$, spectrum at $q_{T}>0$, $r\equiv(q_{T})^{2}/(Q_{i}{\mathcal{T}})$:

[TABLE]

In the first three cases the overall $\theta$ function cuts off the final PDF integral at

[TABLE]

and the Mellin kernel is regular up to and including $z_{\mathrm{cut}}$. In the last case we instead find a singularity at $z=z_{\mathrm{cut}}$, i.e., the subtraction from $\mathcal{V}_{\eta}$ now acts directly on the PDF integral. In either case we have exploited that terms proportional to $\delta(1-z)$ are cut off since $r>0$, so we could replace $\tilde{P}_{qj}^{(0)}$ by $P_{qj}^{(0)}$ throughout. We have also set $\mu={\mathcal{T}}\,({\mathcal{T}_{\mathrm{cut}}}$) in the boundary condition of $\mathcal{V}_{\eta}$ on the left-hand side, which can always be achieved by a shift in $\mu$. This ensures that the right-hand side depends only on the dimensionless parameters $r$ and $\eta$, up to an overall dimensionful prefactor. It is straightforward to check that for $\eta\to 0$ (at fixed order), the above results reduce to cumulants (spectra) of $\Delta I^{(1)}_{qj}$ itself.