UvA-DARE (Digital Academic Repository) Joint two-dimensional resummation in qT and 0-jettiness at

: We consider Drell-Yan production pp ! Z=(cid:13) (cid:3) ! ‘ + ‘ (cid:0) with the simultaneous measurement of the Z -boson transverse momentum q T and 0-jettiness T 0 . Since both observables resolve the initial-state QCD radiation, the double-di(cid:11)erential cross section in q T and T 0 contains Sudakov double logarithms of both q T =Q and T 0 =Q , where Q (cid:24) m Z is the dilepton invariant mass. We simultaneously resum the logarithms in q T and T 0 to next-to-next-to-leading logarithmic order (NNLL) matched to next-to-leading (cid:12)xed order (NLO). Our results provide the (cid:12)rst genuinely two-dimensional analytic Sudakov resummation for initial-state radiation. Integrating the resummed double-di(cid:11)erential spectrum with an appropriate scale choice over either T 0 or q T recovers the corresponding single-di(cid:11)erential resummation for the remaining variable. We discuss in detail the required e(cid:11)ective (cid:12)eld theory setups and their combination using two-dimensional resummation pro(cid:12)le 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 o(cid:11) 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.


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 -1 -

JHEP03(2019)124
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], Njettiness (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]. 1 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]. JHEP03(2019)124 Figure 1. The Drell-Yan cross section double-differential in the transverse momentum q T of the Z boson and the 0-jettiness event shape T at NNLL+NLO. For better visibility, the spectrum is plotted with respect to log 10 q T and log 10 T . On the two side walls we show the corresponding single-differential spectra in q T and T obtained by integrating the double-differential spectrum up to T cut = 100 GeV and q cut T = 100 GeV, respectively.
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 , 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.

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 → Z/γ * (→ + − ), 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 T (aka beam thrust) [38,46,47,53], defined as 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 where n µ a = (1, +ẑ) and n µ b = (1, −ẑ) are lightlike vectors along the beam axisẑ. 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 Our setup applies equally well to other definitions of 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 , T Q, where Q Λ 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 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 T via either of the projections onto q µ a and q µ b in eq. (2.1). This implies that their characteristic transverse momentum is ∼ T . Third, ISR with typical energy ∼ Q can contribute to T as long as it is collinear to either of the incoming beams, such that its contribution to T in eq. (2.1) is small. These collinear emissions then have a typical transverse momentum ∼ √ QT . The factorization and resummation structure of the cross section for q T , T Q depends on the parametric hierarchy between these scales. There are three relevant parametric regimes [51], which are illustrated in figure 2 and are discussed in the following.
In the first (blue) regime, T q T ∼ √ QT , soft emissions with transverse momentum ∼ T and collinear emissions with transverse momentum ∼ √ QT both contribute to the 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 SCET I . It has the same RG structure as the single-differential T spectrum, with q T acting as an auxiliary variable. The SCET I regime is discussed in more detail in section 2.2.
In the opposite (green) regime, T ∼ q T √ QT , both soft and collinear emissions have transverse momentum ∼ q T and thus contribute to q T . On the other hand, only soft radiation at central rapidities contributes to T , while the contribution from collinear radiation is suppressed. This regime is described by SCET II , whose RG structure is analogous to that of the single-differential q T spectrum, with T as the auxiliary variable. The SCET II regime is discussed in more detail in section 2.3.
Third, the intermediate (orange) regime in the bulk, T q T √ QT , shares features with both boundary cases. As in the SCET I regime, central soft radiation contributes to T , while as in the SCET II 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 SCET I and SCET II . 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 2) 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 α s (m Z ) = 0.118 and five active quark flavors.

SCET
In this regime, both soft and collinear modes are constrained by T , while only collinear modes can contribute to q T , whose characteristic transverse momentum √ QT coincides parametrically with q T . The scaling of the relevant EFT modes reads n a -collinear: p µ ∼ T , Q, QT ∼ q 2 T Q , Q, q T , n b -collinear: p µ ∼ Q, T , QT ∼ Q, q 2 T Q , q T , soft: p µ ∼ T , T , T , (2.4) in terms of lightcone coordinates defined by (with n ≡ n a ,n ≡ n b ) This leads to the following factorization formula for the cross section [46,59], The hard function H κ (Q, µ) describes the short-distance scattering that produces the lepton pair through the off-shell γ * or Z. In addition to Q, it depends on the partonic channel κ ≡ {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.6). The beam functions B q (t, x, k T , µ) describe extracting a quark (or antiquark) from the proton with momentum fraction x, virtuality t, and transverse momentum k T . The momentum fractions are directly related to Q and Y , (2.8) The t and k T encode the contribution of the collinear radiation to the T and q T measurement, as captured by the measurement δ functions on the last line of eq. (2.6). For t ∼ k 2 T Λ 2 QCD , these beam functions can be matched onto PDFs [46,59,60], The soft function S κ (k, µ) 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.6) separates the physics at the canonical SCET I scales By evaluating the ingredients at their natural scale and evolving them to a common scale, all logarithms of T /Q ∼ µ I S /µ I H ∼ (µ I B /µ I H ) 2 ∼ (µ I S /µ I B ) 2 are resummed. The hard and soft function in eq. (2.6) are the same as in the single-differential 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 SCET I resummation, with no large logarithms of q T appearing in the cross section as long as q T ∼ √ QT is satisfied. We stress that eq. (2.6) 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 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 T .
The factorization of the double-differential spectrum in eq. (2.6) (and in the following sections) does not account for effects from Glauber gluon exchange. scattering, they are expected to enter at O(α 4 s ) (N 4 LL ) [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 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 T ∼ q 2 T /Q 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, Here we abbreviated dσ ≡ dσ/(dQ dY dq T dT ), and dσ FO denotes the fixed-order cross section in full QCD. The scale subscripts on the right-hand side indicate whether dσ I is RG evolved using the SCET I resummation scales µ I , with their precise choices given below, or whether it is evaluated with all scales set to a common fixed-order scale µ FO . By construction, dσ I evaluated at common scales µ FO exactly reproduces the singular limit of dσ FO , such that the term in square brackets in eq. (2.11) is a pure nonsingular power correction at small T , which we can simply add to the resummed cross section. In the left panel of figure 3, we explicitly check that this is satisfied at fixed O(α s ), and numerically assess the size of the power corrections. We compare the full QCD result (solid orange) to the SCET I singular limit (dashed blue) as a function of T , while keeping q T = √ QT fixed -8 -

JHEP03(2019)124
to ensure that all classes of power corrections in eq. (2.7) uniformly vanish as T → 0. This is indeed satisfied, as the difference (dotted green) vanishes like a power. For T ∼ Q, the SCET I singular contribution and the power corrections are of the same size, implying that the resummation must be turned off to not upset the 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 µ I B ≡ µ I B (T ) and µ I S ≡ µ I S (T ) transition from their canonical values eq. (2.10) at small T to a common high scale for large T , schematically, As a result, the first and third term in eq. (2.11) exactly cancel in this limit, so the matched result reproduces dσ FO as desired.
For the concrete choices of µ I B , µ 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 SCET I -like jet vetoes in ref. [66] and used for the T resummation in Geneva [50]. The profile scales are chosen as with the profile function f I run given by [67] f (2.14) Based on figure 3, we take (x 1 , x 2 , x 3 ) = (0.2, 0.5, 0.8) for the transition points towards the fixed-order region x ∼ 1. In addition, eq. (2.14) turns off the resummation in the nonperturbative region x 2x 0 , where we set x 0 = 1 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 µ FO itself we use µ FO = Q as the central scale. Our central scale choices are illustrated as solid lines in the right panel of figure 3.
Perturbative uncertainties. We estimate perturbative uncertainties in dσ match I by considering two different sources. The first uncertainty contribution ∆ I is inherent to the SCET I resummation. It is estimated by varying the individual SCET I scales while keeping µ FO fixed, effectively probing the tower of higher-order logarithms that are being resummed. For this we use the profile scale variations [66]  where α = β = 0 corresponds to the central scale choice in eq. (2.13), and the variation factor is defined as (2.16) It approaches a factor of two in the resummation region at small x and reduces to unity toward the fixed-order regime at x = x 3 , where the resummation is turned off. The estimate for ∆ I is obtained by computing dσ I match for each of the four profile scale variations and taking the maximum absolute deviation from the central result. These variations are also indicated in the right panel of figure 3. 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, ∆ FO , we consider common variations of µ FO up and down by a factor of two in all pieces of eq. (2.11). Since µ FO enters all µ I scales as a common overall factor, they inherit the same variation, which keeps all resummed logarithms invariant. Hence, the µ FO variation effectively probes the effect of missing higher-order corrections in the fixed-order contributions. The final uncertainty estimate for dσ match I is obtained by adding both contributions in quadrature, The matched result dσ match I 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 ∼ √ QT .

SCET
In this regime, both soft and collinear emissions are constrained by q T . Only soft radiation is constrained by the T measurement, while collinear radiation at transverse momenta ∼ q T √ QT is not affected by it. The relevant EFT modes scale as n a -collinear: In this case, the cross section factorizes as [51] dσ -10 -

JHEP03(2019)124
The factorization receives power corrections of the form The hard function H κ (Q, µ) is the same as in eq. (2.6). In SCET II an additional regulator is required to handle rapidity divergences, for which we use the η regulator of refs. [68,69] as implemented at two loops in ref. [70], with ν the corresponding rapidity renormalization scale. The B q (ω, k T , µ, ν) are the same transverse momentum-dependent beam functions as in the single-differential q T spectrum. The large momentum components ω in eq. (2.20) are given by and we suppress the trivial dependence of the beam function on E cm . For k 2 T Λ 2 QCD , the beam functions satisfy a matching relation similar to eq. (2.9) [69,[71][72][73][74], The double-differential soft function S κ (k, k s , µ, ν) encodes the contribution of soft radiation to both T and q T . The RG consistency of the cross section implies that its µ and ν RGEs do not depend on T . Hence, the overall RG structure of the double-differential cross section is equivalent to the single-differential q T spectrum, with T acting as an auxiliary measurement. The factorization in eq. (2.20) separates the physics at the canonical SCET II invariantmass and rapidity scales 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 q 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 [72,[78][79][80]. The Fourier transform from q T to b T turns the vectorial convolutions in eq. (2.20) into simple products at b T ≡ | b T |. The canonical SCET II scales in b T -space are then given by 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 Q, eq. (2.25) 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 γ i ν , driving the ν running of the soft (or beam) function at fixed µ, is itself perturbatively renormalized at its intrinsic scale µ 0 and requires resummation when µ 0 = µ. Specifically, in the exponent of the b T -space rapidity evolution factor we havẽ where all logarithms of µ/µ 0 are resummed inside [See eq. (4.26) in ref. [77] for the analogous expressions in momentum space.] The canonical choice of µ 0 that eliminates all large logarithms in the fixed-order boundary conditioñ By choosing µ 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 ∼ 1/Λ QCD . 3 The mismatch to the full result can in principle be captured by a nonperturbative modelγ i ν,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γ i ν,np = 0 for simplicity. (We similarly ignore nonperturbative effects in the SCET II beam and soft function boundary conditions.) Our concrete choice of µ 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 ∼ T Q by an additive matching, Here the subscript µ II indicates that we evaluate dσ II at the SCET II resummation scales µ II (given below) in b T space, and take a numerical inverse Fourier transform in the end. The subscript µ FO indicates that it is instead evaluated at common fixed-order scales µ FO , which can be done directly in momentum space. Analogous to the discussion for SCET I , 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 4, which shows that the difference (green dotted) between the full QCD JHEP03(2019)124 result (solid orange) and the SCET II singular result (dashed blue) indeed vanishes like a power as q T → 0 along the line of fixed T = q T . Approaching q T ∼ T ∼ 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.24) as a function of the target q T value, schematically, We note that µ 0 does not need to asymptote to µ FO towards large q T because its effect on the matched result is already turned off as ν II S → ν II B . In this limit, the first and last term in eq. (2.29) exactly cancel, leaving the fixed-order result dσ 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 , (2.31) where f II run is a hybrid profile function given by f II run (x, y) = 1 + g run (x)(y − 1) .

JHEP03(2019)124
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 (2.33) As a result, for q T ≤ x 1 Q, the slope is unity yielding the canonical resummation, while for q T ≥ 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 4. We use the same transition points (x 1 , x 2 , x 3 ) = (0.2, 0.5, 0.8) as for SCET I , which is supported by figure 4. 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 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 µ 0 (or µ II S , or µ II B ) approaches the Landau where b 0 /b max Λ QCD ensures that all scales are canonical for small b T ≈ b * , but remain perturbative for large b T where b * → b max , as shown in the inset in the right panel of figure 4. In practice we pick b 0 /b max = 1 GeV , (2.35) in keeping with our choice of nonperturbative turn-off parameter in the SCET I case. The functional form of eq. (2.34) is the same as in the standard b * 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 * only affects the scales, so it essentially serves the same purpose as the x 0 nonperturbative cutoff in the SCET I scales in eq. (2.14). By contrast, the standard b * prescription corresponds to a global replacement of b T by b * , including the measurement itself. For the single-differential q T spectrum, this global replacement induces power corrections O(b 2 T /b 2 max ) 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 * prescription does in fact not work. This is because substituting b * for b T in the physical measurement renders Fourier integrals of the double-differential SCET II 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 -14 -

JHEP03(2019)124
only depend on x = b T T . Substituting b * for b T makes them asymptote to a constant for any given T , which upsets their required asymptotic behavior ∼ 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 T .
Perturbative uncertainties. To estimate the resummation uncertainty for dσ match II , we adopt the set of profile scale variations introduced for the SCET II -like jet veto in ref. [67]. They are given by , 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 SCET II regime. We add two independent variations of b 0 /b max = {0.5 GeV, 2 GeV} to probe the uncertainty in our nonperturbative prescription. The SCET II resummation uncertainty ∆ 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 SCET I , ∆ FO is estimated by overall variations of µ FO by a factor of two, which is inherited by all SCET II scales, so it probes the fixed-order uncertainties while leaving the resummed logarithms invariant. The total uncertainty estimate for dσ match II is then obtained as The matched result dσ match II in eq. (2.29) provides a prediction for the double-differential spectrum that covers the part of phase space where T ∼ 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 5 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 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.  Figure 5. The single-differential q T spectrum at NLL (blue) and NNLL+NLO (orange), using the q T resummation method described in the text. The bands indicate ∆ II ⊕ ∆ FO . In the right panel, the uncertainties are shown as percent differences relative to the central result at each order.

SCET
This regime is characterized by the presence of intermediate collinear-soft modes that contribute both to the q T and the T measurement, which uniquely fixes their scaling.
Central soft modes only contribute to T as in SCET I , while the energetic collinear modes only contribute to q T as in SCET II , The collinear-soft modes have the same virtuality as the collinear modes, p 2 ∼ q 2 T , but live at more central rapidity e |y| ∼ q T /T , which is small compared to the rapidity e |y| ∼ Q/q T of the collinear modes. Hence, the two have a SCET II -like relation and become a single collinear mode in the SCET I limit q T ∼ √ QT . At the same time, the collinear-soft and soft modes have a SCET I -like relation, being separated in virtuality, and become a single soft mode in the SCET II limit T ∼ q T . In this way, SCET + is able to connect the SCET I and SCET II regimes. This is similar to the collinear-soft mode originally introduced in ref. [28], which instead connected two SCET I theories.
The cross section in SCET + factorizes as [51] dσ which holds up to power corrections The hard function is the same as before. The beam functions are the q T -dependent ones from SCET II , while the soft function is the T -dependent one from SCET I . The new ingredient is the double-differential collinear-soft function S κ (k, k T , µ, ν), which encodes the contributions of the collinear-soft modes to both q T and 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 SCET I and SCET II cross sections [51], which precisely reflects the relation between the EFT modes described above. Expanding the SCET I double-differential beam function in the limit q T √ QT , it factorizes into the SCET II beam function and the collinear-soft function, The ν dependence of the two terms on the right-hand side must cancel, while their µ 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.18). Similarly, expanding the SCET II double-differential soft function in the limit T q T , it factorizes into the SCET I soft function and the two n a -collinear-soft and n b -collinear-soft functions, Since the left-hand side does not depend on ω 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 ω a ω b = Q a Q b = Q 2 . In addition, boost invariance at the level of the collinear-soft matrix element implies that d + a S κ ( + a , k T , µ, ν) can only depend on the product + a ν (and analogously for − b ). 4 Hence, we can rewrite

JHEP03(2019)124
where in the first step we changed variables from ± a,b to k + a = ω a In the second step we performed the rapidity evolution from ν a,b ≡ Q a,b ν/ω a,b back to a common ν at fixed µ [see eq. (B.46)], for which the rapidity evolution factors exactly cancel because ln The SCET + factorization in eq. (2.39) fully disentangles the physics at the following SCET + canonical energy and rapidity scales: As for SCET II , 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 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., The logarithms of the first ratio appear in the double-differential SCET I beam function in the limit q T √ QT , and are resummed in SCET + by the additional ν evolution in the refactorization in eq. (2.41). Similarly, logarithms of the second ratio appear in the doubledifferential SCET II soft function in the limit T q T , and are resummed in SCET + by the additional µ evolution in eq. (2.42). 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.

Outer space
We now briefly discuss the outer phase-space regions left blank in figure 2. The region above the SCET II regime is characterized by the hierarchy q T T √ QT , while the region to the right of the SCET I regime corresponds to T √ QT q T . Both regions are power suppressed.
As we have discussed in section 2.3, only the soft function contributes to T in SCET II , as the collinear contribution is power suppressed. However, for q T T , even the soft contribution to T becomes power suppressed. In particular, for a single real emission at fixed O(α s ), the region T > q T is kinematically forbidden both in SCET II 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 O(q 2 T /T 2 ). Equivalently, expanding the SCET II factorization of the double-differential cross section in the limit q T T reduces it to the single-differential q T spectrum with an overall δ(T ), which we exploit in our numerical implementation, cf. eq. (B.33). Physically this means that by integrating the double spectrum in SCET II up to some T cut 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 O(Λ 2 QCD /T 2 ), with T setting the scale of the second hard scatter producing the back-to-back jets.
Similarly, in the limit √ QT q T , even the contribution from collinear radiation to q T becomes power suppressed in SCET I [cf. eq. (B.24)], and at leading power we recover the single-differential T spectrum with an overall δ(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 √ QT q T is that very energetic forward radiation with energy ∼ q 2 T /T can theoretically contribute [51], pushing the hard scale up to q 2 T /T 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 2, corresponding to the main SCET I , SCET II , and SCET + regimes.

Structure of power corrections
An important feature of our EFT setup is that the factorized cross section in SCET + differs from the ones in SCET I and SCET II only by a subset of the power corrections it receives relative to the full QCD result, This is illustrated in figure 6, 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 σ I and σ + share a logarithmic singularity with respect to 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 µ + scales to be equal to their µ I counterparts, such that any large logarithms inside the refactorized beam function in eq. (2.41) are treated at fixed order. We write dσ + µ I to indicate that dσ + 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 -19 -JHEP03(2019)124 Figure 6. Venn diagram of power corrections to the factorized double-differential spectrum. SCET I (blue) and SCET II (green) each capture a set of power corrections that is expanded away in the SCET + factorization (red) and the opposite boundary regime. A third class of power corrections to the overall soft-collinear limit is captured by the fixed-order calculation in full QCD (gray).
compare dσ + µ I to dσ I µ I , with our choices for µ I as given in section 2.2, i.e., including the whole set of all-order terms from the T resummation in both of them. This comparison is shown in figure 7 for representative choices of fixed T and q T at NNLL. We can clearly read off a power-like behavior of the difference dσ I − dσ + µ I (dotted green) as either q T → 0 for fixed T (left panel) or T → ∞ for fixed q T (right panel). This also provides a nontrivial check on our implementation of σ I and σ + . This comparison in figure 7 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 µ FO . Here, SCET I takes on the role of the full theory, while SCET + provides the singular limit, and the comparison is performed at common scales µ I . Similarly, both σ II and σ + have a common singular structure as q T /Q → 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 µ + scales to be equal to their µ II counterparts, which treats the large logarithms in the refactorized double-differential soft function in eq. (2.42) at fixed order. We denote this choice of scales by dσ + µ II , with scale setting in b T space and the inverse Fourier transform understood as in section 2.3. In figure 8 we compare dσ + µ II to dσ II µ II at NNLL as a function of T at fixed q T (left) and vice versa (right). It is clear that even when evaluated at its intrinsic scales, dσ II µ II (solid orange) exhibits an unresummed singularity as T /q T 1, which, as expected, is captured by dσ + µ 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 T q T is correctly captured by SCET II , 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 6, the complete infrared structure of the double-differential spectrum for q T Q and T Q, i.e., for any hierarchy between q T and T , is described by adding the SCET I and SCET II cross sections and removing the overlap between the two by subtracting the SCET + cross section, In figure 9 we numerically check this relation at fixed O(α s ), which requires setting all scales equal to a common µ FO . We plot the comparison as a function of q T along lines of fixed T /q T = 1/2 (left) and 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.4) (dashed blue), as evident from the power-like behavior of their difference (dotted green) as q T , T → 0. This singular/nonsingular comparison is qualitatively different from the structure of power corrections in either SCET I or SCET II alone, which we already verified in figure 3 and figure 4. Because SCET I and SCET II both involve an additional expansion about a specific hierarchy between q T and T , they incur power corrections O(T 2 /q 2 T ) or O(q 2 T /(QT )), respectively. Accordingly, they only recover the singular limit of full QCD when approaching it along specific lines in the (q T , T ) plane. This is different from figure 9, where the combined expression in eq. (3.4) (dashed blue) describes the singular limit q T , T → 0 along an arbitrary line of approach, with the ratio q T /T effectively controlling the "admixture" of power corrections O(q 2 T /Q 2 ) and O(T /Q), respectively. We have verified that also for other fixed ratios of q T and 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 dσ I +dσ II −dσ + for use as a double-differential subtraction term to treat infrared divergences in fixed-order calculations.

Matching formula
The structure of power corrections discussed in the previous section, together with the allorder resummation shared between SCET + and SCET I or SCET II , suggests the following matching formula to describe all regions of the double-differential spectrum: The only ingredient in this matching formula we have not yet discussed is dσ + µ + , for which all ingredients in the SCET + factorization are evaluated at the SCET + scales µ + , such that the full RGE of SCET + is in effect. In the following we describe the requirements on µ + to ensure the best possible prediction across phase space. Our precise construction of µ + 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, µ + is given by the canonical SCET + scales in eq. (2.46). As for µ II , these scales are set in b T space, followed by an inverse Fourier transform.
As we approach the SCET I region, the resummation inside the refactorization of the beam function in eq. (2.41) must be turned off, In addition we can identify the soft scales in SCET I and SCET + because the soft functions are identical, These relations must hold for every value of the b T argument of the scale. Similarly, as we approach the SCET II region, the scales inside the refactorized soft function eq. (2.42) must become equal and we can identify the scales of the common beam function in SCET II and SCET + , 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 SCET I are -23 -

JHEP03(2019)124
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 which does not coincide with any scale on either boundary. It is instructive to explicitly consider the behavior of eq. (3.5) on the SCET I and SCET II phase-space boundaries, as well as in the fixed-order region. By construction, for any choice of µ + scales satisfying eqs. (3.6) and (3.7) we have (3.11) It follows that This mostly coincides with the result in eq. (2.11) of matching dσ I to the fixed-order result dσ FO , and is guaranteed to capture all large logarithms of T /Q captured by the SCET I RGE. It improves over eq. (2.11) by also resumming logarithms of q T /Q in the power corrections O(T 2 /q 2 T ), encoded in σ II − σ + µ II . This is not a numerically large effect and cannot be exploited to achieve the resummation of 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 and consequently Finally, in the fixed-order region, all µ + , µ I , and µ II scales become equal to µ FO . Thus as desired, the matched prediction reduces to the fixed-order result,

Profile scales
In this section, we describe our choice of the central µ + scales for the various ingredients in the SCET + factorized cross section, taking into account the transition to the SCET I and -24 -JHEP03(2019)124 SCET II 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 SCET I and SCET II .
We start by summarizing the canonical scales for SCET I , SCET II , 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 Its definition is carefully chosen such that a = 1 when the SCET I parametric relation is exactly satisfied, q T = √ T Q, and a = 2 on the SCET II boundary of phase space, q T = T . As illustrated in the left panel of figure 10, the canonical SCET + region lies at intermediate a ∼ 1.5. The requirements on the SCET + scales were given in eqs. (3.6) and (3.7) for the transition to SCET I , and in eqs. (3.8) and (3.9) for SCET II . To satisfy these requirements, we take weighted products of scales on the boundary and in the bulk, schematically, (3.17) The weights in the exponent are given by helper functions that depend on a, as illustrated in the right panel of figure 10. They satisfy for any a and are given explicitly in eq. (3.22) below. The helper functions ensure that the appropriate scales are used in each region, e.g., h 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: The most nontrivial of these cases is ν S , which must be equal to the overall ν in the SCET I 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 SCET II boundary. We note that µ + S also requires a distinct treatment in the bulk to ensure that the hierarchy µ + S < µ + 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 The profile function f I run was introduced for the transition between SCET I and fixedorder QCD in eq. (2.14), and similarly for the hybrid profile f II run in eq. (2.32) and the nonperturbative b * (b T ) prescription in eq. (2.34). These functions turn off the resummation of logarithms involving q T (b T ) and 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 -26 -

JHEP03(2019)124
bulk scales all assume their canonical values for q T , T Q as given in table 1, and asymptote to µ FO when simultaneously taking q T , T → Q. The beam function scales in the bulk can simply be associated with their SCET II counterparts and only require a transition towards the SCET I boundary, In our numerical implementation, we choose the helper functions h I,II,+ as where the polynomials governing the interpolation between zero and one are . (3.23) The transition points a 1,...,6 determine the transition between the different regions, as can be seen from the helper functions in figure 10: for values a 3 ≤ a < a 4 , the exact canonical SCET + scales are selected, implying that the resummation of logarithms of both q T and T is fully turned on. For lower values a 1 ≤ a < a 3 , the additional q T resummation is smoothly turned off and for a < a 1 , SCET I scales are used so that only logarithms of T are resummed. Conversely, for higher values of the regime parameter a 4 ≤ a < a 6 , the resummation of T logarithms is smoothly turned off. At values a 6 ≤ a, SCET II scales are selected by the helper functions, and the additional resummation of logarithms of 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 ≥ a 6 = 2, we fully recover SCET II resummation and faithfully describe the kinematic edge at q T ∼ T by preserving the O(1) cancellation between σ + µ II and the SCET II nonsingular contribution visible at a ∼ 2 in the left panel of figure 8. (In both figures 7 and 8, corresponding values of a are indicated on the horizontal axis at the top of the panels.) On the other hand, from figure 7 we observe that power corrections from SCET I are smaller and tend to set in at values of a lower than the naively expected a = 1. E.g., an O(1) cancellation between σ + µ I and the SCET I nonsingular only is in effect around a ∼ 0.5 in the right panel of figure 7, leaving more room for slowly turning off the SCET + resummation down towards a 1 = 0.5. This is expected because the SCET I nonsingular encodes the suppression of collinear recoil beyond the naive phase-space boundary at a ∼ 1 (q T ∼ √ QT ) that is washed out by the PDFs, unlike the sharp kinematic edge at q T ∼ T encoded in the SCET II nonsingular. For simplicity we set a 3 = a 4 for our central prediction,

JHEP03(2019)124
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.

Perturbative uncertainties
In this section we describe how we assess perturbative uncertainties by varying the scales entering the matched prediction in eq. (3.5). Following the same approach as for SCET I and SCET II on their own (see sections 2.2 and 2.3), we distinguish between resummation uncertainties and a fixed-order uncertainty. The fixed-order uncertainty ∆ FO is estimated by varying µ FO up and down by a factor of two, i.e., by setting µ FO = {Q/2, 2Q}. Since all scales (in any piece of the matching formula) include an overall factor of µ FO , the ratios between the various scales remain unchanged and the same logarithms are resummed. The fixed-order uncertainty ∆ 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.5). To probe the tower of logarithms of T /Q predicted by the SCET I RGE, we perform variations of µ I B and µ I S parametrized by α and β as in eq. (2.15). This directly affects the resummed power corrections dσ I − dσ + µ I captured by SCET I . In addition, however, dσ + µ + near the SCET I boundary also undergoes variations because for large h I , the SCET + scales in eqs. (3.19) and (3.21) strongly depend on their SCET I counterparts and inherit their variations. Our setup thus ensures that in (or near) the SCET I region, variations probing resummed logarithms of T /Q are properly treated as correlated between the SCET + cross section and the SCET I matching correction. When referring to the matched prediction in eq. (3.5), we take ∆ I to be the maximum deviation of dσ match from its central value under these correlated variations of α, β.
In complete analogy, we define ∆ II as the maximum deviation under correlated variations of µ II as described in section 2.3. These variations act on both dσ II − dσ + µ II and dσ + µ + , where now the SCET + scales inherit variations from µ II near the SCET II boundary (where h II is large). As a result, ∆ II probes an all-order set of logarithms of (b 0 /b T )/Q predicted and resummed by the SCET II 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 ), we perform the following two variations of the SCET + bulk scales, where γ = 0 corresponds to the central scales in eq. (3.20). Similarly to the role of β in the SCET I variations [see eq. (2.15)], making the strength of the γ variations depend on the ratio q T /T ensures that the hierarchy µ S < µ S implied by the SCET + power counting is not upset by variations, counting b 0 /b T ∼ q T . We note that the third independent bulk scale ν + S,bulk does not require independent variation because it only enters through rapidity logarithms of ν + B /ν + S , which are already being probed by variations of ν + B inherited from SCET II . 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 ∆ + . The total uncertainty assigned to the matched prediction is then given by adding all contributions in quadrature, (3.26)

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 T as an example. There are two equivalent quantities of interest in this case, namely the spectrum dσ/dT with respect to T and the cumulant cross section σ(T cut ) with a cut on T . The two quantities are related by 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 T to directly predict the spectrum, or in terms of the cumulant variable T cut to predict the cross section integrated up to T 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 In the first term under the integral in the cumulant cross section, all scales µ entering the resummed and matched prediction depend on the integration variable T . Because our setup only reliably predicts the spectrum away from the nonperturbative region, we choose -29 -

JHEP03(2019)124
to integrate the resummed spectrum with differential scale setting up from some small cutoff T np , and include an "underflow" contribution given by the second term under the integral. For the underflow contribution for T ≤ T np , the spectrum is evaluated at fixed scales corresponding to T 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 In this case, the scales in the cumulant cross section depend on T cut , and not the integration variable T , so the integral up to T 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 µ i in µ with respect to T . Cumulant scale setting ensures that for T cut → 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, Thus for cumulant scale setting, the spectrum has the correct (fixed-order) normalization. However, the additional derivatives of the scales in eq. (3.29) tend to produce artifacts in the spectrum if the profile functions µ i (T ) used to interpolate between the resummation region T Q to the fixed-order region T ∼ 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. [50,65,91,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 -30 -

JHEP03(2019)124
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 T , there are in principle four quantities of interest, namely the double-differential spectrum dσ/dq T dT , the single-differential spectra dσ(q cut T )/dT and dσ(T cut )/dq T with a cut on the other variable, and the double cumulant σ(q cut T , T 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 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 (T ) replaced by q cut T (T 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., σ(q cut T , T cut ) will by construction recover the inclusive fixedorder cross section when lifting both cuts, while dσ(q cut T )/dT and dσ(T cut )/dq 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 , 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 (T ) to a prediction for the cumulant up to q cut T (T cut ) using the analogue of eq. (3.28). The only nontrivial new procedure is computing the double cumulant directly from (q T , T ) scales, where we need to account for an overlap in underflow contributions as . The distinction between differential or cumulant scale setting is only relevant for q T versus q cut T 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 np T = T np = 1 GeV, and implement the integrals in eqs. In figures 11 to 13, 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 µ in terms of q T rather than q cut T and/or T rather than T cut . In figure 11, we show the double cumulant cross section, for which our default is to use scales in terms of q cut T and T cut . The horizontal reference line indicates the inclusive fixed-order cross section. In figure 12 we show the T spectrum with a cut on q T , for which our default scales are in terms of q cut T and T , and the converse for figure 13. In figures 12 and 13 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 13 slightly overshoots the inclusive result beyond the phase-space boundary T cut 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 cut T scale setting are mutually compatible, i.e., their uncertainty bands (very nearly) overlap, as long as the scale setting with respect to T is done the same way in both cases. This can be seen from the right panel of figure 11 by contrasting the default (q cut T , T cut ) scales (solid orange) and (q T , T cut ) scales (dotted green). Similarly, in figure 12 we find that the default (q cut T , T ) scales (solid orange)   Figure 13. The q T spectrum with a cut on T as a function of T cut for q T = 15 GeV (left) and as a function of q T for T cut = 100 GeV (right). The bands indicate the total perturbative uncertainty ∆ total , see section 3.4. The colors correspond to different scale setting prescriptions (default: solid orange); see the text for detail. and (q T , 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 cut T 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 run is nonzero. E.g. for a high T cut = 100 GeV, which is also a good proxy for the inclusive q T spectrum, the two prescriptions fully agree in the canonical region q cut T ≤ 20 GeV (see the right panel of figure 11). This is responsible for the good overall agreement because most of the cross section is concentrated in the canonical region.

JHEP03(2019)124
The comparison of T versus T 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 +NLO). The mismatch is most striking between the default scales (solid orange) and (q T , T ) scales (dashed blue) in figures 11 and 13, 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 T alone, so that the methods developed for the single-differential case in refs. [50,92] 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 µ I S ∼ (µ I B ) 2 /µ I H ∼ T in SCET I 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 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 14, where we repeat the left panels of figures 11 and 13 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 cut T , T ) scales (dotted green, left panel) covering the inclusive fixed-order cross section and the result from (q T , 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 T , e.g. if large bins in either observable are considered, the double-differential profile setup given in this paper, using (q cut T , T ) or (q T , T cut ) scales as appropriate, will be completely sufficient.

Results
In this section we present our results for Drell-Yan production pp → Z/γ * → + − at the LHC, with a simultaneous measurement of the transverse momentum q T of the lepton pair and the 0-jettiness event shape T . The center-of-mass energy is taken to be E cm = 13 TeV. We assume that in addition, the invariant mass Q of the lepton pair is measured, and write pp → Z for short if Q = m Z , and pp → Z * otherwise. The subsequent decay and the contribution from the virtual photon are included in either case.
To obtain numerical results for the SCET I , SCET II , and SCET + contributions, we have implemented all pieces of the relevant double-differential factorized cross sections to O(α 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 α 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 GeV in q T and 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 , T ) plane and for selected slices along lines of constant q T or T . We demonstrate that our prediction smoothly interpolates between the SCET I and SCET II boundary theories, i.e., we show that our matching formula in eq. (3.5) 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 dσ(q cut T )/dT and dσ(T cut )/dq T with a cut on the other variable, and show how they recover the inclusive single-differential T and q T spectra for large values of q cut T and T cut , respectively.

Double spectrum and comparison with boundary theories
To highlight the necessity of the simultaneous resummation of large logarithms of both q T and T , we start by showing results for the double spectrum (the cross section doubledifferential in q T and T ) where only some of the logarithms are resummed. These results are shown as surface plots in figure 15, where we plot the double-differential spectrum with respect to log 10 q T and log 10 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 T cut = 100 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 T , with a cut at q cut T = 100 GeV. 5 The top left panel of figure 15 shows the spectrum evaluated at fixed O(α s ), without any resummation. The double-differential fixed-order spectrum diverges logarithmically for small T at any value of q T , while its projections onto the single-differential spectra in q T and T feature double-logarithmic singularities. Notably, the double-differential spectrum has a sharp kinematic edge along q T = 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 η, which contributes at most T = k T e −|η| ≤ k T = q T . Due to the vectorial nature of q T , however, back-to-back emissions can populate the region 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 15, we show the result of resumming logarithms of T using the SCET I matched result in eq. (2.11). The resummation is performed at NNLL and is matched to full NLO, which we refer to as NNLL T +NLO. As discussed in section 2.2, this prediction is valid as long as the parametric relation T q T ∼ √ QT is satisfied. This corresponds to the SCET I phase-space boundary (blue) in figure 2, running from the region of small T and intermediate q T towards the fixed-order region where q T ∼ T ∼ Q. It is clear that away from its region of validity, the NNLL T +NLO result contains unresummed logarithms of q T because at any point in T , the prediction diverges for very small q T . In particular, power corrections of O(T 2 /q 2 T ) are only captured by the fixed-order matching. They become O(1) as one approaches the diagonal T = q T (the green line in figure 2), and encode the phase-space boundary at q T ∼ 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 SCET I cannot resolve the boundary because it arises from the dynamics at central rapidities. The projections onto the rear walls highlight that only T is resummed. The single-differential q T spectrum still diverges as q T → 0, while the T spectrum features a physical Sudakov peak.
In the bottom panel of figure 15, 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 SCET II matched result in eq. (2.29). We denote this order by NNLL q T +NLO. This result is valid for T ∼ q T √ QT , i.e., around the SCET II phase-space boundary (green) in figure 2, where we find the onset of a Sudakov peak from the q T resummation and a smooth kinematic 5 We refer the reader to section 3.5 for the precise way we perform these integrals.  Figure 16. The double-differential Drell-Yan cross section at NNLL+NLO, at Q = m Z (top) and Q = 300 GeV (bottom), with respect to log 10 q T and log 10 T . On the rear walls we show the result of integrating the double spectrum over either variable up to a cut at 100 GeV. The contour plots indicate total perturbative uncertainties relative to the cross section, ∆ total = ∆ + ⊕∆ I ⊕∆ II ⊕∆ FO . The contour plots are left blank in the region where dσ/(dQ d log 10 q T d log 10 T ) is less than 3% of its peak height.
suppression towards T q T . However, the NNLL q T +NLO result diverges for smaller values of T . This is due to unresummed logarithms of T in both the factorized cross section in SCET II and terms of O(q 2 T /(QT )) 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 T .
Our final results for the Drell-Yan double spectrum are shown in figure 16, as given by the fully matched prediction in eq. (3.5). 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 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 2), is smoothly suppressed beyond, and shifts towards higher values of q T and T for Q = 300 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 T +NLO results in figure 15 over T and q T , respectively. The contour plots in figure 16 show the total perturbative uncertainties ∆ total as percent deviations from the central result for the double spectrum. As described in section 3.4, ∆ total combines estimates of all sources of resummation uncertainty in the prediction.
In figure 17, we break down the uncertainty for the Drell-Yan double-differential spectrum at Q = m Z into its contributions from SCET I , SCET II and SCET + resummation uncertainties, respectively. As expected, the SCET I resummation uncertainty dominates in the SCET I region of phase space, and similarly for SCET II . 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 NNLL T , we take slices of the surface plots and overlay them in figure 18, keeping q T (left) or T (right) fixed. The solid red curve corresponds to the matched and fully resummed cross section in eq. (3.5), with the uncertainty band given by the total perturbative uncertainty ∆ total , see eq. by ∆ I total and ∆ II total , which only probe a subset of higher-order terms as predicted by the respective RGE, see eqs. (2.18) and (2.37). The SCET I prediction features an unphysical sharp edge at T = q T , cf. the middle panel of figure 15, and for this reason is cut off at T = 0.8 q T .
All panels in figure 18 show that our final prediction smoothly interpolates between the SCET I and SCET II boundary theories, both for the central values and for the uncertainties. Specifically, the matched prediction tends towards SCET I (SCET II ) for small (large) values of T and large (small) values of q T . In the left column one clearly sees that SCET II only captures logarithms of T at fixed order, leading to a diverging spectrum as T → 0, while the complete NNLL result features a physical Sudakov peak. Conversely, the SCET I result diverges as q T → 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  Figure 19. Slices of the double-differential Drell-Yan cross section at q T = 15 GeV (left) and T = 5 GeV (right). The solid red, dashed blue, and dotted green curves are identical to the central results in figure 18. The solid blue and green curves depict the SCET I and SCET II limits of our fully resummed result, given in eqs. (3.12) and (3.14). The thin vertical lines indicate the transition points a i described in section 3.3. intrinsic SCET + resummation is turned off. The reason for this is that even in these limits, the matched cross section in eq. (3.5) improves over the matched SCET I and SCET II cross sections in eq. (2.11) and eq. (2.29) by an additional resummation of power-suppressed terms, cf. eqs. (3.12) and (3.14). 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 19, but for reference include the case where σ + in the matched prediction is evaluated at µ I (solid blue) or µ 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 SCET I 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.

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 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 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 20, we show some more detailed results for the single-differential spectra with an additional cut, where the left panel shows dσ(q cut T )/dT as a function of T for various values of q cut T , and the right panel shows dσ(T cut )/dq T as a function of q T for various values of T cut . By increasing the value of the cut, they can be seen to approach the inclusive single-differential spectra (orange solid),  Figure 20. The single-differential T (left) and q T (right) spectrum with a cut on the other variable at NNLL+NLO. The different curves represent different values of the cut. The solid orange lines correspond to the inclusive single-differential spectrum obtained by lifting the cut. with which they must agree when sending q cut T → ∞ or T cut → ∞, 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 10 GeV push the peak to lower values and suppress the tail, where the q T spectrum is somewhat more resilient to cuts on T than vice versa. Intermediate cuts ∼ 10-15 GeV keep the peak and mostly lead to a suppression in the tail, while the effect of cuts 40 GeV is almost negligible in the q T and T ranges of interest.

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 T . Both T and q T probe the initial state radiation, leading to Sudakov double logarithms of 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 nextto-leading fixed order. We accomplish this resummation by using SCET I and SCET II to describe the regions T q T ∼ √ T Q and T ∼ q T √ 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 SCET II . 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 SCET I and SCET II boundaries. that scales and scale variations are still, to the extent possible, inherited from the singledifferential resummation of 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 SCET I , SCET + and SCET II , we introduced profile scales in terms of a regime parameter a, designed such that a = 1 for SCET I and a = 2 for SCET II . 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 T and q T yields the correct resummed single-differential cross sections after integrating over either 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 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 T , and would benefit from the simultaneous resummation of q T . Indeed, our NNLL results clearly indicate that only resumming the logarithms of either 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.

A Plus distributions and Fourier transform
We use the following standard plus distributions with dimensionless arguments,

JHEP03(2019)124
In intermediate steps we need a two-dimensional plus distribution originally defined in ref. [51], Our shorthands for distributions with dimensionful arguments in one spatial dimension are In terms of the second class of (power-like) plus distributions, we further define which have a group property, assuming identical boundary condition µ, Shifting the boundary condition of V a (k, µ) from µ to µ is also straightforward, We use the conventions from appendix C of ref. [77] for logarithmic plus distributions in two integer spatial dimensions, with k 2 Our convention for the two-dimensional Fourier transform also follows ref. [77], Here we make the mass dimension of df /d p T = df /(dp x dp y ) explicit. If f is azimuthally symmetric, i.e., if for the azimuthal integral can be performed, leaving 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

JHEP03(2019)124
where J 1 (x) is the first-order Bessel function of the first kind. Fourier transforms of L n ( p T , µ) can be found in table 5 of ref. [77], and are most conveniently expressed in terms of We expand the β function of QCD as The coefficients in the MS scheme are, up to three loops [93,94], We work with n f = 5 light flavors. The cusp and noncusp anomalous dimensions are expanded as 3) The coefficients of the MS cusp anomalous dimension to three loops are [95][96][97] 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, where we have already used that γ i ν 0 = 0. For our choice of regulator, the two-loop boundary condition is given by [70]  . The hard matching coefficient for qq → Z/γ * is renormalized as The coefficients of the quark noncusp anomalous dimension up to two loops are Beam functions and PDFs. In b T space, the TMD beam function is renormalized as We include a tilde to indicate thatγ q B is related to the SCET II beam function, even though it does not depend on b T . Its coefficients through two loops are [70] γ q B 0 = 6C F , γ q B 1 = C F (2 − 24ζ 3 )C A + (3 − 4π 2 + 48ζ 3 )C F + 1 + The resummed rapidity anomalous dimensionγ i ν (b T , µ) is given in eq. (2.26). The double-differential beam function satisfies the same RGE as the inclusive SCET I beam function, 6 µ d dµ B q (t, x, k T , µ) = dt γ q B (t − t , µ) B q (t , x, k T , µ) , γ q B (t, µ) = −2Γ q cusp (α s ) L 0 (t, µ 2 ) + γ q B [α s (µ)] δ(t) . (B.11) The coefficients of the SCET I quark beam anomalous dimension are [60,98] γ q B 0 = 6C F , γ q B 1 = C F 146 9 − 80ζ 3 C A + (3 − 4π 2 + 48ζ 3 )C F + 121 9 + 2π 2 3 β 0 . (B.12) 6 We note that ref. [51] incorrectly did not distinguish between γ i B (αs) andγ i B (α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.18) and the nonvanishing two-loop noncusp coefficient in our eq. (B.19).

JHEP03(2019)124
We also require the one-loop coefficients of the PDF anomalous dimension, (B.13) Note that we expand P ij (α s , z) in α s /(4π). The one-loop coefficients are P (0) qq (z) = 2C F θ(z)P qq (z) , P (0) qg (z) = 2T F θ(z)P qg (z) , (B.14) in terms of the standard color-stripped one-loop QCD splitting functions where we again use a tilde on the µ anomalous dimension to indicate that it relates to the SCET II soft function. The RGE of the collinear-soft function in b T space reads µ d dµS i (k, b T , µ, ν) = dk γ i S (k − k , µ, ν)S i (k , b T , µ, ν) , The soft and collinear-soft noncusp anomalous dimension coefficients are only nonzero starting at two loops and can be inferred from consistency, Hard process. The Born cross section for qq → Z/γ * → + − is given by where Q q is the quark charge in units of |e|, v ,q and a ,q are the standard vector and axial couplings of the leptons and quarks, and m Z and Γ Z are the mass and width of the Z boson. The one-loop Wilson coefficient C V,A qq (Q 2 , µ) from matching the quark current in QCD onto SCET was computed in refs. [99,100]. This leads to the following hard function [46], δ iq δ jq +δ iq δ jq |C V,A qq (Q 2 , µ)| 2 , where Re denotes the real part.
Beam functions. The one-loop matching coefficients for the single-differential SCET II beam function in b T space are given by [70,101] I qj (ω, b T , z, µ, ν) (B.22) As for its anomalous dimension, we use a tilde to indicate that these boundary conditions are part of the SCET 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 I quark beam function as I qj (t, z, k T , µ) = δ( k T ) I qj (t, z, µ) + ∆I qj (t, z, k T , µ) , (B.24) where I qj (t, z, µ) is the matching coefficient for the inclusive quark beam function [46,60], with the finite terms in this case given by 2T F (1 − z) 2 + z 2 , j = g .
The ∆I qj piece can be interpreted as a correction over the limit t k 2 T , where recoil from collinear radiation is power suppressed and the double-differential beam function becomes proportional to δ( k T ). Specifically, it scales as ∆I qj (t, z, k T , µ) and by construction satisfies d 2 k T ∆I qj (t, z, k T , µ) = 0 . (B.28) At one loop it can be extracted from the full calculation of I qj (t, z, k T , µ) [59,102] and has the compact form ∆I qj (t, z, k T , µ) = α s (µ) 4π ∆I (1) qj (t, z, k T ) + O(α 2 s ) ,

∆I
(1) The second line is regular in t because the term in square brackets vanishes as t → 0. After accumulating over the transverse plane up to q cut T > 0, we have d 2 k T θ(q cut T − | k T |) ∆I The one-loop collinear-soft function in b T space is [51] S i (k, b T , µ, ν) = δ(k) + α s (µ) 4π It is again convenient to decompose the double-differential soft function computed in ref. [51] into separate pieces with distinct power counting, S i (k, k T , µ, ν) = δ(k) S i ( k T , µ, ν) + ∆S i (k, k T , µ, ν) . (B.33) Here S i ( k T , µ, ν) is the standard single-differential q T soft function, which in b T space at one loop is given by [69] S i (b T , µ, ν) = 1 + α s (µ) 4π − Γ i where B x (a, b) is the incomplete Beta function, 2. spectrum at T > 0, cumulant up to q cut T > 0, r ≡ (q cut T ) 2 /(Q i T ): 3. cumulant up to T cut , spectrum at q T > 0, r ≡ (q T ) 2 /(Q i T cut ): (B.51) 4. spectrum at T > 0, spectrum at q T > 0, r ≡ (q T ) 2 /(Q i T ): In the first three cases the overall θ function cuts off the final PDF integral at z < z cut ≡ 1 1 + r , (B.53) -52 -

JHEP03(2019)124
and the Mellin kernel is regular up to and including z cut . In the last case we instead find a singularity at z = z cut , i.e., the subtraction from V η now acts directly on the PDF integral.
In either case we have exploited that terms proportional to δ(1 − z) are cut off since r > 0, so we could replaceP (0) qj by P (0) qj throughout. We have also set µ = T (T cut ) in the boundary condition of V η on the left-hand side, which can always be achieved by a shift in µ. This ensures that the right-hand side depends only on the dimensionless parameters r and η, up to an overall dimensionful prefactor. It is straightforward to check that for η → 0 (at fixed order), the above results reduce to cumulants (spectra) of ∆I (1) qj itself.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.