Event-Based Transverse Momentum Resummation

TL;DR

This paper introduces an automated framework for precise transverse momentum resummation in electroweak processes, enabling detailed analysis of small transverse momentum observables with high accuracy.

Contribution

It presents a fully differential, automated resummation method at NNLL accuracy matched to NLO fixed order, integrated into MadGraph5_aMC@NLO for electroweak final states.

Findings

Successful implementation for Z and W boson production

Good agreement with experimental data for transverse momentum distributions

Enhanced analysis capabilities for small transverse momentum observables

Abstract

We have developed a framework for automated transverse momentum resummation for arbitrary electroweak final states based on reweighting tree-level events. It is fully differential in the kinematics of the electroweak final states, which facilitates a straightforward analysis of arbitrary observables in the small transverse momentum region. We have implemented the resummation at next-to-next-to-leading logarithmic accuracy and match to next-to-leading fixed-order results using the event generator MadGraph5_aMC@NLO. Results for $Z$ and $W$ boson production with leptonic decay as well as $WZ$ production are presented. We compare to experimental measurements for the transverse momentum and the angular observable $\phi^*$.

Figures

Equations

$${\varphi }^{\ast }:=\tan \left(\frac{\pi -\Delta \varphi }{2}\right)\sin ({\theta }^{\ast }),\text{with}\cos ({\theta }^{\ast }):=\tanh \left(\frac{\Delta \eta }{2}\right).$$

$$q=q_1+q_2+\cdots +q_N.$$

$$q^{\mu }={\xi }_1{p}_1+{\xi }_2{p}_2+\mathcal{O}(q_{\perp }),$$

$$d\sigma =$$

$$\times \frac{1}{4\pi }{\int }_{-\infty }^{\infty }{d}^2{x}_{\perp }{e}^{-i{q}_{\perp }\cdot x_{\perp }}{\left(\frac{x_T^2{Q}^2}{b_0^2}\right)}^{-F_{ij}(x_{\perp },\mu )}{B}_i({\xi }_1,x_{\perp },\mu )B_j({\xi }_2,x_{\perp },\mu ),$$

$$d{\sigma }^{\mathrm{LO}}$$

$${\mathcal{H}}_{ij}({\hat{p}}_1,{\hat{p}}_2,q_1,...,q_N,\mu )=1+\frac{{\alpha }_s(\mu )}{4\pi }{\mathcal{H}}_{ij}^{(1)}({\hat{p}}_1,{\hat{p}}_2,q_1,...,q_N,\mu )+\mathcal{O}({\alpha }_s^2).$$

$${\mathcal{H}}_{q\bar{q}}^{(1)}=-2C_F{\ln }^2\frac{Q^2}{{\mu }^2}+6C_F\ln \frac{Q^2}{{\mu }^2}+h_0({\hat{p}}_1,{\hat{p}}_2,q_1,...,q_N).$$

$$h_0({\hat{p}}_1,{\hat{p}}_2,q_1,...,q_N)=2C_0({\hat{p}}_1,{\hat{p}}_2,q_1,...,q_N,{\mu }_{\mathrm{Mad}})+C_F\left[\frac{{\pi }^2}{3}+2{\ln }^2\frac{Q^2}{{\mu }_{\mathrm{Mad}}^2}-6\ln \frac{Q^2}{{\mu }_{\mathrm{Mad}}^2}\right].$$

$${\mathcal{H}}_{q\bar{q}}({\hat{p}}_1,{\hat{p}}_2,q_1,...,q_N,\mu )=U(Q^2,{\mu }_h,\mu ){\mathcal{H}}_{q\bar{q}}({\hat{p}}_1,{\hat{p}}_2,q_1,...,q_N,{\mu }_h),$$

$$L_{\perp }=\ln \frac{x_T^2{\mu }^2}{b_0^2},$$

$$B_i({\xi }_i,x_{\perp },\mu )=e^{h_i(L_{\perp },a_s)}{\bar{B}}_i({\xi }_i,x_{\perp },\mu ),$$

$$\frac{d}{d\ln \mu }{h}_i(L_{\perp },a_s)=C_i{\gamma }_{\mathrm{cusp}}{L}_{\perp }-2{\gamma }^i(a_s)$$

$$e^{g_i({\eta }_i,L_{\perp },a_s)}={\left(\frac{x_T^2{Q}^2}{b_0^2}\right)}^{-F_{ij}(L_{\perp },a_s)}{e}^{h_i(L_{\perp },a_s)}{e}^{h_j(L_{\perp },a_s)}.$$

$${\eta }_i\equiv {\eta }_i(Q^2,\mu )=\frac{C_i{\alpha }_s(\mu )}{\pi }\ln \frac{Q^2}{{\mu }^2}\sim 1,$$

$$q_{\ast }^2=Q^2\exp \left(-\frac{\pi }{C_i{\alpha }_s(q_{\ast })}\right).$$

$${\bar{B}}_i(\xi ,x_{\perp },\mu )=\sum _j{\int }_{\xi }^1\frac{dz}{z}{\bar{I}}_{i\leftarrow j}(z,x_{\perp },\mu ){\varphi }_j(\xi /z,\mu ),$$

$${\bar{I}}_{i\leftarrow j}(z)=\delta (1-z){\delta }_{ij}-a_s\left[{\mathcal{P}}_{i\leftarrow j}^{(1)}(z)\frac{L_{\perp }}{2}-{\mathcal{R}}_{i\leftarrow j}(z)\right]+\mathcal{O}(a_s^2).$$

$$\Delta {\bar{I}}_{i\leftarrow j}(z)=a_s^2\left({\mathcal{D}}_{i\leftarrow j}(z)-2{\beta }_0{\mathcal{P}}_{i\leftarrow j}^{(1)}(z)\right)\frac{L_{\perp }^2}{8},$$

$${\mathcal{D}}_{i\leftarrow j}(z)=\sum _k{\mathcal{D}}_{i\leftarrow k\leftarrow j}(z)=\sum _k{\int }_z^1\frac{du}{u}{\mathcal{P}}_{i\leftarrow k}^{(1)}(u){\mathcal{P}}_{k\leftarrow j}^{(1)}(z/u).$$

$${\bar{B}}_i(\xi ,x_{\perp },\mu )$$

$$\displaystyle+a_{s}^{2}\left({\cal D}_{i\leftarrow j}(z)-2\beta_{0}{\cal P}_{i\leftarrow j}^{(1)}(z)\right)\frac{L_{\perp}^{2}}{8}\bigg{]}\,\phi_{j}(\xi/z,\mu)$$

$$\equiv B_i^{(0)}(\xi ,\mu )+a_s{B}_i^{(1)}(\xi ,\mu )-a_s\frac{L_{\perp }}{2}{B}_i^{(2)}(\xi ,\mu )$$

$$+a_s^2{L}_{\perp }^2\left(-\frac{{\beta }_0}{4}{B}_i^{(2)}(\xi ,\mu )+\frac{1}{8}{B}_i^{(3)}(\xi ,\mu )\right),$$

$${\mathcal{M}}_n(Q,\mu ,q_T)=\frac{1}{4\pi }\int d^2{x}_{\perp }{e}^{-i{q}_{\perp }\cdot x_{\perp }}{e}^{g_F({\eta }_F,L_{\perp },a_s)}{L}_{\perp }^n.$$

$${\int }_{-\infty }^{\infty }{d}^2{x}_{\perp }{e}^{-i{q}_{\perp }\cdot x_{\perp }}=2\pi {\int }_0^{+\infty }d{x}_T{x}_T{J}_0(x_T{q}_T).$$

$${\mathcal{M}}_n(Q,\mu ,q_T)=-\frac{1}{\pi }\text{mbox}Im{\int }_0^{+\infty }d{x}_T{x}_T{K}_0(x_T{q}_T)e^{g_F({\eta }_F,L_{\perp },a_s)}{L}_{\perp }^n.$$

$${\mathcal{F}}_{ij}(Q$$

$$=$$

$$-\frac{a_s}{2}{\mathcal{M}}_1(Q,\mu ,q_T)\left[B_i^{(0)}({\xi }_1,\mu )B_j^{(2)}({\xi }_2,\mu )+B_i^{(2)}({\xi }_1,\mu )B_j^{(0)}({\xi }_2,\mu )\right]$$

$$\displaystyle+\frac{a_{s}^{2}}{4}{\cal M}_{2}(Q,\mu,q_{T})\Bigg{[}-\beta_{0}B^{(0)}_{i}(\xi_{1},\mu)B^{(2)}_{j}(\xi_{2},\mu)-\beta_{0}B^{(2)}_{i}(\xi_{1},\mu)B^{(0)}_{j}(\xi_{2},\mu)$$

$$\displaystyle\qquad\qquad+\frac{B^{(0)}_{i}(\xi_{1},\mu)B^{(3)}_{j}(\xi_{2},\mu)}{2}+\frac{B^{(3)}_{i}(\xi_{1},\mu)B^{(0)}_{j}(\xi_{2},\mu)}{2}+B^{(2)}_{i}(\xi_{1},\mu)B^{(2)}_{j}(\xi_{2},\mu)\Bigg{]}\,.$$

$$w={\left(\frac{{\alpha }_s(\mu )}{{\alpha }_s({\mu }_{\mathrm{Mad}})}\right)}^k\frac{{\mathcal{H}}_{ij}({\hat{p}}_1,{\hat{p}}_2,q_1,...,q_N,\mu ){\mathcal{F}}_{ij}({\xi }_1,{\xi }_2,q_T,\mu )}{{\varphi }_i({\xi }_1,{\mu }_{\mathrm{Mad}}){\varphi }_j({\xi }_2,{\mu }_{\mathrm{Mad}})}{w}^0,$$

$${\xi }_1{p}_1+{\xi }_2{p}_2=q.$$

Full text

Event-Based Transverse Momentum Resummation

Thomas Becher and Monika Hager

Albert Einstein Center for Fundamental Physics, Institut für Theoretische Physik,

Universität Bern, Sidlerstrasse 5, 3012 Bern, Switzerland

We have developed a framework for automated transverse momentum resummation for arbitrary electroweak final states based on reweighting tree-level events. It is fully differential in the kinematics of the electroweak final states, which facilitates a straightforward analysis of arbitrary observables in the small transverse momentum region. We have implemented the resummation at next-to-next-to-leading logarithmic accuracy and match to next-to-leading fixed-order results using the event generator MadGraph5\Q_aMC@NLO. Results for $Z$ and $W$ boson production with leptonic decay as well as $WZ$ production are presented. We compare to experimental measurements for the transverse momentum and the angular observable $\phi^{*}$.

1 Introduction

Since the emission of particles with large transverse momentum $q_{T}$ is suppressed by the strong coupling $\alpha_{s}(q_{T})$, most of the cross section at hadron colliders arises from events with low-$q_{T}$ radiation. This is true in particular for the production of electroweak bosons, i.e. $Z$’s, $W^{\pm}$’s and Higgs bosons. With its large data sets, the LHC can measure transverse momentum spectra of such bosons with exquisite precision and these results are used, for example, to determine the $W$ mass, as first achieved at the LHC in [1]. For this determination, the region of small $q_{T}$ is especially important.

The suppression of radiation by the coupling constant becomes ineffective at low transverse momentum, since it gets compensated by large logarithms of the ratio of the transverse momentum to the invariant mass of the electroweak boson. The all-order structure of these enhanced contributions was first understood by Collins, Soper and Sterman (CSS) [2], who showed that the cross section factorizes in (transverse) position space into a product of a hard function which encodes the virtual contributions to the electroweak boson production process and collinear functions describing the QCD emissions at low transverse momentum. The hard function depends on the electroweak process under consideration, while the collinear functions are universal and only distinguish quark-induced from gluon-induced processes. The result of CSS has been implemented by different authors and has also been rederived in the context of Soft-Collinear Effective Theory (SCET) [3, 4, 5] (see [6, 7, 8] for reviews). The result in SCET [9, 10] makes it clear that the process involves two distinct sources of large logarithms: i.) logarithms due to the different scales associated with the hard process and the radiation and ii.) logarithms due to the rapidity difference in the low-$q_{T}$ emissions from the partons flying along the beams to the left and right (these rapidity logarithms were not addressed in earlier work on transverse momentum resummation within SCET [11, 12, 13]). The first kind of logarithms are resummed by standard renormalization-group (RG) methods, while the second class is either exponentiated directly, using the collinear anomaly formalism [9], or resummed via a dedicated rapidity RG [14, 10]. While the factorization holds in position space, recently also methods to resum directly in momentum space have been developed [15, 16]. A number of computer codes for the resummation are available, e.g. RESBOS[17], CuTe [18], DYRes [19], MATRIX [20, 21] and RadISH [22]. For single-boson processes, the resummation is now performed up to next-to-next-to-next-to-logarithmic (NNNLL) accuracy [18, 22, 23, 24].

In the present paper, we present an efficient and flexible framework which achieves the following three goals:

It performs the resummation for arbitrary electroweak final states. 2. 2.

It computes any hadronically inclusive observable dominated by low-$q_{T}$ radiation and can take into account experimental cuts on the final-state leptons. 3. 3.

It automatically matches the resummed predictions to fixed-order results in kinematic regions were $q_{T}$ becomes large.

In order to achieve this flexibility, we make use of MadGraph5\Q_aMC@NLO [25], to compute the process-specific parts of the resummed cross section and supply it with the universal ingredients needed to achieve the resummation. Since we work at low $q_{T}$, we are close to Born level kinematics and can use the tree-level event generator to produce the leptonic final state. The automated one-loop code included in MadGraph5\Q_aMC@NLO is used to compute the virtual corrections to the hard scattering process. These results are then combined with the resummation factors and the universal collinear functions, which are tabulated and interpolated using PDF codes. More specifically, we start with tree level events which we reweight and boost to obtain resummed events. By analyzing these resummed events, we are able to impose cuts and extract arbitrary leptonic distributions. To perform the matching, we rely on the NLO fixed-order implementation of MadGraph5\Q_aMC@NLO.

The approach of reweighting fixed-order results to perform resummation was pioneered in [26], using MCFM [27]. Our implementation of the reweighting follows [28], in which we used MadGraph5\Q_aMC@NLO for automated jet-veto cross section resummation. Compared to this earlier work, transverse-momentum resummation involves a number of complications, which include the Fourier inversion back to momentum space and the necessity to account for recoil effects. The framework for transverse momentum resummation we use was developed in [29, 18]. This was implemented earlier in the CuTe code [18], which was however restricted to the inclusive $q_{T}$ spectrum of single bosons. Our new event-based framework extends it in the ways enumerated above and also introduces a novel efficient method to perform the matching and switch off the resummation at larger $q_{T}$. Our current implementation computes quark induced processes, has NNLL accuracy and matches at $\mathcal{O}(\alpha_{s})$ to fixed order. Extending it to higher accuracy requires two-loop ingredients which are not universally known, but could be implemented by hand for single-boson processes and for those diboson processes where they are available.

An important example of a kinematic quantity which correlates with the dilepton transverse momentum is the variable $\phi^{*}$ introduced in [30, 26, 31]. It is defined using the directions of the final-state leptons from the decay $Z\to\ell^{+}\ell^{-}$ as follows

[TABLE]

Here $\Delta\phi$ is the opening angle of the leptons in the azimuthal plane and $\Delta\eta=\eta^{-}-\eta^{+}$ the difference in their pseudorapidity. Since only angular measurements are needed to determine $\phi^{*}$, this quantity can be obtained even more precisely than $q_{T}$, which also requires lepton-energy measurements. Once $q_{T}$ approaches zero, the two leptons align back-to-back in the azimuthal plane and $\phi^{*}$ approaches zero. Indeed, computing the double differential cross section in $q_{T}$ and $\phi^{*}$, we observe a strong correlation among the two variables. We will compare to the $\phi^{*}$ measurements of ATLAS [32] in our paper, to illustrate our method in practice.

Our paper is organized as follows. In Section 2 we will review the factorization formula and the ingredients needed for NNLL accuracy. The implementation is discussed in Section 3, which provides details on the treatment of recoil effects, event generation, the matching to fixed order and the structure of the codes used to perform the resummation. In Section 4, we give numerical results for different processes, validate our results against CuTe and compare to experimental data. Conclusions and an outlook are presented in Section 5.

2 Factorization at low transverse momentum

We consider the scattering of protons with momenta $p_{1}$ and $p_{2}$ producing any number of massive electroweak bosons ($W^{\pm}$, $Z$, $H$) with momenta $q_{i}$, possibly decaying to leptons or photons, accompanied by hadronic radiation with total momentum $p_{X}$. The center-of-mass energy of the collision is $s=(p_{1}+p_{2})^{2}$ and the total electroweak momentum is

[TABLE]

We will analyze the cross section in the region where the transverse momentum $q_{\perp}$ is much smaller than the invariant mass $Q^{2}=q^{2}$ of the electroweak final state, and we use the notation $q_{T}=\sqrt{-q_{\perp}^{2}}$ to denote the positive scalar quantity associated with it. If we neglect the small transverse momentum, we can write the electroweak momentum as

[TABLE]

where the momenta $\hat{p}_{1}=\xi_{1}p_{1}$ and $\hat{p}_{2}=\xi_{2}p_{2}$ correspond to the large light-cone momentum components along the beam directions. Expanding the cross section around $q_{T}=0$ one obtains the factorization formula [9, 29, 18]

[TABLE]

which is equivalent to the CSS result [2]. We sum over partonic channels and integrate over the momentum fractions $\xi_{1}$ and $\xi_{2}$ of the partons entering the hard scattering process. We have introduced the abbreviation $b_{0}=2e^{-\gamma_{E}}$, where $\gamma_{E}$ is the Euler-Mascheroni constant. The formula, whose ingredients will be discussed below, involves a Fourier convolution over the transverse separation $x_{\perp}$ and holds up to terms suppressed by powers of $q_{T}^{2}/Q^{2}$. The cross section $d\sigma$ is inclusive in the hadronic radiation but completely differential in the electroweak momenta $q_{1},\dots,q_{N}$. To compute a specific cross section, such as the transverse momentum spectrum, one imposes suitable constraints on these momenta and integrates (4) over the electroweak phase space.

The structure of the cross section (4) is shown graphically in Figure 1 and is similar to the leading order cross section which reads

[TABLE]

Compared to the Born level result, the cross section (4) involves two additional ingredients. First of all, the resummed result involves a Fourier convolution with two beam functions ${B}_{i}(\xi,x_{\perp},\mu)$ instead of a convolution with ordinary parton distribution functions (PDFs) $\phi_{i}(\xi,\mu)$. The beam functions describe the soft and collinear QCD emissions which accompany the incoming parton, see Figures  1 and 2. We will discuss these functions and the associated Fourier integral over the transverse separation $x_{\perp}$ in more detail below. Let us note that for gluon-induced processes, such as Higgs production, two beam function structures arise. In this case the factorization formula involves a sum of two products of beam functions rather than just a product [33, 18]. However, the second structure first arises at NNNLL and is thus not relevant in the present paper.

Secondly, the resummed result also includes the virtual corrections to the Born level process. These are part of the hard function ${\cal H}_{ij}$, which is given by the loop contribution to the process, after subtracting its divergences in $\overline{\rm MS}$ renormalization. We write the expansion of the hard function in the form

[TABLE]

The one-loop hard function for quark-induced processes takes the form

[TABLE]

The $\mu$ dependence is universal since it is driven by the anomalous dimension of the operator with a single collinear quark field for each beam direction. All nontrivial information about the process resides in the scale independent piece $h_{0}$. For $Z$ boson production we have $h_{0}=C_{F}(-16+7\pi^{2}/3)$. For more complicated processes, we use MadGraph5\Q_aMC@NLO to compute the one-loop corrections, as described in detail in [28]. Specifically, running the code at an arbitrary reference scale $\mu_{\rm Mad}$, the hard function is related to the finite part $C_{0}$ of the virtual contribution obtained from MadGraph5\Q_aMC@NLO as follows:

[TABLE]

We observe that (7) suffers from large logarithms when $\mu^{2}\ll Q^{2}$, while the beam functions will involve large logarithms for $\mu^{2}\gg q_{T}^{2}$. To avoid this problem, we solve the RG equation of the hard function to evolve it to low values of $\mu$ at which the beam function is free of large logarithms. The result then takes the form

[TABLE]

and we choose the starting scale of the evolution to be $\mu_{h}\sim Q$. The analytical expression for the evolution factor $U(Q^{2},\mu_{h},\mu)$ is given in Appendix A.1.

Let us now discuss the Fourier integral. Despite the fact that it describes low-energy dynamics, the integral depends on the large scale $Q^{2}$ through the collinear anomaly [9]. This dependence exponentiates in (4) and is driven by the anomaly exponent $F_{ij}$, that was derived to two loops in [9] and has now even been determined at $\mathcal{O}(\alpha_{s}^{3})$ in [34, 35]. The beam functions ${B}_{i}$ are given by a convolution of a perturbative part, describing collinear and soft emissions at small transverse momentum, with the usual PDFs. The beam functions are illustrated in Figure 2 and will be discussed in detail below.

In perturbation theory, the functions ${B}_{i}$ are polynomials in the logarithm

[TABLE]

and it is useful to follow [29] and factor out their double logarithmic dependence by rewriting

[TABLE]

where we have introduced the abbreviation $a_{s}=\alpha_{s}(\mu)/4\pi$. The double-logarithmic exponent $h_{i}(L_{\perp},a_{s})$ is defined as the solution of the RG equation

[TABLE]

with boundary condition $h_{i}(0,a_{s})=0$. For quark-induced processes, we have $C_{i}=C_{F}$, while $C_{i}=C_{A}$ in the gluon case. The functions ${\bar{B}}_{i}$ are single logarithmic and it is convenient to combine the double logarithmic part with the anomaly into a single exponent

[TABLE]

While the beam functions are flavor dependent because they contain the PDFs, the exponents in this equation only depend on the color representation of the partons entering the hard scattering, i.e. they only differ between the quark case (fundamental representation, $g_{i}=g_{F}$) and the gluon induced process (adjoint representation, $g_{i}=g_{A}$). We list the exponent $g_{i}$ in the appendix in (A4), it was first given in [29]. The exponent depends on the variable

[TABLE]

which captures the large anomaly logarithms.

The Fourier integral in the factorization formula has some remarkable properties at very low transverse momentum. One would naively expect that the relevant scale for the integral tends to zero as $q_{T}\to 0$ but this is not the case, as was noted by Parisi and Petronzio [36] already before the all-order factorization of the cross section was fully understood. For very low $q_{T}$, the Fourier factor becomes ineffective and it is instead the Sudakov double logarithms inside the exponent $g_{i}$ which regularize the integration of the transverse separation. Analyzing the corresponding Gaussian integral, one finds that the associated scale $q_{*}$ is given by the value of $\mu$ at which $\eta_{i}$ becomes equal to one [29]

[TABLE]

For $Z$ production $q_{*}\approx 1.88\,{\rm GeV}$. In our numerical work, we therefore use $\mu=q_{T}+q_{*}$ as the default choice for the factorization scale. A consequence of the appearance of the dynamical scale $q_{*}$ is that the logarithm $L_{\perp}$, which usually counts as an $\mathcal{O}(1)$ quantity, must be counted as $L_{\perp}\sim\frac{1}{\sqrt{\alpha_{s}}}$ for $q_{T}\to 0$ [29]. One must therefore resum terms of the form $\alpha_{s}^{n}L_{\perp}^{2n}$, which now count as $\mathcal{O}(1)$. In (11), we have achieved this resummation by pulling out the factor $h_{i}$ from the beam functions and exponentiating it. The exponent $g_{i}(\eta_{i},L_{\perp},a_{s})$ given in (A4) contains all necessary terms to achieve $\mathcal{O}(\alpha_{s})$ accuracy also in the counting relevant for $q_{T}\to 0$. The appearance of the dynamical scale $q_{*}$ can also be understood from a momentum space perspective. Instead of soft radiation recoiling against the weak boson, the typical radiation for $q_{T}\to 0$ consists of QCD emissions at a scale $q_{*}$ recoiling against each other. This is the physical picture which underlies the momentum space formalism proposed in [16].

As depicted in Figure 2, the transverse-position dependent beam function ${\bar{B}}_{i}$ factorizes into a perturbative kernel $\bar{I}_{i\leftarrow j}$ describing the soft and collinear emissions at low transverse momentum, with the PDFs

[TABLE]

For NNLL resummation, we need the one-loop result for $\bar{I}_{i\leftarrow j}$ which takes the form [9, 29]

[TABLE]

The logarithmic piece is proportional to the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) splitting functions ${\cal P}_{i\leftarrow j}^{(1)}$ at one loop. For completeness, these are listed in Appendix A.3, together with the remainder functions ${\cal R}_{i\leftarrow j}$. As discussed above, for $q_{T}\to 0$, we must count $L_{\perp}\sim\frac{1}{\sqrt{\alpha_{s}}}$. To achieve uniform accuracy over the entire low $q_{T}$ region, we must also include the leading logarithmic piece of the two-loop beam functions

[TABLE]

where

[TABLE]

Sample contributions to ${\cal P}_{i\leftarrow j}^{(1)}$ and ${\cal D}_{i\leftarrow k\leftarrow j}$ for quark-induced processes are depicted in Figure 3. The right diagram in the figure shows that the quark flavor before and after radiation can differ in the two-loop terms. Explicit results for all relevant functions are listed in the appendix in (A8).

The complete beam function at NNLL accuracy is thus a second order polynomial in the logarithm $L_{\perp}$

[TABLE]

where the coefficients $B^{(m)}_{i}(\xi,\mu)$ are functions of the renormalization scale $\mu$ and $\xi$, the fraction of the incoming momentum which enters the hard process after the soft and collinear emissions. To be able to work efficiently with these functions, we tabulate them and use a PDF code for their interpolation.

With the coefficients $B^{(m)}_{i}(\xi,\mu)$ at hand, the Fourier integral reduces to a set of integrals involving the $n$-th power of a logarithm

[TABLE]

Since the logarithm $L_{\perp}$ only depends on $x_{T}^{2}=-x_{\perp}^{2}$, we can rewrite $q_{\perp}\cdot x_{\perp}=-x_{T}\,q_{T}\cos\phi$ and integrate over the azimuthal angle $\phi$, which yields

[TABLE]

Due to the oscillatory nature of the Bessel function $J_{0}(x_{T}q_{T})$, the numerical convergence is slow. It can be improved by using the identity $J_{0}(x_{T}q_{T})=\frac{2}{\pi}\mbox{Im}K_{0}(-ix_{T}q_{T})$ and then performing a Wick rotation $x_{T}\rightarrow ix_{T}$, which leads to

[TABLE]

We compute these integrals on the fly when running our code. Alternatively, one could implement the approximate analytical form developed in [37].

Expressed in terms of the integrals ${\cal M}_{n}$ and coefficients $B^{(m)}_{i}$, the final form of the Fourier integral, as implemented in our code, is

[TABLE]

Putting together the Fourier part ${\cal F}_{ij}$ with the RG-evolved hard function ${\cal H}_{ij}$ given in (9) we obtain the resummed cross section (4).

3 Event-based resummation

The basic method for the automated computation of the resummed cross section involves the following steps. We first generate events in the Les Houches Event File (LHEF) [38] format using the MadGraph tree-level event generator. We then use a script written during our earlier work on jet veto resummation [28] to compute the loop correction for each tree-level event and store this information in the event file. The event files are then processed using our code, which reads in the flavors $i$ and $j$ of the incoming partons, and their momentum fractions $\xi_{1}$ and $\xi_{2}$, as well as the $Q^{2}$ for each event. For a given value of $q_{T}$, the code then computes the function ${\cal F}_{ij}$ and constructs the RG-evolved hard function $\mathcal{H}_{ij}$ using the loop correction provided with each event and the RG evolution factor. The cross section at a given $q_{T}$ is then obtained as a weight factor

[TABLE]

where $w^{0}$ was the original weight of the tree-level event which was generated with the factorization and renormalization scales set equal to a reference value $\mu_{\rm Mad}$. The denominator is needed to remove the PDFs, which are replaced by the beam functions. The exponent $k$ is the power of $\alpha_{s}$ of the Born level process. For the quark-induced electroweak vector-boson processes we consider here $k=0$.

The procedure we just outlined is enough to produce a resummed transverse momentum spectrum for a given process, but we would also like to add experimental cuts on the momenta of the leptons produced in the decays of the vectors bosons and to compute related quantities such as $\phi^{*}$, which are constructed from lepton momenta. To be able to do so, we will construct a sample of events with different $q_{T}$, which we then analyze at the end. Before discussing how to generate these events, we must ensure that the hadronic recoil is transmitted to the electroweak final state.

3.1 Recoil effects

In the derivation of the factorization formula for the cross section at small transverse momentum $q_{T}$, one systematically expands in small momentum components. In particular, one drops the small transverse momentum of the partons entering the hard scattering process producing the electroweak bosons. In the factorization formula (4), the hard partons then have tree-level kinematics with

[TABLE]

Expanding away the small transverse momenta is appropriate for the computation of the QCD corrections associated with the large scale $Q^{2}$. It is also useful because the hard part of the process is then given by the tree-level amplitude $d\sigma^{0}_{ij}$ times corrections factors, allowing us to generate this part using a tree-level generator. Due to the expansion, momentum is no longer conserved exactly. We now have a mismatch between the electroweak kinematics, which has zero transverse momentum, and the hadronic part, in which the beam functions generate hadronic emissions at a low transverse momentum $q_{T}$. To correct for this, given a hadronic momentum $p_{X}^{\perp}=-q_{\perp}$ parametrized as

[TABLE]

we boost the entire tree-level event such that its total transverse momentum becomes $q_{\perp}^{\mu}$. The total cross section is invariant under this transformation, but the tree-level process now has two incoming partons with small transverse momenta. In our reweighting, we use the momentum fractions of the partons before the boost to determine the momentum fractions $\xi_{1}$ and $\xi_{2}$ for the beam functions. Doing so, we again neglect small momentum components but the advantage of proceeding this way is that the electroweak final state has the correct transverse momentum. This gives us access to the transverse-momentum distribution of individual final-state particles. The paper [19] has discussed different schemes for implementing recoil effects, which differ by terms suppressed by $q_{T}^{2}/Q^{2}$. These power suppressed terms are not captured by the resummation formula, but we will match to the fixed-order results to account for them up to $\mathcal{O}(\alpha_{s})$, see below.

3.2 Sampling of $q_{T}$ values

As indicated above, we want to generate a sample of events with different transverse momenta. The most natural way of doing this, would be to distribute the events according to the cross section, i.e. to compute

[TABLE]

Inverting this relation one obtains $q_{T}(z)$ and can then use a random number $z\in[0,1]$ to generate $q_{T}$ values. Proceeding in this way would yield events with equal weight, but a disadvantage of the procedure is that one would obtain only few events at larger $q_{T}$ values where the cross section is small.

In order to have a sample which also covers the region of larger $q_{T}$ values, we instead generate weighted events by sampling the $q_{T}$ values uniformly, i.e. we generate a random number $z$ and set

[TABLE]

Imposing a maximum $q_{T}$ value is necessary in any case because the resummed results for the cross section becomes unphysical at large values $q_{T}\gtrsim Q$. The value of $q_{\rm max}$ must be large enough to cover the entire region where the resummation is relevant. Choosing $q_{\rm max}\approx Q$ is clearly large enough. Using even larger values would not affect the final result once the matching to fixed order is performed (see Section 3.3 below), but would make the event generation inefficient. Writing $\Delta q\coloneqq dq_{T}/dz=q_{\rm max}$, the cross section integral takes the form

[TABLE]

In a MC evaluation of the above integral with $N$ events, each event thus contributes a weight

[TABLE]

Equivalently, we can assign a cross section

[TABLE]

to each event. In the practical implementation, we start with MadGraph tree-level events, generate $q_{T}$ according to (29) and a random angle $\phi\in(0,2\pi)$ to obtain the transverse momentum vector (27). Then we boost the event in the LHEF as discussed above and compute the event weight using (25). The boosted momenta and the event weight are then written back into the event file. In a final step, we analyze the resummed events, impose cuts and read out the observable of interest.

3.3 Matching to fixed order

Our resummed result captures logarithms which arise at small transverse momentum but expands away contributions which are suppressed by powers of $q_{T}^{2}/Q^{2}$. At larger transverse momentum these become more and more relevant and should be included. In order to obtain a result which covers all transverse momentum values, we combine our result with the fixed-order prediction. The labelling of fixed-order results is not uniform in the literature. We will use the term NLO to denote the $\mathcal{O}(\alpha_{s})$ result, so that the LO prediction is a $\delta$-function term at $q_{T}=0$. To avoid double-counting, the NLO-expanded NNLL-result must be subtracted from the sum,

[TABLE]

The first term on the right-hand side of (33) is our resummed result, the second term the fixed-order NLO result obtained from MadGraph5\Q_aMC@NLO, and the last term the resummed result expanded to NLO. The combination of the two latter terms is called the matching correction $\Delta\sigma$. The NLO-expansion of the resummed result can be obtained using the same reweighting method as for the resummed result; the relevant formula is given in Appendix A.4. The result of this naive matching procedure is shown in Figure 4.

While formally correct, the matched result (33) suffers from two problems. First of all, we do not recover the pure fixed-order result, even at very large $q_{T}$, because the resummed result includes higher-order terms in $\alpha_{s}$. Formally they are beyond the accuracy of the computation and can be kept, but since they are based on the $q_{T}\rightarrow 0$ limit, they can induce unphysical behavior at large $q_{T}$. Indeed, naively keeping those terms one ends up with a negative cross section at $q_{T}\gtrsim Q$. We should therefore switch off the resummation at large transverse momentum.

The second, more immediately visible problem concerns the other end of the spectrum. Both the fixed-order result and the NLO expansion of the resummed result diverge for $q_{T}\to 0$. The difference goes to zero, but numerically the cancellation is imperfect which leads to large numerical noise that renders the matched result useless for very small $q_{T}$. The numerical problems are especially visible because the resummed leading-power cross section is Sudakov suppressed for very small $q_{T}$. Of course the matching correction is not needed in this region and it can even be problematic to include it, because it contains (power suppressed) unresummed large logarithms. In the following, we will improve our matching scheme to solve both of these problems.

To eliminate the numerical noise at small $q_{T}$, we simply switch off the matching correction for very low $q_{T}<q_{0}$, where $q_{0}$ is a cutoff of the order of a few GeV. The cutoff $q_{0}$ is chosen large enough to avoid the numerical noise from the incomplete cancellation and small enough that the neglected matching correction, which parametrically scales as $q_{0}^{2}/Q^{2}$, is within the scale uncertainty of the resummed result. Both conditions are fulfilled for the choice $q_{0}=5\,{\rm GeV}$, which we adopt as our default value.

To switch off the resummation at large $q_{T}$, we introduce a transition function

[TABLE]

with $\lambda=\Delta\sigma/\sigma^{\text{matched}}$, where $\Delta\sigma$ is the matching correction and $\sigma^{\text{matched}}$ the naively matched cross section (33). We use $a=4$, $b=8$ and $C_{i}=C_{F}=4/3$ for the quark induced processes discussed here. The resulting functional form is plotted in Figure 5. The plot shows that we start switching off the resummation when the power suppressed terms amount to 20% of the result and switch if off completely once they are larger than 40% of the total cross section. While the value of $\lambda$ is affected by the numerical noise at low transverse momentum, this does not present a problem, since $t(\lambda)$ is equal to one in the region of low $\lambda$, see Figure 5.

With the cutoff and the transition function in place, formula (33) gets replaced by

[TABLE]

For low values of $q_{T}$, the function $t(\lambda)=1$ up to power corrections so that we reproduce expression (33) up to the fact that we switch the matching off at very small $q_{T}<q_{0}$. To obtain the matching correction $\Delta\sigma$, we evaluate the NLO result for the cross section in MadGraph5\Q_aMC@NLO for the observable under consideration, imposing $q_{T}>q_{0}$. One then subtracts from this the expanded resummed result imposing the same cutoff. For large values of $q_{T}$, we have $t(\lambda)\to 0$ so that the first term vanishes and we go back to the fixed-order result.

There are various other prescriptions to switch off resummation. One can eliminate the logarithms at large $q_{T}$ by suitably modifying their arguments with power suppressed contributions, or one can use special choices of the renormalization scales to achieve the same. An advantage of working with a transition function is that this approach is simple and transparent. In [19] the transition to fixed order was based on the value of $q_{T}$. Using instead the size of the power corrections as a measure is useful because it immediately generalizes to other observables such as $\phi^{*}$ or the lepton momentum distribution. In Figure 6 we plot the expectation value of $q_{T}^{2}/Q^{2}$, evaluated with the resummed cross section before matching, the size of the matching corrections and the value of the transition function for $q_{T}$, $\phi^{*}$ and $p_{T}^{\ell}$, the lepton momentum distribution. One observes that there is a good correlation between the quantity $\lambda$, which tracks the size of the power corrections, and the expectation value $\langle q_{T}^{2}/Q^{2}\rangle$ in the region of low transverse momentum. So the power corrections exhibit the expected scaling behavior.

Some care is required when computing expectation values using the resummed events since the resummed cross section becomes negative at large $q_{T}\gtrsim Q$ if the matching corrections are not included. We therefore restricted the sample to events with $q_{T}<Q$ when computing the expectation values for Figure 6. While the expectation value of $q_{T}^{2}/Q^{2}$ and $\lambda$ display similar behavior, we prefer to use $\lambda$ since it does not require any additional computations beyond the ingredients of (35). In Figure 7, we show the matched result based on the improved formula (35). One observes that the numerical noise at small $q_{T}$ is gone. The plot on the right shows the transition from the resummed result to the fixed-order case which takes place between $q_{T}$ values of $50-70\,{\rm GeV}$. Our results reduce to the fixed-order predictions for higher values of $q_{T}$.

3.4 Implementation

Let us briefly summarize the relevant steps to obtain a resummed cross section for $pp\to Y+X$, where $Y$ is the electroweak final state under consideration and $X$ the hadronic part of the final state which we assume to have low transverse momentum.

Use MadGraph5\Q_aMC@NLO to produce tree level events for the process $pp\to Y$, given as an event file in LHEF format. 2. 2.

Use the code virt_reweighter.py to compute the virtual correction for each tree-level event. Write this information back into the LHEF. 3. 3.

Run our code qT_reweighter.py. This code generates a transverse momentum for each event, boosts the electroweak particles to transmit the recoil and computes the resummed cross section at the given transverse momentum. To this end it assembles the hard function using the virtual correction computed in the previous step, computes the necessary Fourier integrals ${\cal M}_{i}$ and combines them with the interpolated beam function coefficients $B_{f}^{(k)}$. The code then writes the boosted vectors and the cross section as a weight back into the event file. In fact, to estimate the uncertainty, the result is computed not only with the default scale choices, but also after varying the scales $\mu$ and $\mu_{h}$ by a factor of two. Furthermore, we not only compute the resummed cross section, but also its fixed-order expansion to be able to perform the matching. The end result of this step is a statistical ensemble of events with transverse momentum containing different weights for different scale choices. 4. 4.

In the next step, we compute the cross section by analyzing the events for the given observable (such as $q_{T}$ or $\phi^{*}$) imposing also the relevant experimental cuts (such as the transverse momentum and rapidity cuts which ATLAS puts on the leptons). At this point, we fill a set of histograms for the observable under consideration, containing the resummed result as well as its NLO expansion for different scale choices. 5. 5.

Next we compute the NLO result using MadGraph5\Q_aMC@NLO and fill two NLO histograms. One is the full NLO result for the observable under consideration, the other one the NLO result with a cut $q_{T}>q_{0}$ needed to perform the matching using (35). 6. 6.

Finally, we combine the ingredients according to (35) to produce the final resummed and matched predictions, together with their scale variation bands.

The first three steps are observable independent and fully automated, while the analysis part in step 4 and the fixed-order computation in step 5 need to be set up for the specific observable under consideration. Our resummation codes can be obtained upon request and step-by-step instructions on how to use them can be found in [39].

4 Numerical results

We now present a few computations made with our code and compare to experimental predictions. Our predictions are based on version 2.6.4 of the MadGraph5\Q_aMC@NLO framework and unless stated otherwise, we adopt the default parameter values of this code. These include ${M_{Z}}=91.188\,{\rm GeV}$, $\alpha_{s}(M_{Z})=0.118$, $\alpha_{\rm EM}=1/132.507$, $G_{F}=1.16639\times 10^{-5}\,{\rm GeV}^{-2}$, and the derived quantities $M_{W}=80.419\,{\rm GeV}$ and $\theta_{W}=0.490912$. We will work with the MMHT 2014 NLO PDF set with $n_{f}=5$ flavors [40]. For the hard scale, we adopt the value $\mu_{h}=Q$, where the value of $Q$ is set dynamically, on an event-by-event basis. For the low scale, we choose $\mu=q_{T}+q_{*}$, where $q_{*}$ was defined in (15). The value of $q_{*}$ is obtained by numerically solving (15) for the value of $Q$ in the event. In our fixed-order computations and for the matching we set the renormalization and factorization scale to the hard scale, $\mu_{f}=\mu_{r}=\mu_{h}$. To estimate the uncertainties of our computation, we individually vary the scales $\mu$ and $\mu_{h}$ by a factor of two around their default values and take the envelope of the variations as our scale uncertainty. As expected, the scale bands are driven by the $\mu$ variation at low $q_{T}$. The $\mu_{h}$ variation becomes dominant at larger values, when we start to switch off the resummation.

The simplest process we can consider is $Z$ production with decay to leptons. We will compare our results to ATLAS measurements of the $q_{T}$ spectrum and the related observable $\phi^{*}$, which is obtained on the basis of angular measurements on the leptons. Before confronting experiment, it is interesting to compare NLL resummation to the predictions at NNLL accuracy. The corresponding spectra are shown in Figure 8. We observe that the scale uncertainties are reduced by about a factor two by going from NLL to NNLL. We also find that NNLL results lie within the NLL uncertainties. Predictions for the inclusive $q_{T}$ spectrum at the same accuracy and using the same resummation formalism were already obtained using the CuTe code [41] in [18]. We have verified that we reproduce these earlier results if we adopt the same value of the scales and compute to the same order in the improved expansion for $q_{T}\to 0$. Version 2 of CuTe includes resummation for the inclusive spectrum up to NNNLL and performs fixed-order matching up to NNLO.111Version 2.0.2 of CuTe incorporates the results for the four-loop cusp anomalous dimension [42] and the three-loop anomaly coefficient [34, 35]. The code thus achieves full NNNLL accuracy.

In Figure 9 we show the spectrum up to this accuracy. While the corrections are small at low $q_{T}$, the higher-order matching corrections at larger transverse momentum become significant, about 20%. One also observes that the fixed-order scale bands from varying $\mu=\mu_{f}=\mu_{r}$ underestimate the size of the corrections. The value of $\mu_{h}$ is mainly important for the normalization of the cross section. Choosing $\mu_{h}^{2}=M_{Z}^{2}$ gives a relatively low value for the NLL cross section, which then increases as one goes to higher orders. Adopting instead $\mu_{h}^{2}=-M_{Z}^{2}$ [43], the NLL result overshoots and the higher orders give negative corrections. Since the hard function arises as an overall factor, the choice of the hard scale plays only a minor role for the spectrum shown in Figure 9. In addition to scale variation, CuTe provides other options which can be used to estimate uncertainties. One can e.g. keep the higher-order corrections in the exponent, or expand them out. One can also vary the order of the improved expansion for $q_{T}\to 0$, which is implemented up to $[\alpha_{s}(q_{*})]^{5/2}$.

What is new compared to CuTe is that our code allows us to implement the ATLAS [32] cuts on the final-state leptons, which are restricted to have $p^{\ell}_{T}>20\,{\rm GeV}$ and pseudorapidity $|\eta|<2.4$. We focus on the $Z$ resonance region and restrict the invariant mass of the lepton system to the region $66\,{\rm GeV}<m_{\ell^{+}\ell^{-}}<116\,{\rm GeV}$. In contrast to our earlier work we are able to directly compare our resummed results to the measurement. Furthermore, we can also compute $\phi^{*}$ since we have access to the full lepton kinematics. The final ingredient for the comparison to ATLAS is the normalization, as our code produces cross sections not spectra. We first compute the fiducial cross section in the region $q_{T}=2-80\,{\rm GeV}$ from our matched result (using default scale choices) and then divide by this number to get the spectrum. The lower bound at $2\,{\rm GeV}$ is imposed to reduce sensitivity to possible non-perturbative effects. We also normalize the experimental result to the measured cross section in this momentum region. The upper bound was chosen because the unmatched resummed cross section turns negative at higher values of $q_{T}$ which would lead to unphysical behavior in the unmatched spectra shown in Figure 8.

In Figure 10 we plot our matched results for the $q_{T}$ and $\phi^{*}$ spectra, with the lepton cuts imposed by ATLAS [32]. The agreement is generally quite good, but at intermediate values we overshoot a little bit and our cross section is too small in the fixed order region at large $q_{T}$. Our fixed-order matching at $\mathcal{O}(\alpha_{s})$ only includes the leading term for $q_{T}\neq 0$ and thus has limited accuracy. The CuTe results shown in Figure 9 show that matching to $\mathcal{O}(\alpha_{s}^{2})$ would bring the cross section into agreement with the data. This is confirmed by [23] who match to the known $\mathcal{O}(\alpha_{s}^{3})$ result [44] and obtain a result which nicely agrees with the experimental data. In reference [23] the resummation is performed up to NNNLL, which leads to an excellent description of the data over the entire momentum range. In the context of the fixed-order computation, let us mention that in the matching scheme (35) with a cutoff $q_{0}$ on the matching corrections, we could extend the matching with some effort to $\mathcal{O}(\alpha_{s}^{2})$. To do so, one would use the MadGraph5\Q_aMC@NLO to perform a NLO computation of $Z+j$ with $p_{T}^{j}>q_{0}$ and also expand the resummed results one order higher in $\alpha_{s}$ to extract $\Delta\sigma$.

As discussed in the introduction, the variable $\phi^{*}$ was constructed as an alternative to $q_{T}$, as it can be measured more precisely. To illustrate their correlation, we show in Figure 11 a density plot of the cross section in $q_{T}$ and $\log_{10}\phi^{*}$. For a given $q_{T}$, there is a maximum possible value of $\phi^{*}$, which is obtained when the two leptons are produced at $\Delta\eta=0$. Determining the minimum $\Delta\phi$ and inserting it into the definition (1), one finds that $\phi^{*}_{\rm max}=q_{T}/Q$. The corresponding relation (for $Q=M_{Z}$) is shown as a dashed red line in Figure 11 and the red area above the line is kinematically excluded. The largest cross section is found near the maximum possible value of $\phi^{*}$ which demonstrates the close correlation among the two observables. In [45] it was observed that the logarithms in the $\phi^{*}$ distribution which arise at NLO can be obtained from the one in the transverse momentum spectrum by the substitution $q_{T}/M_{Z}\to 2\phi^{*}$, in agreement with our findings.

Since it is interesting in the context of the $W$-mass determination, we also show in Figure 12 the matched result for the lepton transverse momentum distribution in $Z$ production, imposing again the ATLAS [32] cuts as described above. Due to the lepton transverse momentum cut, this distribution starts at $p_{T}^{\ell}=20\,{\rm GeV}$. For this observable, resummation effects are especially important near the endpoint of the tree-level result at half of the mass of the produced boson. Indeed, for $p_{T}^{\ell}\approx M_{V}/2$ the distribution is dominated by low-$q_{T}$ events, while the matching becomes important at higher values of $p_{T}^{\ell}$, see Figure 6. The lepton momentum spectrum is much easier to measure than the transverse momentum of the weak boson, especially for the $W$ where one has to reconstruct the missing energy to obtain the boson momentum. To our knowledge no LHC measurements were presented for $p_{T}^{\ell}$, so that we cannot directly compare to data. The right-hand plot in Figure 12 shows the charged-lepton momentum distribution in $W^{+}$ production at $\sqrt{s}=7\,{\rm TeV}$ as in [1], imposing the same cuts on the charged lepton as in the $Z$ boson case.

As discussed earlier, predictions for $q_{T}$ and for $\phi^{*}$ at NNLL accuracy have been presented before and in the recent paper [23] even NNNLL results were given. The advantage of our implementation is its flexibility to perform transverse momentum resummation for arbitrary massive electroweak final states. As a first example, we have performed resummation for $W^{\pm}Z$ production. The transverse momentum spectrum of $Z$ bosons in $W^{\pm}Z$ production was measured in [46] and used to put limits on anomalous triple gauge boson couplings. In Figure 13 we show results for the diboson as well as the $Z$ boson transverse momentum spectrum. Resummation for diboson production has been studied earlier in the papers [47, 48, 49, 50, 21].

5 Conclusions

An important benefit of factorization — obtained using effective field theory or by other means — is universality: the same low-energy matrix elements typically arise in many different processes. The prime example in a collider context is provided by the parton distribution functions which capture the low-energy physics of arbitrary hard-scattering processes. The same universality is present for electroweak boson production at low transverse momentum. The accompanying QCD radiation is process independent and described by a Fourier convolution of two beam functions. In this paper, we have made use of this universality to automate transverse momentum resummation for arbitrary electroweak final states. To this end we have used MadGraph5\Q_aMC@NLO to generate the hard process together with its virtual corrections. The universal beam functions are tabulated similar to PDFs and then added to the hard process from the event generator using reweighting techniques. Our results automate the resummation to NNLL as well as the $\mathcal{O}(\alpha_{s})$ matching to fixed order and allow us to compute cross sections for arbitrary electroweak production processes with leptonic decays. Our event-based framework allows us to impose experimental cuts on the leptons which arise in the decay of the electroweak bosons and to have access to a variety of variables such as $\phi^{*}$.

There are a number of additions and improvements which could be incorporated into our framework in the future. First of all, it would be useful to also implement gluon-induced processes up to NNLL so that one could also study Higgs-boson production and associated processes. This would also allow one to study gluon-induced diboson production, which, although loop suppressed, can be sizeable due to the large gluon PDFs. Another improvement would be to extend the fixed-order matching at $q_{T}>0$ to $\mathcal{O}(\alpha_{s}^{2})$ by implementing the fixed-order expansion of the resummed result to one order higher. Finally, one could extend the resummation to NNNLL. However, this would require the two-loop hard functions, which would have to be supplied by hand for the cases where they are known. For diboson production, these functions were computed in [51, 52, 53] and are provided in numerical form through the VVamp code [54].

In the present paper, we have compared our results to the $q_{T}$ and $\phi^{*}$ distributions for $Z$ production measured by ATLAS and we have also computed sample diboson observables to demonstrate that our method also works for more complicated final states. In the future, it would be interesting to study diboson processes in more detail. Their sensitivity to new physics has been pointed out a long time ago [55, 56] and they are measured increasingly precisely at the LHC [57, 58, 59, 60, 61, 62, 63]. In this context our method would be quite useful. For example, it was proposed to impose a cut on $q_{T}$ in the process $pp\to WZ\to\ell\nu\bar{\ell}^{\prime}\ell^{\prime}$ to increase the sensitivity to new physics [64]. This specific cut is imposed to preserve the amplitude zero present in the Standard Model [65], but more generally these types of cuts are useful in order to prevent radiative corrections from washing out the operator structure one wants to probe with the measurements. Such a cut of course leads to Sudakov logarithms which can be resummed using our method. Since our code works with any model implemented into MadGraph5\Q_aMC@NLO, we can perform this resummation also for new-physics models or in Standard Model Effective Field Theory, which parametrizes the possible deviations in a model independent manner. The Feynman rules for this effective theory have been automated in a way suitable for the MadGraph5\Q_aMC@NLO framework [66, 67, 68]. We look forward to applications of our framework in the Standard Model and beyond!

Acknowledgments: This research is supported by the Swiss National Science Foundation (SNF) under grants 200020_182038 and 200020_165786. We thank Stephan Caspar, Lorena Rothen and Ding Yu Shao for discussions and Marcel Balsiger for beta testing our code and a thorough reading of the manuscript.

Appendix A Ingredients of the cross section at NNLL

For completeness we provide here all the ingredients for the resummed cross section as well as its fixed-order expansion needed to perform the matching.

A.1 Evolution of the hard function

The RG evolution factor $U(Q^{2},\mu_{h},\mu)$ in (9) needed to evaluate the hard function at a low scale has the form

[TABLE]

where $C_{q}=C_{F}$ for quark-induced processes and $C_{g}=C_{A}$ for the gluon case. The exponent of the evolution factor involves the single-logarithmic functions

[TABLE]

with $r=\alpha_{s}(\mu)/\alpha_{s}(\mu_{h})$, and $a_{\Gamma}(\mu_{h},\mu)$, which is given by the same expression after replacing the anomalous dimension by the cusp anomalous dimension, $\gamma^{i}_{n}\to\Gamma_{n}$. The anomalous dimension coefficients $\gamma^{q}_{n}$ and $\gamma^{g}_{n}$ for quark and gluon-induced processes, the cusp anomalous dimensions, as well as the $\beta$-function are listed in the appendix of [69]. The cusp anomalous dimension governs the Sudakov integral

[TABLE]

A.2 Anomaly exponent

In (13) in the main text we have combined the anomaly exponent with the double logarithmic part of the beam functions into the exponent [29]

[TABLE]

In addition to the various anomalous dimensions and the quantity $\eta_{i}$ defined in (14), the exponent involves the anomaly coefficient

[TABLE]

As discussed in the main text, the quantity $L_{\perp}$ has to be counted as $L_{\perp}\sim\frac{1}{\sqrt{\alpha_{s}}}$ for very small $q_{T}$, which explains the presence of higher-order terms multiplying powers of $L_{\perp}$ in the exponent (A4). These terms ensure that we include terms up to $\mathcal{O}(\alpha_{s})$ also in the modified power counting relevant for $q_{T}\to 0$.

A.3 Beam functions

We list here the ingredients of the perturbative part $\bar{I}_{q\leftarrow i}$ in (17) of the beam functions for quark induced processes. At one-loop level, these include the standard one-loop Altarelli-Parisi kernels

[TABLE]

which multiply the logarithm $L_{\perp}$ at one loop, as well as the remainder functions

[TABLE]

obtained in [9]. To correctly treat the region of very low transverse momentum, we further need the convolutions (18) of Altarelli-Parisi kernels, which multiply $L_{\perp}^{2}$ at two-loop order. The results for these quantities are [29]

[TABLE]

A.4 NLO expansion of the resummed cross section

To $\mathcal{O}(\alpha_{s})$, the expansion of the resummed cross section for $q_{T}>0$ for quark-induced processes is given by

[TABLE]

The above result arises from the expansion of ${\cal F}_{ij}$ which starts at $\mathcal{O}(\alpha_{s})$ for $q_{T}>0$. For the hard function we can thus use the leading order result ${\cal H}_{ij}=1$ at this accuracy. In our code, we implement the expanded result as a weight factor exactly as we did for the resummed result, see (25).