Generic dijet soft functions at two-loop order: correlated emissions

TL;DR

This paper introduces a systematic algorithm and a new software package for calculating two-loop soft functions in quantum chromodynamics, applicable to various collider observables, and confirms existing results while providing new predictions.

Contribution

The authors develop a universal parametrisation method and a software tool, SoftSERVE, enabling efficient two-loop calculations of soft functions for different collider observables.

Findings

Confirmed existing two-loop soft function results.

Provided new predictions for C-parameter and angularities.

Demonstrated applicability to both SCET-1 and SCET-2 soft functions.

Abstract

We present a systematic algorithm for the perturbative computation of soft functions that are defined in terms of two light-like Wilson lines. Our method is based on a universal parametrisation of the phase-space integrals, which we use to isolate the singularities in Laplace space. The observable-dependent integrations can then be performed numerically, and they are implemented in the new, publicly available package SoftSERVE that we use to derive all of our numerical results. Our algorithm applies to both SCET-1 and SCET-2 soft functions, and in the current version it can be used to compute two out of three NNLO colour structures associated with the so-called correlated-emission contribution. We confirm existing two-loop results for about a dozen $e^+e^-$ and hadron-collider soft functions, and we obtain new predictions for the C-parameter as well as thrust-axis and broadening-axis angularities.

Figures

Tables

Table 1. Table 2 : One- and two-emission measurement functions for the e + ​ e − superscript 𝑒 superscript 𝑒 e^{+}e^{-} event shape C-parameter, threshold and transverse-momentum resummation in Drell-Yan production.

Observable	$n$	$f(y_k,t_k)$	$F(a,b,y,t_k,t_l,t_{kl})$
C-parameter	$1$	$\frac{1}{1+y_k}$	$\frac{ab}{a\left(a+b\right)+\left(1+ab\right)y}+\frac{a}{a+b+a\left(1+ab\right)y}$
Threshold resum.	$-1$	$1+y_k$	$1+y$
$p_T$ resum.	0	$-2i\left(1-2t_k\right)$	$-2i\sqrt{\frac{a}{\left(1+ab\right)\left(a+b\right)}}\left(b\left(1-2t_k\right)+1-2t_l\right)$

Equations

$$S(\tau ,\mu )=\frac{1}{N_c}\sum _{i\in X}\mathcal{M}(\tau ;\{k_i\})\mathrm{Tr}|⟨X|T[S_n^{†}(0)S_{\bar{n}}(0)]|0⟩{|}^2,$$

$$\int d^{d}p{\left(\frac{\nu }{n\cdot p+\bar{n}\cdot p}\right)}^{\alpha }\delta (p^2)\theta (p^0),$$

$$d\sigma =H\cdot J_n\otimes J_{\bar{n}}\otimes S,$$

$$\displaystyle\mathcal{M}(\tau;\{k_{i}\})=\exp\big{(}-\tau\,\omega(\{k_{i}\})\,\big{)}\,,$$

$${\omega }^C(\{k_i\})=\sum _i\frac{k_i^{+}{k}_i^{-}}{k_i^{+}+k_i^{-}},$$

$${\omega }^{DY}(\{k_i\})=\sum _i(k_i^{+}+k_i^{-}),$$

$${\omega }^{p_T}(\{k_i\})=-2i\sum _i|k_i^{\perp }|\cos {\theta }_i.$$

$${\mathcal{M}}_0(\tau )=1$$

$$y_k=\frac{k_{+}}{k_{-}},k_T=\sqrt{k_{+}{k}_{-}},t_k=\frac{1-\cos {\theta }_k}{2},$$

$$\displaystyle\mathcal{M}_{1}(\tau;k)=\exp\big{(}-\tau\,k_{T}\,y_{k}^{n/2}\,f(y_{k},t_{k})\,\big{)}\,,$$

$$p_T=\sqrt{(k_{+}+l_{+})(k_{-}+l_{-})},y$$

$$t_k=\frac{1-\cos {\theta }_k}{2},t_l=\frac{1-\cos {\theta }_l}{2},t_{kl}=\frac{1-\cos {\theta }_{kl}}{2}.$$

$$\displaystyle\mathcal{M}_{2}^{corr}(\tau;k,l)=\exp\big{(}-\tau\,p_{T}\,y^{n/2}\,F(a,b,y,t_{k},t_{l},t_{kl})\,\big{)}\,,$$

$$F_A(a,b,y,t_k,t_l,t_{kl})$$

$$F_B(a,b,y,t_k,t_l,t_{kl})$$

$$p_T=\sqrt{k_T^2+\left(\sqrt{\frac{y_k}{y_l}}+\sqrt{\frac{y_l}{y_k}}\right)k_T{l}_T+l_T^2},y=\sqrt{y_k{y}_l}\frac{\sqrt{y_k}{k}_T+\sqrt{y_l}{l}_T}{\sqrt{y_l}{k}_T+\sqrt{y_k}{l}_T}.$$

$${\mathcal{M}}_2^{corr}(\tau ;k,l)\text{xlongrightarrow}{k}^{\mu }\to 0{\mathcal{M}}_1(\tau ;l)$$

$$F(a,0,y_l,t_k,t_l,t_{kl})=f(y_l,t_l).$$

$${\mathcal{M}}_2^{corr}(\tau ;k,l)\text{xlongrightarrow}{k}^{\mu }\parallel l^{\mu }{\mathcal{M}}_1(\tau ;k+l),$$

$$F(1,b,y_l,t_l,t_l,0)=f(y_l,t_l).$$

$$\omega (\{k\})=k_{+}^s{k}_{-}^{(1-s)}{f}_c(z_k,t_k),$$

$$\omega(\{k\})=\big{(}\sqrt{y_{k}}\,k_{T}\big{)}^{s}\,\Big{(}\frac{k_{T}}{\sqrt{y_{k}}}\Big{)}^{(1-s)}\,f_{c}(0,t_{k})=k_{T}\,y_{k}^{s-1/2}\,f_{c}(0,t_{k})\,,$$

$$n=2s-1=\frac{1}{p}-1.$$

$$S_n(x)=\mathbf{P}\text{exp}\left(i{g}_s{\int }_{-\infty }^{0}dsn\cdot A_s(x+sn)\right),$$

$$\int d^{d}p{\left(\frac{\nu }{n\cdot p+\bar{n}\cdot p}\right)}^{\alpha }\delta (p^2)\theta (p^0),$$

$$S_0(\tau ,\nu )=1$$

$$\displaystyle+\left(\frac{Z_{\alpha}\alpha_{s}}{4\pi}\right)^{2}(\mu^{2}\bar{\tau}^{2})^{2\epsilon}\;\bigg{\{}(\nu\bar{\tau})^{\alpha}\,S_{RV}(\epsilon,\alpha)+(\nu\bar{\tau})^{2\alpha}\,S_{RR}(\epsilon,\alpha)\bigg{\}}+\mathcal{O}(\alpha_{s}^{3})\,,$$

$$S_R(ϵ,\alpha )=\frac{(4\pi e^{{\gamma }_E}{\tau }^2{)}^{-ϵ}{\bar{\tau }}^{-\alpha }}{(2\pi {)}^{d-1}}\int d^{d}k\delta (k^2)\theta (k^0)\frac{|{\mathcal{A}}_R(k){|}^2}{(n\cdot k+\bar{n}\cdot k{)}^{\alpha }}{\mathcal{M}}_1(\tau ;k),$$

$$|{\mathcal{A}}_R(k){|}^2=\frac{64{\pi }^2{C}_F}{n\cdot k\bar{n}\cdot k},$$

$$k^{\mu }$$

$$\int d^{d}k\delta (k^2)\theta (k^0)$$

$$=\frac{{\pi }^{1/2-ϵ}}{\Gamma (1/2-ϵ)}{\int }_0^{\infty }d|{\vec{k}}_{\perp }|d{k}_{+}d{k}_{-}|{\vec{k}}_{\perp }{|}^{1-2ϵ}\delta (k_{+}{k}_{-}-|{\vec{k}}_{\perp }{|}^2){\int }_{-1}^{1}d\cos {\theta }_k{\sin }^{-1-2ϵ}{\theta }_k,$$

Full text

1 Introduction

Soft functions are an essential ingredient of QCD factorisation theorems. They describe the low-energy contribution to a scattering process, which is usually easier to compute than the analogous hard process in full QCD. Due to the eikonal form of the soft interactions, the soft functions can be represented by a vacuum matrix element of Wilson lines that point along the directions of the energetic, coloured particles in the scattering process.

As long as the underlying scale of the soft interactions is large enough, the soft functions can be calculated order-by-order in perturbation theory. At next-to-leading order (NLO) the calculation involves one-loop virtual and single real-emission contributions, which can be computed with standard techniques. Starting at NNLO and beyond, the singularity structure of the individual contributions becomes intricate and the divergences in the phase-space integrals overlap. The calculation of NNLO soft functions is needed for high-precision resummations and has attracted considerable attention in the past years Belitsky:1998tc ; Becher:2005pd ; Kelley:2011ng ; Monni:2011gb ; Hornig:2011iu ; Li:2011zp ; Kelley:2011aa ; Becher:2012za ; Ferroglia:2012uy ; Becher:2012qc ; vonManteuffel:2013vja ; Ferroglia:2013awa ; Czakon:2013hxa ; vonManteuffel:2014mva ; Boughezal:2015eha ; Echevarria:2015byo ; Luebbert:2016itl ; Gangal:2016kuo ; Li:2016tvb ; Campbell:2017hsw ; Wang:2018vgu ; Li:2018tsq ; Dulat:2018vuy ; Angeles-Martinez:2018mqh . Very recently, first results for N3LO soft functions have been presented Li:2016ctv ; Moult:2018jzp .

Whereas most of these calculations were performed analytically on a case-by-case basis, a systematic approach that exploits the universal structure of the soft functions is currently missing. The purpose of our work is to fill this gap, and as a first step we focus on soft functions that arise in processes with two massless, coloured, hard partons. The soft functions for these processes can be written in the form

[TABLE]

where $S_{n}$ and $S_{\bar{n}}$ are soft Wilson lines, and $n^{\mu}$ and $\bar{n}^{\mu}$ denote the directions of the hard partons with $n^{2}=\bar{n}^{2}=0$. For concreteness, we assume that the hard partons are in a back-to-back configuration ($n\cdot{\bar{n}}=2$), and that the Wilson lines are in the fundamental colour representation. The definition in (1) contains a trace over colour indices as well as a function $\mathcal{M}(\tau,\{k_{i}\})$, which specifies what is measured on the soft radiation $X$ with parton momenta $k_{i}$ for the observable under consideration. We will also see later that it is irrelevant whether $n^{\mu}$ and $\bar{n}^{\mu}$ are incoming or outgoing directions up to the order we consider, NNLO. Our method therefore equally applies to dijet observables in $e^{+}e^{-}$ annihilation, single-jet observables in deep-inelastic scattering, and zero-jet observables at hadron colliders. For convenience, we will refer to all of these cases with two massless, coloured, hard partons as dijet soft functions in the following.

The key observation of our analysis is that the soft matrix element in the definition (1) is universal, i.e. independent of the considered dijet observable. It is therefore possible to isolate the implicit divergences in the phase-space integrals with a universal parametrisation, and to compute the observable-dependent coefficients in an expansion in the dimensional regulator $\epsilon=(4-d)/2$ numerically. The dependence of the soft function on the observable is thus entirely confined to the measurement function $\mathcal{M}(\tau,\{k_{i}\})$, which acts as a weight factor for the numerical integrations. We will discuss the specific form we assume for the measurement function in the following section, where we will also learn that it is crucial to understand its properties in the singular limits of the matrix element.

The goal of our analysis thus consists in devising an algorithm that allows for an automated calculation of dijet soft functions to NNLO in the perturbative expansion. At NNLO the double real-emission contribution consists of three colour structures, which are often referred to as correlated ($C_{F}T_{F}n_{f}$, $C_{F}C_{A}$) and uncorrelated ($C_{F}^{2}$) emissions. As the phase-space parametrisations that are needed to factorise the divergences are different in both cases, we will concentrate in this work on correlated emissions, leaving the uncorrelated emissions for a future study BRT . For observables that obey the non-Abelian exponentiation (NAE) theorem Gatheral:1983cz ; Frenkel:1984pz , a dedicated calculation of the uncorrelated-emission contribution is in fact not needed, and we can therefore present complete NNLO results for a number of $e^{+}e^{-}$ and hadron-collider soft functions already in this work. This is, however, not true for observables that violate the NAE theorem, like jet-veto or grooming observables, which we will address in BRT (preliminary results can be found in Bell:2018jvf ).

As explained earlier, we aim at a numerical evaluation of bare dijet soft functions in an expansion in the dimensional regulator $\epsilon$. It is, however, well known that the phase-space integrals for certain soft functions suffer from rapidity divergences that are not regularised in dimensional regularisation (DR). This typically arises whenever the soft radiation is constrained to have small transverse momenta. Several prescriptions for the regularisation of the rapidity divergences have been proposed in the literature (see e.g. Chiu:2009yx ; Becher:2011dz ; Chiu:2012ir ; Echevarria:2015byo ; Li:2016axz ), and in this work we will use a variant of the analytic regulator introduced in Becher:2011dz . Specifically, this results in a modification of the generic $d$-dimensional phase-space measure of the form,

[TABLE]

where $\alpha$ is the rapidity regulator. The rapidity scale $\nu$ is introduced on dimensional grounds, similar to the renormalisation scale $\mu$ in conventional DR. The rapidity divergences then show up as poles in $1/\alpha$, and the renormalised soft function $S(\tau,\mu,\nu)$ depends on two scales $\mu$ and $\nu$. In the context of Soft-Collinear Effective Theory (SCET) Bauer:2000ew ; Bauer:2000yr ; Bauer:2001yt ; Beneke:2002ph , these soft functions are classified as SCET-2 observables, whereas those functions $S(\tau,\mu)$ that are not sensitive to the rapidity scale $\nu$ (and are well-defined in DR) refer to SCET-1 observables.

With the dimensional and the rapidity regulator in place, the bare soft functions can be evaluated in a double expansion in $\epsilon$ and $\alpha$. The main result of our analysis is an integral representation of a generic dijet soft function, in which all divergences are factorised. After introducing standard plus-distributions, the expansion in the various regulators can be performed and the coefficients of this expansion can be evaluated numerically. For the numerical integrations, we developed a new stand-alone program called SoftSERVE, which uses the Divonne integrator of the Cuba library Hahn:2004fe . The code contains a number of refinements to improve the convergence of the numerical integrations, which we will not discuss in detail in this work, but which are explained in the user manual of SoftSERVE. The SoftSERVE package is publicly available at https://softserve.hepforge.org/.

Although the main objective of our work is the computation of bare dijet soft functions, we go one step further and extract the ingredients that are needed in practical applications of resummation within SCET. To this end, we assume that the renormalised soft function obeys a multiplicative renormalisation group equation (RGE) in Laplace space, which allows us to define the (non-cusp) soft anomalous dimension $\gamma^{S}_{i}$ for SCET-1 soft functions and the collinear anomaly exponent $d_{i}$ for SCET-2 soft functions. In a previous study Bell:2018vaa , we derived integral representations for these quantities at the two-loop level for the same class of dijet soft functions we consider in the present work. Using SoftSERVE, which provides a script for the automated extraction of the resummation ingredients, it is thus possible to cross check the results of Bell:2018vaa , and to in addition obtain the finite (non-logarithmic) term $c^{S}_{i}$ of the renormalised soft function for SCET-1 observables and the bare soft remainder function $W^{S}_{i}$ for SCET-2 observables (both in Laplace space). According to the standard counting of logarithms for Sudakov problems (see Table 1), the two-loop expressions of $\gamma^{S}_{i}$ and $d_{i}$ are needed at next-to-next-to-leading logarithmic (NNLL) accuracy, while the two-loop constants $c^{S}_{i}$ and $W^{S}_{i}$ enter at the NNLL*′* level. SoftSERVE thus allows for increased logarithmic accuracy of SCET resummations, as was shown for the $e^{+}e^{-}$ event-shape angularities in Bell:2018gce , where the improvement was from NLL*′* to NNLL*′* (using preliminary results for the angularity soft function that were published in Bell:2015lsf ).

The outline of this paper is as follows: in Section 2 we define more precisely which dijet soft functions are amenable to our algorithm and we define the general properties as well as the specific form we assume for their measurement functions. In Section 3 we outline the technical aspects of the bare soft function calculation, and in Section 4 we specify the form we assume for the RGEs of both SCET-1 and SCET-2 soft functions. In Section 5 we examine several extensions of our formalism which are relevant, e.g., for multi-differential observables and soft functions that are defined in Fourier space. In Section 6 we briefly discuss the numerical implementation of our algorithm in SoftSERVE, and in Section 7 we present sample results for $e^{+}e^{-}$ and hadron-collider soft functions. All of our numerical results were generated using SoftSERVE, and the explicit examples we consider in Section 7 illustrate both the versatility and the usage of our code (template files for all soft functions considered in this work are provided in the SoftSERVE package). Most of these results are in fact already available in the literature at NNLO accuracy, and they hence provide strong cross-checks of our code, while also allowing us to study its numerical performance. Moreover, we obtain new predictions for the C-parameter, as well as thrust-axis and broadening-axis angularities. We finally conclude in Section 8 and provide some technical details of our analysis in the Appendix.

2 Measurement function

2.1 General considerations

We are concerned with soft functions that arise in processes with two massless, coloured, hard partons. A typical factorisation theorem of a dijet observable takes the form

[TABLE]

where the symbol $\otimes$ denotes a convolution in some kinematic variables, and $H$ is a hard function that contains the virtual corrections to the Born process at the large scale $Q$ of the scattering process. The jet functions $J_{n}$ and $J_{\bar{n}}$ encode the effects from collinear emissions into the directions $n^{\mu}$ and $\bar{n}^{\mu}$ of the hard partons, and the soft function $S$ describes the low-energetic interactions between the two jets (for incoming partons the collinear functions are often called beam functions). As the soft, long-wavelength partons cannot resolve the inner structure of the jets, the soft function only sees the directions of the hard partons as well as their colour charges. This is reflected in the definition (1) of the soft function, where the Wilson lines depend on the direction and the colour representation of the associated hard partons. For concreteness, we adopt the notion for $e^{+}e^{-}$ dijet observables in the following, and assume that the Wilson lines are given in the fundamental colour representation (the notation can be generalised to other processes by means of the colour-space formalism Catani:1996vz , and we will discuss some examples for hadron-initiated process below). We further assume that the hard partons are in a back-to-back configuration ($n\cdot\bar{n}=2$, along with $n^{2}=\bar{n}^{2}=0$), which is appropriate for both $e^{+}e^{-}$ and hadron-collider kinematics.

The soft function in the factorisation theorem (3) has a double-logarithmic evolution in the renormalisation scale $\mu$ and, possibly, also the rapidity scale $\nu$. In order to make the associated divergences explicit, we find it convenient to consider an integral transformation, which turns the convolution in (3) into a product. Apart from avoiding distribution-valued quantities, this considerably simplifies the solution of the associated RGEs. For many observables this is achieved by a Laplace transformation. Denoting the corresponding Laplace variable by $\tau$, we write the generic measurement function in the definition (1) of the soft function in the form

[TABLE]

where $\omega(\{k_{i}\})$ is a function of the final-state momenta $k_{i}$ that is specific to the observable. We thus assume that the distributions can be resolved by a single Laplace transformation, which implies that the soft function is differential in a single kinematic variable. We further assume that the Laplace variable has dimension 1/mass and that the measurement cannot distinguish between the two jets, i.e. the function $\omega(\{k_{i}\})$ is supposed to be symmetric under $n\leftrightarrow\bar{n}$ exchange. In addition, we allow for a non-trivial azimuthal dependence of the observable around the jet (or beam) axis. In other words, there may exist an external reference vector $v^{\mu}$ that singles out a direction in the plane transverse to the jet (or beam) direction. We finally impose two technical restrictions on the function $\omega(\{k_{i}\})$, namely its real part must be positive and it must be independent of the dimensional and rapidity regulators $\epsilon$ and $\alpha$.

Before we turn to some examples, let us recapitulate the assumptions that underlie our approach:

1. (A1)

Dijet factorisation theorem: We assume that the soft function is embedded in a factorisation theorem of the form (3), which refers to a process with two massless, colour-charged, hard partons. The hard partons can be in the initial or final state, and they are supposed to be in a back-to-back configuration ($n^{2}=\bar{n}^{2}=0$, $n\cdot\bar{n}=2$). The soft function for such dijet observables has a double-logarithmic evolution in the renormalisation scale $\mu$ and, possibly, also the rapidity scale $\nu$. 2. (A2)

Measurement function: We assume that the measurement function can be written in the form (4). Typically, this is achieved by taking a Laplace (or Fourier) transform of a momentum-space soft function, with $\tau$ being the associated Laplace (or Fourier) variable. In order to ensure that the phase-space integrals converge, we require that \Re\big{(}\omega(\{k_{i}\})\big{)}>0. More specifically, the function $\omega(\{k_{i}\})$ is allowed to vanish only for configurations with zero weight in the phase-space integrations, and it is furthermore assumed to be independent of the regulators $\epsilon$ and $\alpha$. 3. (A3)

Mass dimension: We assume that the variable $\tau$ has dimension 1/mass, and the function $\omega(\{k_{i}\})$, which only depends on the final-state momenta $k_{i}$, must therefore have the dimension of mass. This requirement could easily be relaxed to any positive mass dimension in the future, although we have not encountered any example that requires such a generalisation so far. 4. (A4)

$n$-$\bar{n}$ symmetry: We assume that the measurement cannot distinguish between the two jets, and the function $\omega(\{k_{i}\})$ is therefore symmetric under the exchange of $n^{\mu}$ and $\bar{n}^{\mu}$. This requirement could again easily be relaxed in the future, at the expense of doubling the number of input functions that need to be provided by the SoftSERVE user. 5. (A5)

Single-differential observables: We assume that the soft function only depends on one variable $\tau$ apart from the renormalisation and rapidity scales $\mu$ and $\nu$. Physically, this implies that the observable is differential in one kinematic variable. This requirement is in fact not strictly imposed in our approach (we will discuss multi-differential observables in Section 5), but we find it instructive to develop our formalism for this simplified class of observables first. 6. (A6)

Azimuthal dependence: Although we allow for a general azimuthal dependence of the observable around the jet/beam axis, we point out that the function $\omega(\{k_{i}\})$ is allowed to depend only on one angle $\theta_{i}$ per emitted particle in the $(d-2)$-dimensional transverse plane. This implies that the measurement is performed with respect to an external reference vector $v^{\mu}$, and the angle $\theta_{i}$ is then introduced as the angle between $\vec{v}_{\perp}$ and $\vec{k}_{i}^{\perp}$ in the plane transverse to the jet/beam direction.222For general $N$-jet soft functions with non-back-to-back kinematics, it was shown that two angles per emitted particle are required in the general case Bell:2018mkk . In addition, the function $\omega(\{k_{i}\})$ may depend on relative angles $\theta_{ij}$ between two emissions, which are defined as the angles between their respective transverse momenta $\vec{k}_{i}^{\perp}$ and $\vec{k}_{j}^{\perp}$.

Conditions (A1) and (A2) can be viewed as the strongest assumptions of our approach, although the generalisation to $N\geq 2$ jet directions with non-back-to-back kinematics is already in progress Bell:2018mkk (which also requires an extension of (A6)). Whereas the generalisation of (A5) will be discussed in Section 5, we already mentioned that assumptions (A3) and (A4) could easily be relaxed in the future. We further point out that our formalism is not limited to observables that obey the NAE theorem. For observables that violate NAE, the uncorrelated-emission contribution becomes non-trivial BRT ; Bell:2018jvf , but our method still allows for the calculation of the correlated-emission contribution, and hence it yields two out of three NNLO colour structures for NAE-violating observables.

Let us now see which type of observables fall into the considered class of soft functions. First, there are $e^{+}e^{-}$ event-shape variables that obey a hard-jet-soft factorisation theorem of the form (3) in the dijet limit. As an example we consider the C-parameter distribution, which was studied within SCET in Hoang:2014wka ; Hoang:2015hka . In an appropriate normalisation, its Laplace-space soft function can be written in the form (4) with

[TABLE]

where the plus- and minus-components represent the projections onto the $n^{\mu}$ and $\bar{n}^{\mu}$ directions with $k_{i}^{+}=n\cdot k_{i}$ and $k_{i}^{-}=\bar{n}\cdot k_{i}$. This function indeed has the dimension of mass, it is symmetric under $n\leftrightarrow\bar{n}$ exchange, and it does not depend on the regulators $\epsilon$ and $\alpha$. It is furthermore strictly positive, except for the trivial configuration with all $k_{i}^{\mu}=0$, which has zero weight in the phase-space integrations. The C-parameter is a single-differential observable with a trivial azimuthal dependence since the measurement is performed with respect to the jet axis itself.

As a second class of observables, we consider threshold resummation at hadron colliders. The classic example is Drell-Yan production, which was factorised in the form (3) using methods from SCET in Becher:2007ty (the collinear functions are the standard parton distribution functions in this case). In position space, the corresponding soft function can be written in the form (1) with a weight factor $\exp(-i\,x\cdot P_{X})$, where $P_{X}^{\mu}$ is the total momentum of the soft emissions. The vector $x^{\mu}$ thus plays the role of the reference vector $v^{\mu}$ in this case, and in the threshold kinematics it can be expanded as $x^{\mu}=(x^{0},\vec{0})$ in the centre-of-mass frame of the collision. As its spatial components vanish, the observable again has a trivial azimuthal dependence around the beam axis. In terms of $\tau=ix^{0}/2$, the position-space soft function can then be expressed in the form (4) with

[TABLE]

and one easily verifies that assumptions (A1)-(A6) are again satisfied for this observable.

We finally consider another class of hadron-collider soft functions, which arise in the context of transverse-momentum resummation. Taking again the Drell-Yan process as an example, the corresponding SCET analysis, which now involves beam functions that describe the effects from energetic initial-state radiation, can be found in Becher:2010tm . The respective position-space soft function can then be written in a similar form as the one for threshold resummation, except that the reference vector $x_{\perp}^{\mu}$ is now purely transverse to the beam direction. It therefore induces a non-trivial azimuthal dependence, which is precisely of the form we anticipated in (A6). Writing $\tau=|x_{\perp}|/2$, one finds that this soft function can again be written in the form (4) with

[TABLE]

With the usual exponential damping factor of a Fourier transform in mind, we can then argue that its real part is positive as required by assumption (A2). The function itself, however, now vanishes for non-trivial kinematic configurations, which single out a complicated hypersurface in the phase-space integrations. These configurations still have zero weight in the phase-space measure – as required by (A2) – but we will see later that our numerical results are less accurate for this observable in comparison with other examples that do not suffer from this problem. As the function $\omega^{p_{T}}(\{k_{i}\})$ is purely imaginary, we will also explain later in Section 5 that the numerical implementation of the transverse-momentum-dependent soft function in SoftSERVE requires special attention. One easily verifies that the remaining assumptions (A3)-(A5) are fulfilled for this observable.

The above examples should not be understood as an exhaustive list of observables that can be treated in our formalism; they should rather help to illustrate what kind of restrictions are imposed by assumptions (A1)-(A6). Other observables relevant e.g. for jet-veto resummation and jet-grooming observables also fall in the considered class of dijet soft functions. We will discuss further examples in Section 7.

2.2 Specific parametrisations

After these general considerations, we now present the specific form we assume for the measurement function in our calculation. At NNLO the measurement is performed on either zero, one, or two emitted partons.

According to (A3), the function $\omega(\{k_{i}\})$ is supposed to have the dimension of mass, and since it only depends on the final-state momenta $k_{i}$, it must vanish if there is no emission. We can therefore write the zero-emission measurement function in the form

[TABLE]

for all observables we consider.

For one emission, we have to find a phase-space parametrisation that allows us to control the implicit soft and collinear divergences in the phase-space integrations. We choose the variables

[TABLE]

where $y_{k}$ is a measure of the rapidity, $k_{T}$ is the magnitude of the transverse momentum, and $t_{k}$ parametrises the azimuthal dependence around the jet axis (with $\theta_{k}=\sphericalangle(\vec{v}_{\perp},\vec{k}_{\perp})$ as described in (A6)). The inverse of this transformation is then given by $k_{-}=k_{T}/\sqrt{y_{k}}$, $k_{+}=\sqrt{y_{k}}\,k_{T}$ and $\cos\theta_{k}=1-2t_{k}$.

In terms of these variables, we will see in the following section that the soft divergence arises in the limit $k_{T}\to 0$. The variable $k_{T}$ is in fact the only dimensionful quantity in this parametrisation, and according to (A3) the function $\omega(\{k\})=\omega(y_{k},k_{T},t_{k})$ must therefore be linear in $k_{T}$. The collinear divergences, on the other hand, emerge in the limits $y_{k}\to 0$ and $y_{k}\to\infty$. The $n$-$\bar{n}$ symmetry from (A4) then allows us to focus on one of the collinear limits, of which we choose the former. It turns out that the function $\omega(y_{k},k_{T},t_{k})$ may vanish or diverge as $y_{k}\to 0$, and that we have to control its scaling in this limit to properly extract the collinear divergence. Taken together, this motivates the following ansatz for the one-emission measurement function:

[TABLE]

where the power $n$ is fixed by the requirement that the function $f(y_{k},t_{k})$ is finite and non-zero in the limit $y_{k}\to 0$. For one emission, the observable is thus characterised by a parameter $n$ and a function $f(y_{k},t_{k})$ that encodes the angular and rapidity dependence.333It follows from (A2) that \Re\big{(}f(y_{k},t_{k})\big{)}>0 and that $f(y_{k},t_{k})$ is assumed to be independent of the regulators $\epsilon$ and $\alpha$.

Finding a suitable phase-space parametrisation for the double-emission contribution is much more involved. On the one hand the divergence structure of the matrix elements is more complicated, and on the other hand the measurement function must be controlled in various singular limits (whereas we only had to consider the limit $y_{k}\to 0$ in (10)). Moreover, we find that different parametrisations are needed for correlated and uncorrelated emissions. In this work we focus on the former, for which we introduce the variables

[TABLE]

The variables $p_{T}$ and $y$ are thus functions of the sum of the light-cone momenta, the quantity $a$ is a measure of the rapidity difference of the emitted partons, and $b$ is the ratio of their transverse momenta.444The variable $p_{T}$ should not be confused with the total transverse momentum of the emitted partons. In general the measurement function now depends on three angles since the emitted partons may not only see the reference vector $v^{\mu}$, but they will also see each other. The angles in (11b) are then introduced as $\theta_{k}=\sphericalangle(\vec{v}_{\perp},\vec{k}_{\perp})$, $\theta_{l}=\sphericalangle(\vec{v}_{\perp},\vec{l}_{\perp})$, and $\theta_{kl}=\sphericalangle(\vec{k}_{\perp},\vec{l}_{\perp})$, and the inverse transformation is now given by $k_{-}=a\,b\,p_{T}/(1+ab)/\sqrt{y}$, $k_{+}=b\,\sqrt{y}\,p_{T}/(a+b)$, $l_{-}=p_{T}/(1+ab)/\sqrt{y}$, $l_{+}=a\,\sqrt{y}\,p_{T}/(a+b)$, and $\cos\theta_{i}=1-2t_{i}$ for $i\in\{k,l,kl\}$.

After using the symmetries under $n\leftrightarrow\bar{n}$ and $k\leftrightarrow l$ exchange, we will see in the following section that the implicit phase-space divergences now arise in four limits:

• •

$p_{T}\rightarrow 0$, which corresponds to the situation in which both partons become soft;

• •

$y\rightarrow 0$, which reflects the fact that one of the partons becomes collinear to the jet direction $n^{\mu}$;

• •

$b\rightarrow 0$, which implies that the parton with momentum $k^{\mu}$ becomes soft;

• •

$a\to 1$ and $t_{kl}\to 0$, which means that the emitted partons become collinear to each other.

The first limit can in fact be treated in analogy to the limit $k_{T}\to 0$ in the one-emission case; since $p_{T}$ is the only dimensionful variable in the parametrisation (11), we know that the function $\omega(\{k,l\})=\omega(p_{T},a,b,y,t_{k},t_{l},t_{kl})$ must be linear in $p_{T}$. Yet, we still have to control the measurement function in the remaining three limits to make sure that we can properly extract the associated divergences.

What helps us in this situation is the underlying assumption in the factorisation theorem (3) that the observable is infrared safe. The one- and two-emission measurement functions are therefore not independent from each other, and we will derive explicit relations between them in the limit where one of the partons becomes soft ($b\to 0$) or both partons merge into a single parton ($a\to 1$ and $t_{kl}\to 0$) below. The last two limits from the above list are thus, as we say, protected by infrared safety, which means that we are guaranteed that the measurement function does not vanish in these limits since it must fall back to the one-emission case (which does not vanish for a generic configuration of the remaining parton). We therefore only have to consider the limit $y\rightarrow 0$ explicitly, which can be treated similarly to the limit $y_{k}\to 0$ in the one-emission case. Our ansatz for the correlated double-emission measurement function therefore reads

[TABLE]

where the function $F(a,b,y,t_{k},t_{l},t_{kl})$ is supposed to be finite and non-zero in the limit $y\to 0$. Interestingly, this is achieved by factorising the same power of the variable $y$ as in the one-emission case – see (10). We explain in Appendix A why this is so, and we address the physical meaning of the parameter $n$ in the next section.555We again demand that \Re\big{(}F(a,b,y,t_{k},t_{l},t_{kl})\big{)}>0 and that $F(a,b,y,t_{k},t_{l},t_{kl})$ is independent of any regulators as required by assumption (A2).

In order to extract the divergences of the bare soft function, we find it convenient to map the phase-space integrations onto a unit hypercube in the variables (9) and (11). While this can easily be achieved by exploiting the $n$-$\bar{n}$ symmetry for one emission, we will see in the following section that this leads to two independent regions in the two-emission case, which we label by the letters “A” and “B”. Our formulae therefore depend on two different versions of the two-emission measurement function, which are defined as

[TABLE]

Further explanations about the origin of these regions and the different representations of the measurement function in region B can be found in Section 3.3.

Before we come back to the explicit examples that we discussed towards the end of the last section, we derive the constraints from infrared safety that we mentioned earlier. To this end, we write the variables $p_{T}$ and $y$ from (11) in terms of $y_{k}$ and $k_{T}$ from (9) and the analogous variables $y_{l}$ and $l_{T}$ for the second emitted parton,

[TABLE]

In the limit in which the parton with momentum $k^{\mu}$ becomes soft, i.e. $k_{T}\to 0$, we thus see that $p_{T}\to l_{T}$ and $y\to y_{l}$. Infrared safety then tells us that the value of the observable should not change under infinitesimally soft emissions,

[TABLE]

which leads to the relation

[TABLE]

We can derive a similar relation between the one- and two-emission measurement functions in the limit in which the two partons with momenta $k^{\mu}$ and $l^{\mu}$ become collinear to each other. In this situation, which implies $y_{k}\to y_{l}$ and $t_{k}\to t_{l}$, we see that $p_{T}\to k_{T}+l_{T}$ and $y\to y_{l}$. As the value of the observable should again not change under collinear emissions,

[TABLE]

we obtain

[TABLE]

Relations (15) and (17) follow from the fundamental assumption that the observable that we factorised in (3) must be infrared safe. These relations are thus expected to hold for all observables we consider, and – as argued before – they guarantee that the measurement function does not vanish in two of the critical limits that we discussed above.

Starting from the observable definitions in (5) – (7), we can now easily derive the measurement functions for the C-parameter, threshold and transverse-momentum resummation in the phase-space parametrisations that we use in our calculation. The result is shown in Table 2, which illustrates that some observables have a non-trivial rapidity dependence, while others are sensitive to the reference vector $\vec{v}_{\perp}$ and therefore depend on the angular variables $t_{k}$ and $t_{l}$. From these expressions, we can verify that the functions $f(y_{k},t_{k})$ and $F(a,b,y,t_{k},t_{l},t_{kl})$ are finite in the limits $y_{k}\to 0$ and $y\to 0$, respectively, as this was the basis for extracting the corresponding values of the parameter $n$. We indeed see that these values can differ among the observables, and we will learn later that the case $n=0$ always corresponds to a SCET-2 observable. While we have so far introduced this parameter on purely technical grounds, we will see in the following section that it is related to the power counting of the momentum modes in the effective theory. Moreover, we can also easily verify that the constraints from infrared safety, (15) and (17), are satisfied for the considered class of observables.

2.3 Interpretation of the parameter $n$

We saw in the previous section that the parameter $n$ is related to the scaling of the observable in the soft-collinear limit, and we will indeed see later that it controls the double logarithmic contributions to the renormalised soft function. We also mentioned that the parameter $n$ allows us to distinguish between SCET-1 and SCET-2 observables, and we would like to understand why the values for the C-parameter ($n=1$) and threshold resummation ($n=-1$) are different, given that both observables are defined within SCET-1.

To this end, we go back to the factorisation theorem (3), which emerges in an effective field theory with hard, collinear, anti-collinear, and soft momentum modes. Denoting the small expansion parameter in the theory by $\lambda$, we associate the following power counting to the momenta $p^{\mu}=(\bar{n}\cdot p,n\cdot p,p_{\perp}^{\mu})$ with

• •

$p_{h}^{\mu}=Q(1,1,1)$        (hard)

• •

$p_{c}^{\mu}=Q(1,\lambda^{2p},\lambda^{p})$   (collinear)

• •

$p_{\bar{c}}^{\mu}=Q(\lambda^{2p},1,\lambda^{p})$   (anti-collinear)

• •

$p_{s}^{\mu}=Q(\lambda,\lambda,\lambda)$       (soft)

where $Q$ is the large scale in the process, and where we have allowed for a generic scaling of the collinear momenta that is controlled by a parameter $p>0$.

The factorisation theorem (3) then tells us that collinear, anti-collinear, and soft modes contribute to the observable at the same power, i.e. the observable $\omega(\{k_{i}\})$ must have the same scaling in $\lambda$ in the three regions.666One is often left with additive observables of the form $\omega(\{k_{i}\})=\omega_{c}\{k_{i}^{c}\})+\omega_{\bar{c}}(\{k_{i}^{\bar{c}}\})+\omega_{s}(\{k_{i}^{s}\})$, which leads to a multiplicative factorisation theorem in Laplace space. We know, however, that the observable scales as $\omega(\{k_{i}\})\sim\lambda$ in the soft region, since it has mass dimension one – see (A3) – and the power counting in the soft region is directly tied to the mass dimension. This can easily be verified for the examples in (5) – (7).

Our goal consists in establishing a relation between the parameter $n$ and the power-counting variable $p$, and for this purpose it is sufficient to focus on a single emission. In this case, the parameter $n$ controls the scaling of the observable in the collinear limit $y_{k}\to 0$ in the soft region. In the collinear region, on the other hand, we can exploit the hierarchy between the light-cone components, $k_{-}\gg k_{+}$, to express the observable in the form

[TABLE]

where the function $f_{c}(z_{k},t_{k})$ encodes an arbitrary dependence on the splitting variable $z_{k}=k_{-}/Q$ and the angular variable $t_{k}=(1-\cos\theta_{k})/2$, and we have furthermore used that the observable has mass dimension one. The parameter $s$ then varies among the observables, and we find $s=1$ for the C-parameter, $s=0$ for threshold and $s=1/2$ for transverse-momentum resummation.

We argued before that the observable must scale as $\lambda$ in the collinear region as well, and since $k_{-}\sim 1$ and $k_{+}\sim\lambda^{2p}$ for collinear momenta, we obtain $s=1/(2p)$. We can extract further information if we express (18) in terms of the variables from (9) and if we consider the soft limit $z_{k}\to 0$,

[TABLE]

which must match the expression in (10) in the collinear limit $y_{k}\to 0$. In particular, we see that the parameter $s$ controls the scaling of the observable with the rapidity-like variable $y_{k}$, which brings us to the desired relation,

[TABLE]

We thus see that the parameter $n$ is directly related to the power counting of the modes in the effective theory via $p_{c}^{\mu}=Q(1,\lambda^{2p},\lambda^{p})$, when the soft scaling is fixed to $p_{s}^{\mu}=Q(\lambda,\lambda,\lambda)$. For the C-parameter, for instance, the collinear modes scale as $p_{c}^{\mu}=Q(1,\lambda,\sqrt{\lambda})$ Hoang:2014wka , which implies that $p=1/2$ and hence $n=1$, which is in line with what we have found in Table 2. For transverse-momentum resummation, on the other hand, the relevant soft and collinear modes have the same virtuality, and so $p=1$ which translates into $n=0$ for a SCET-2 observable. However, the third example from our list appears to be peculiar, since $n=-1$ requires that $p\to\infty$, and the relevant collinear modes should therefore scale as $p_{c}^{\mu}=Q(1,\lambda^{2p},\lambda^{p})\to Q(1,0,0)$. We recall, though, that the collinear functions for threshold resummation are the standard parton distribution functions, and the power counting of the collinear modes is therefore not related to the threshold parameter $\lambda=1-M^{2}/\hat{s}$, but rather to the non-perturbative scale $\Lambda_{QCD}$ in this case. In other words, the relevant collinear modes scale as $p_{c}^{\mu}=Q(1,\varepsilon^{2},\varepsilon)$ with $\varepsilon=\Lambda_{QCD}/Q$ for this observable, and since $\varepsilon\ll\lambda$ this is indeed consistent with $p_{c}^{\mu}\to Q(1,0,0)$.

3 Calculation of the bare soft function

The definition (1) of what we call a generic dijet soft function depends on a measurement function $\mathcal{M}(\tau;\{k_{i}\})$, whose explicit form we discussed extensively in the previous section, as well as a matrix element of soft Wilson lines. The Wilson line associated e.g. with an incoming quark that travels in the $n^{\mu}$ direction is defined as

[TABLE]

where $\mathbf{P}$ is the path-ordering symbol, $A_{s}^{\mu}(x)=A_{s}^{\mu,A}(x)\,T^{A}$ is the soft gluon field and $T^{A}$ are the generators of SU(3) in the fundamental representation. The Wilson line associated with an outgoing quark has a similar representation, except that the integration now runs from [math] to $+\infty$. This subtle difference leads to the opposite sign in the $i\varepsilon$-prescription of the associated eikonal propagators, which – to the considered order in the perturbative expansion – is only relevant for the NNLO real-virtual interference. However, as we will see later, it turns out that the corresponding squared matrix element does not depend on the sign of this prescription, and our formulae therefore equally apply to soft functions with incoming and outgoing light-like directions. The Wilson lines associated with anti-quarks, moreover, are anti-path-ordered and their definition can be found e.g. in Chay:2004zn .

At leading order in the perturbative expansion, the soft matrix element in the definition (1) is trivial. Together with the form (8) of the zero-emission measurement function, this implies that the soft function is normalised to one at leading order. At higher orders, the bare soft function is subject to various divergences, which we control by a dimensional regulator $\epsilon=(4-d)/2$ and a rapidity regulator $\alpha$ that is needed only for SCET-2 observables. The latter is introduced on the level of the phase-space integrals as

[TABLE]

which is in the spirit of Becher:2011dz , except that our version respects the $n$-$\bar{n}$ symmetry that we assume on the level of the observable, see (A4). In this regularisation, the purely virtual corrections are scaleless and vanish at every order in perturbation theory. We are thus left with a single real-emission contribution at NLO, and with mixed real-virtual and double real-emission corrections at NNLO. The bare soft function can hence be written in the generic form

[TABLE]

where $\alpha_{s}$ is the renormalised strong coupling constant in the $\overline{\text{MS}}$ scheme, which is related to the bare coupling $\alpha_{s}^{0}$ via $Z_{\alpha}\alpha_{s}\,\mu^{2\epsilon}=e^{-\epsilon\gamma_{E}}(4\pi)^{\epsilon}\alpha_{s}^{0}$ with $Z_{\alpha}=1-\beta_{0}\alpha_{s}/(4\pi\epsilon)$ and $\beta_{0}=11/3\,C_{A}-4/3\,T_{F}n_{f}$. We furthermore introduced the rescaled variable $\bar{\tau}=\tau e^{\gamma_{E}}$ for convenience.

In the following we in turn address the computation of the single real-emission correction $S_{R}(\epsilon,\alpha)$, the mixed real-virtual interference $S_{RV}(\epsilon,\alpha)$, and the double real-emission contribution $S_{RR}(\epsilon,\alpha)$ for a generic dijet soft function.

3.1 Single real emission

In the normalisation (23), the single real-emission correction takes the form

[TABLE]

where $|\mathcal{A}_{R}(k)|^{2}$ denotes the corresponding soft matrix element and $\mathcal{M}_{1}(\tau;k)$ is the one-emission measurement function. At NLO the matrix element receives contributions from the four cut diagrams in Figure 1, where the double lines represent eikonal propagators associated with the $n^{\mu}$ and $\bar{n}^{\mu}$ Wilson lines. We find that the first two diagrams yield equal contributions, while the latter two vanish because they are proportional to $n^{2}=0$ or $\bar{n}^{2}=0$. The NLO squared matrix element is then given by

[TABLE]

where we suppressed the $i\varepsilon$-prescription of the eikonal propagators since it is irrelevant at this order.

We next decompose the gluon momentum in terms of light-cone coordinates

[TABLE]

with $k_{-}=\bar{n}\cdot k$, $k_{+}=n\cdot k$, and $n\cdot k_{\perp}=\bar{n}\cdot k_{\perp}=0$, along with $k^{2}_{\perp}=-\vec{k}_{\perp}^{2}$. As the one-emission measurement function only depends on one angle in the $(d-2)$-dimensional transverse space, see (A6), the phase-space measure can be simplified as

[TABLE]

where $\theta_{k}=\sphericalangle(\vec{v}_{\perp},\vec{k}_{\perp})$ is measured with respect to the external reference vector $v^{\mu}$. We switch to the parametrisation (9) and use the explicit form (10) of the one-emission measurement function to obtain

[TABLE]

with $\bar{t}_{k}=1-t_{k}$. As the $k_{T}$-dependence is universal among the considered class of observables, this integration can also be performed explicitly. We further use the $n$-$\bar{n}$ symmetry of the observable, which implies $y_{k}\to 1/y_{k}$ in the given parametrisation, to map the $y_{k}$-integration over the interval $[1,\infty]$ to an integral over $[0,1]$. We then arrive at the following master formula for the computation of the single real-emission correction:

[TABLE]

which is valid for arbitrary dijet soft functions that fall into the considered class of observables and which are characterised by the parameter $n$ and the function $f(y_{k},t_{k})$.

Upon expanding in the regulators $\alpha$ and $\epsilon$, the result exposes divergences, whose origin can be more clearly identified in (29). First, there is a soft singularity that arises in the limit $k_{T}\to 0$ and gives rise to the factor $\Gamma(-2\epsilon-\alpha)$. Second, the $y_{k}$-integral in (29) diverges in the collinear limits $y_{k}\to 0$ and $y_{k}\to\infty$, i.e. when the gluon is emitted into the $n^{\mu}$ or $\bar{n}^{\mu}$ directions. Due to the $n$-$\bar{n}$ symmetry, we can focus on one of these limits, and from (30) we finally read off that the collinear divergence is not regularised in dimensional regularisation for $n=0$, which is precisely the SCET-2 case we identified in Section 2.3.

For $n\neq 0$, on the other hand, the rapidity regulator $\alpha$ can be set to zero, and the expansion of the SCET-1 soft function starts with a $1/\epsilon^{2}$ pole, whose coefficient is controlled by the parameter $n$. This can be seen explicitly if we rewrite the divergent rapidity factor in terms of distributions according to

[TABLE]

As the function $f(y_{k},t_{k})$ is by construction finite and non-zero in the limit $y_{k}\to 0$, the remaining integrations in the expansion of (30) are well-defined and suited for a numerical integration. We in fact already presented the SCET-1 NLO master formula in Bell:2015lsf , and an earlier derivation along similar lines – although less general – was given in Hoang:2014wka . A similar NLO formula, valid also for non-back-to-back configurations, was presented in Kasemets:2015uus .

For SCET-2 soft functions with $n=0$, it is evident from (30) that the $y_{k}$-integration produces a $1/\alpha$ pole. It is in this case important that the $\alpha$-expansion is performed before the $\epsilon$-expansion, since the $\alpha$-regulator is supposed to regularise rapidity divergences only. The expansion of the factor $\Gamma(-2\epsilon-\alpha)$ therefore generates $1/\epsilon$ and $\alpha/\epsilon^{2}$ terms, which yield $1/(\alpha\epsilon)$ and $1/\epsilon^{2}$ poles on the level of the bare SCET-2 soft function, whose coefficients are related since they descend from the same Gamma function.

3.2 Real-virtual interference

The mixed real-virtual contribution is structurally identical to the single real-emission term. We now start from

[TABLE]

where the only difference is the soft matrix element $|\mathcal{A}_{RV}(k)|^{2}$, which can be calculated from the first diagram in Figure 2. Interestingly, the one-loop correction now depends on the $i\varepsilon$-prescription of the eikonal propagators on the amplitude level, but this dependence drops out in the interference with the Born diagram (see also Catani:2000pi ; Kang:2015moa ). One finds

[TABLE]

which again resembles the NLO matrix element (26), except that its expansion now starts with a $1/\epsilon^{2}$ pole, which is to be multiplied with the $1/(\alpha\epsilon)$ and $1/\epsilon^{2}$ poles of the subsequent phase-space integrations. The very fact that the matrix element (35) does not depend on the rapidity regulator $\alpha$ – which is implemented only on the level of the phase-space integrals in (34) – is a key advantage of the regularisation prescription from Becher:2011dz .

The subsequent calculation then follows along the same lines outlined in the previous section, and the master formula for the computation of the real-virtual interference becomes

[TABLE]

which we again already presented in the SCET-1 case in Bell:2015lsf .

3.3 Double real emissions

For the double real-emission contribution, we start from

[TABLE]

where $\mathcal{M}_{2}(\tau;k,l)$ is the two-emission measurement function. The respective soft matrix element now follows from the two-particle cut diagrams in Figure 2, which give rise to three colour structures – $C_{F}^{2}$, $C_{F}C_{A}$ and $C_{F}T_{F}n_{f}$ – of which the latter two are covered in this paper. The corresponding squared matrix elements are given by

[TABLE]

where we again suppressed the $i\varepsilon$-prescription of the propagators since it is irrelevant for the subsequent calculation. In comparison to (26) and (35), we observe that the singularity structure of the double real-emission contribution is much more complicated, and that it gives rise to overlapping divergences that are encoded e.g. in $(k_{-}+l_{-})$. The propagator $2k\cdot l=k_{-}l_{+}+k_{+}l_{-}-2|\vec{k}_{\perp}|\,|\vec{l}_{\perp}|\cos\theta_{kl}$, moreover, depends on the angle $\theta_{kl}$ between the two emissions in the transverse plane. Nevertheless, our basic strategy for the evaluation of the double real-emission correction is the same as before: We switch to the parametrisation (11), use the explicit form (12) of the two-emission measurement function, perform the observable-independent integrations and use symmetry arguments to map the integration domain onto a unit hypercube. However, the last two steps require us to find a suitable parametrisation for the angular integrations and to understand the implications of the $n$-$\bar{n}$ and $k$-$l$ symmetries, which we will address in the next two sections, before we present the master formula for the computation of the double real-emission contribution.

3.3.1 Angular parametrisation

According to assumption (A6), the two-emission measurement function in general depends on three angles, $\theta_{k}=\sphericalangle(\vec{v}_{\perp},\vec{k}_{\perp})$, $\theta_{l}=\sphericalangle(\vec{v}_{\perp},\vec{l}_{\perp})$, and $\theta_{kl}=\sphericalangle(\vec{k}_{\perp},\vec{l}_{\perp})$, and we would like to perform the integration over the remaining angles in the $(d-2)$-dimensional transverse space explicitly. This is similar in spirit to (28), where we retained the dependence on $\theta_{k}$, the only angle that arises in the one-emission measurement function.

To do so, we parametrise the vectors in the transverse space as

[TABLE]

which is illustrated in Figure 3 (for convenience we show unit vectors $\hat{k}_{\perp}$, $\hat{l}_{\perp}$, and $\hat{v}_{\perp}$ that point into the $\vec{k}_{\perp}$, $\vec{l}_{\perp}$, and $\vec{v}_{\perp}$ directions, respectively). In this parametrisation, the angular part of the phase-space measure becomes

[TABLE]

We thus have singled out a three-dimensional subspace spanned by the angles $\{\theta_{kl},\theta_{l},\theta_{5}\}$, in which the angle $\theta_{k}$ is given by

[TABLE]

This choice is of course arbitrary, but it is convenient since the matrix element depends on the angle $\theta_{kl}$ through the propagator $2k\cdot l$. In order to resolve the corresponding divergence, we want to keep the expression of said propagator simple, and $\theta_{kl}$ should therefore be one of the integration variables. The angles $\theta_{k}$ and $\theta_{l}$, on the other hand, only enter the calculation through the measurement function, which at the end of the day represents a weight factor for the numerical integrations. We therefore do not mind that the analytic expression of the measurement function becomes complicated once we express $\theta_{k}$ in terms of the integration variables through (41).

We next map the integration domain onto the unit hypercube by substituting as usual $\cos\theta_{i}=1-2t_{i}$ for $i\in\{k,l,kl,5\}$. In terms of $\bar{t}_{i}=1-t_{i}$, we then have

[TABLE]

and we arrive at

[TABLE]

which is almost the final expression, except that the $t_{5}$-integration suffers from spurious divergences that arise in the limits $t_{5}\to 0$ and $t_{5}\to 1$. These divergences are clearly unphysical, and they indeed cancel once they are combined with the prefactor $1/\Gamma(-\epsilon)$. They simply arise because we are resolving more angles than exist in four space-time dimensions.

Yet the $t_{5}$-divergences forbid a naive $\epsilon$-expansion on the integrand level, and we must therefore treat them in our formalism as if they were regular divergences. To this end, we first disentangle the two divergences by splitting the integration domain at $t_{5}=1/2$, and we subsequently rescale the two contributions as $t_{5}\rightarrow t_{5}^{\prime}/2$ and $t_{5}\rightarrow 1-t_{5}^{\prime}/2$, respectively. This yields for both cases

[TABLE]

where we now have to pay attention that we integrate over two copies of the actual integrand, one with the substitution $t_{5}\rightarrow t_{5}^{\prime}/2$ and the second one with $t_{5}\rightarrow 1-t_{5}^{\prime}/2$. But the integrand only depends implicitly on the variable $t_{5}$ through relation (42), and we therefore simply have to sum over two contributions in which the angle $\theta_{k}$, and hence the variable $t_{k}=(1-\cos\theta_{k})/2$, is resolved as777Notice that this definition of $t_{k}^{\pm}$ differs from the one we used in Bell:2018vaa .

[TABLE]

3.3.2 Symmetry considerations

We find it convenient to further map the entire integration domain onto the unit hypercube, and one can see in (44) that this has already been achieved for the angular integrations. We therefore only have to consider remappings that involve the remaining variables $\{p_{T},y,a,b\}$ in the parametrisation (11), which are a priori all defined on the interval $[0,\infty]$.

The idea is again similar in spirit to what we have seen in the single-emission case. There we arrived at the representation (29), in which the integration over the variables $k_{T}$ and $y_{k}$ both run from [math] to $\infty$. After performing the integration over $k_{T}$ analytically, we used the $n$-$\bar{n}$ symmetry to map the $y_{k}$-integration onto the unit interval. In the present case, the $p_{T}$-integration can similarly be performed analytically since the $p_{T}$-dependence is universal among the considered class of dijet soft functions – see (12). We then split the integrations over $y$, $a$ and $b$ at the value one, and substitute $y\to 1/y$, $a\to 1/a$, and $b\to 1/b$ to map the $[1,\infty]$ intervals onto $[0,1]$. Explicitly, this leads to eight different contributions

[TABLE]

where $\mathcal{I}(a,b,y)$ symbolically represents the integrand (after $p_{T}$-integration), which implicitly depends on the angular variables $t_{kl}$, $t_{l}$, and $t_{5}^{\prime}$ that we introduced in the previous section. Our goal thus consists in exploiting the symmetries under $n\leftrightarrow\bar{n}$ and $k\leftrightarrow l$ exchange to reduce the number of independent integrations.

We first consider the $n$-$\bar{n}$ symmetry, which is satisfied on the level of the observable because of (A4), and which is also respected by the form (22) that we use for the rapidity regulator. It is easy to see that under $n\leftrightarrow\bar{n}$ exchange

[TABLE]

Obviously, the measurement cannot distinguish between the two emitted partons, and the integrand is therefore also symmetric under $k\leftrightarrow l$ exchange, which implies

[TABLE]

In order to illustrate how we can make use of these symmetry considerations, let us for the moment focus on observables which do not depend on the angles $\theta_{k}$ and $\theta_{l}$. As the matrix element (38) does not depend on these angles either, the integrand in (46) is of the form $\mathcal{I}(a,b,y,t_{kl})$.888Recall that $\mathcal{I}(a,b,y)$ is a short-hand notation for $\mathcal{I}(a,b,y,t_{kl},t_{l},t_{5}^{\prime})$ in (46). We can then exploit the $n$-$\bar{n}$ and $k$-$l$ symmetries to reduce the integration to two regions with

[TABLE]

where the form of the second term is not unique, as we show now. This reduction is illustrated in Figure 4, where the effect of the symmetry transformations is shown for selected regions of the integration domain in figures (a) and (b). If plotted as eight stacked cubes in the three-dimensional $\{a,b,y\}$-space, the two symmetries ultimately enforce that the result of the integration in each of the four cubes marked in blue in figure (c) is the same. The eight cubes thus fall into two groups of four each, and the integration reduces to the form shown in (49), where the first term corresponds to the blue cube that is marked with dashes in figure (c). The second term, on the other hand, represents one of the three white cubes adjacent to this cube, and we see that it can be recovered by inverting one of the variables $a$, $b$, or $y$, each corresponding to one of the adjacent white cubes.

Up to this point we have assumed that the measurement function does not depend on the angles $\theta_{k}$ and $\theta_{l}$. In the general case, the preceding discussion still caries through, except that the $k$-$l$ symmetry now also exchanges the angles $\theta_{k}$ and $\theta_{l}$. We may therefore expect that we end up with four different regions in this case, since the symmetry transformation in figure (b) exchanges the role of the angular integrations. We are, however, always free to rename $\theta_{k}\leftrightarrow\theta_{l}$ in the transformed regions, which brings us back to (49). We should also point out that the symmetry considerations shall be exploited on the level of the full solid angle measure $d\Omega_{d-2}^{(k)}\,d\Omega_{d-2}^{(l)}$, i.e. before we choose which angle will be expressed through the integration variables (see the discussion in the previous section). In our setup, we essentially define the angle $\theta_{k}$ as the one we want to express in terms of the integration variables via (45).

In conclusion we find that the integration domain in the variables $a$, $b$, and $y$ can be mapped onto the unit hypercube, and in doing so we obtain two contributions which we denote by the letters “A” and “B”. Region A corresponds to the blue cube in Figure 4(c) that is highlighted by the dashed lines, and the corresponding integrand in (49) is just the original integrand $\mathcal{I}(a,b,y,t_{kl},t_{l},t_{5}^{\prime})$. Region B refers to any of the adjacent white cubes in Figure 4(c), which can also be mapped onto the unit hypercube by inverting either of the variables $a$, $b$, or $y$.

3.3.3 Master formula

We have now assembled all ingredients necessary to derive the master formula for the double real-emission correction. Starting from the representation (37), we first introduce light-cone coordinates and switch to the parametrisation (11). After inserting the explicit form (12) of the two-emission measurement function, we may then perform the integration over the dimensionful variable $p_{T}$ explicitly. We further use the symmetry arguments that we discussed in the previous section to map the integration domain onto the unit hypercube, and we resolve the angular phase-space measure as described in Section 3.3.1. We then obtain

[TABLE]

with $t_{k}^{\pm}$ from (45) and the colour factor is given by $C^{(n_{f})}=C_{F}T_{F}n_{f}$ or $C^{(C_{A})}=C_{F}C_{A}$. The integration kernels read

[TABLE]

As explained in the previous section, the master formula of the double real-emission correction consists of two contributions with measurement function

[TABLE]

in region A, whereas the one in region B is obtained by inverting either of the variables $a$, $b$, or $y$ with

[TABLE]

We stress that the three representations of this function need not be identical, since the symmetry arguments only guarantee that their integrals in (50) are equal, but not necessarily the integrands. One is therefore free to derive the measurement function in region B using any of the expressions on the right-hand side of (53).

From (50) we can analyse the divergence structure of the double real-emission correction. First, we find an explicit divergence that is encoded in the factor $\Gamma(-4\epsilon-2\alpha)$, which is associated with the limit $p_{T}\to 0$, i.e. the configuration where both emitted partons become soft. Second, we observe that the $y$-integral diverges in the limit $y\to 0$, which reflects the fact that one of the partons becomes collinear to the jet direction $n^{\mu}$. Similar to the single-emission case, this divergence yields a $1/\alpha$ pole for SCET-2 observables with $n=0$. Third, we identify an overlapping divergence in the limit in which the two emitted partons become collinear to each other, i.e. $a\to 1$ and $t_{kl}\to 0$. Fourth, the $C_{F}C_{A}$ contribution displays an additional divergence in the limit $b\to 0$, which implies that the parton with momentum $k^{\mu}$ becomes soft (due to the $k$-$l$ symmetry, the configuration with $l^{\mu}\to 0$ is mapped onto the same constraint). Finally, we observe that the expression diverges in the limit $t_{5}^{\prime}\to 0$, which is an unphysical divergence that is cancelled by the prefactor $1/\Gamma(-\epsilon)$.

For $n\neq 0$, the expansion of the $C_{F}T_{F}n_{f}$ structure thus starts with a $1/\epsilon^{3}$ divergence, whereas the $C_{F}C_{A}$ term has a $1/\epsilon^{4}$ pole because of the additional $b\to 0$ singularity. Except for the overlapping divergence, the singularities are in fact already factorised and can easily be isolated using an expansion in terms of plus-distributions in analogy to (31). As we explained in detail in Section 2.2, the functions $F_{A/B}(a,b,y,t_{k},t_{l},t_{kl})$ are by construction finite and non-zero in the singular limits999We note that the limit $t_{5}^{\prime}\to 0$ is indirectly protected by infrared safety, which can be seen as follows. The two-emission measurement function reduces to the one-emission function in the limit $b\to 0$ – see (15) – and it is therefore finite and independent of the angle $\theta_{5}$ in this limit. If, for finite values of $b$, the limit $t_{5}^{\prime}\to 0$ were to cause the observable to vanish or diverge, it would mean that the combined limit $t_{5}^{\prime}\to 0$ and $b\to 0$ was discontinuous. This is, however, not allowed because it would enable us to infer the presence of an infinitesimally soft second emission., and the remaining integrations are therefore well-defined upon an expansion in the dimensional regulator $\epsilon$. For SCET-2 soft functions with $n=0$, on the other hand, the $y$-integration generates a $1/\alpha$ pole which leads to additional $1/(\alpha\epsilon^{2})$ divergences for the $C_{F}T_{F}n_{f}$ and $1/(\alpha\epsilon^{3})$ divergences for the $C_{F}C_{A}$ colour structures.

In order to isolate the overlapping divergence, one could apply a standard sector decomposition strategy Binoth:2000ps , but we prefer to resolve it by means of an additional substitution,

[TABLE]

This change of variables matches a unit hypercube in the variables $(a,t_{kl})$ onto another unit hypercube in the new variables $(u,v)$, and it maps the overlapping divergence that arises in the limit $a\to 1$ and $t_{kl}\to 0$ onto the line $u\to 0$. The critical propagator then takes a particularly simple form

[TABLE]

and it turns out that the singularity from the limit $u\to 0$ is completely factorised, at the expense of increasing the complexity of the integrand. This is, however, a minor price to pay since the integrations are eventually performed numerically, and the above substitution does not worsen the convergence of the numerical integrations (see Section 6 for further refinements we implement to improve the numerical convergence).

4 Renormalisation

While the main objective of our work is the computation of bare dijet soft functions, we also extract the anomalous dimensions and matching corrections that are needed for resummations within SCET. This requires us to make some additional assumptions about the structure of the underlying RG equations. For many observables – including all examples that we discuss in Section 7 – the soft function renormalises multiplicatively in Laplace (or Fourier) space. We therefore focus on this particular class of observables in this section, leaving other soft functions that renormalise directly in momentum (or cumulant) space – like certain jet-veto observables – for a future study BRT .

4.1 SCET-1 observables

For SCET-1 observables with $n\neq 0$, we can set the additional regulator $\alpha=0$, and the expansion of the bare soft function takes the generic form

[TABLE]

where $x_{i}$ and $y_{i}$ are the NLO and NNLO coefficients at order $1/\epsilon^{i}$, respectively. The former are obtained by expanding the master formula (30) of the single real-emission contribution, while the latter are given by the sum of the real-virtual interference (36) and the double real-emission correction (50). It should be understood that the coefficients carry colour factors; these are given by $C_{F}$ in the case of the $x_{i}$, while the $y_{i}$ are sums of three different numbers multiplying the colour factors $C_{F}^{2}$, $C_{F}T_{F}n_{f}$, and $C_{F}C_{A}$. The correlated emission formulae provided in this paper yield the $C_{F}T_{F}n_{f}$ and $C_{F}C_{A}$ contributions, whereas the calculation of the $C_{F}^{2}$ correction is not covered in this work. For soft functions that obey the NAE theorem, this contribution is however proportional to the square of the one-loop correction and our results in Section 7 are therefore complete for this particular class of observables.

We now assume that the soft function renormalises multiplicatively in Laplace space, $S=Z_{S}S_{0}$, and that the renormalised soft function fulfils the RGE

[TABLE]

where $\Gamma_{\mathrm{cusp}}(\alpha_{s})$ is the cusp anomalous dimension and $\gamma^{S}(\alpha_{s})$ denotes the (non-cusp) soft anomalous dimension. The parameter $n\neq 0$ in the RGE is related to the power counting of the modes in the effective theory, as we explained in detail in Section 2.3. We find it convenient to define the non-cusp anomalous dimensions with a prefactor $2/n$, similar to the conventions we used in Bell:2018vaa . Expanding the anomalous dimensions as

[TABLE]

and using $Z_{\alpha}=1-\beta_{0}\alpha_{s}/(4\pi\epsilon)+\mathcal{O}(\alpha_{s}^{2})$, one can show that the RGE is solved to two-loop order by

[TABLE]

with $L=\ln(\mu\bar{\tau})$. The $Z$-factor $Z_{S}$ satisfies the same RGE (57), and its explicit solution to two-loop order is given by

[TABLE]

The universal expansion coefficients appearing in (59) and (60) read

[TABLE]

Equipped with this knowledge, we can extract the non-cusp soft anomalous dimension directly from the $1/\epsilon$ coefficients of the bare soft function using the relations

[TABLE]

while the non-logarithmic coefficients of the renormalised soft function follow from the finite terms via

[TABLE]

A strong check of our calculation is provided by the requirement that the higher poles $1/\epsilon^{j}$ with $j=2,3,4$ must vanish in the product of the $Z$-factor and the bare soft function.

4.2 SCET-2 observables

The SCET-2 case is slightly more complicated, owing to the double expansion in the regulators $\alpha$ and $\epsilon$. Starting from (23), the expansion of a bare SCET-2 soft function takes the generic form

[TABLE]

where $x^{i}_{j}$, $y^{i}_{j}$, and $z^{i}_{j}$ label the $1/(\alpha^{i}\epsilon^{j})$ coefficients of the single real-emission, double real-emission, and real-virtual interference term, respectively. We recall that the rapidity regulator is implemented on the level of the phase-space integrals, which explains the different powers of $(\nu\bar{\tau})$ in the NNLO correction. The coefficients again carry colour factors given by $C_{F}$ for the $x^{i}_{j}$, $C_{F}C_{A}$ for the $z^{i}_{j}$, and with $y^{i}_{j}$ containing three contributions proportional to $C_{F}^{2}$, $C_{F}T_{F}n_{f}$, and $C_{F}C_{A}$, of which the latter two are covered in this paper. The $C_{F}^{2}$ contribution is, on the other hand, again proportional to the square of the NLO correction for soft functions that obey the NAE theorem.

In the following we adopt the notation of the collinear anomaly approach Becher:2010tm ; Becher:2011pf to extract the relevant quantities for resummations in SCET-2. The formalism is equivalent to the rapidity renormalisation group (RRG) advocated in Chiu:2012ir , and we briefly comment on the translation into the RRG framework – including some subtleties about the choice of the rapidity regulator – at the end of this section.

In the collinear anomaly language, the bare soft function can be written in the form

[TABLE]

where we made the $\nu$-dependence explicit and we suppressed the terms divergent in $\alpha$, which cancel between the soft and collinear functions. The bare collinear anomaly exponent $\mathcal{F}_{0}(\tau)$ controls the logarithmic dependence on the rapidity scale $\nu$, and the bare soft remainder function $W_{0}^{S}(\tau)$ collects the terms that are not associated with the rapidity divergences. The latter is in fact meaningless in the collinear anomaly framework without knowledge about the corresponding collinear remainder function $W_{0}^{C}(\tau)$, since only their product obeys a well-defined RG equation in the $\overline{\text{MS}}$ scheme Becher:2011pf . As the collinear remainder function is not known in the chosen regularisation scheme for most of the observables we consider in Section 7, we disregard the soft remainder function $W_{0}^{S}(\tau)$ and focus on the collinear anomaly exponent $\mathcal{F}_{0}(\tau)$ in the following.

Following the procedure described in Becher:2012qc , we can extract the bare anomaly exponent from the bare SCET-2 soft function $S_{0}(\tau,\nu)$. This extraction is in fact subtle since the anomaly exponent is related to the coefficient of the logarithm $\ln(\nu^{2}\bar{\tau}^{2})$ – see (65) – rather than the associated $1/\alpha$ divergences Becher:2012qc . In terms of the expansion coefficients of the bare soft function from (64), we find that the bare anomaly exponent takes the form

[TABLE]

Owing to its place in the exponent, the anomaly coefficient renormalises additively in Laplace space, $\mathcal{F}_{0}=\mathcal{F}+Z_{\mathcal{F}}$, and the renormalised anomaly exponent satisfies the RGE

[TABLE]

which to two-loop order is solved by

[TABLE]

where again $L=\ln(\mu\bar{\tau})$ and the expansion coefficients of the cusp anomalous dimension and the beta function can be found in (61). The $Z$-factor $Z_{\mathcal{F}}$ satisfies a similar RGE as the anomaly coefficient, and its explicit form to two-loop order reads

[TABLE]

We can then extract the non-logarithmic terms of the renormalised anomaly coefficient (68) using the relations

[TABLE]

The cancellation of $1/\epsilon^{j}$ divergences with $j=1,2,3$ in the renormalised anomaly exponent then provides another strong check of our calculation.

We finally translate our findings into the RRG framework from Chiu:2012ir . Here the renormalisation is implemented directly on the level of the soft function rather than the anomaly exponent, $S=Z_{S}S_{0}$, and the Z-factor $Z_{S}$ absorbs both $1/\epsilon$ and $1/\alpha$ divergences according to a modified $\overline{\text{MS}}$ prescription. Furthermore, the renormalised soft function satisfies the RRG equation

[TABLE]

where

[TABLE]

which is solved by

[TABLE]

The solution can be compared to (65) in the collinear anomaly approach, bearing in mind that a similar relation holds among the renormalised quantities in this case. With the all-order solution to the RGE (67),

[TABLE]

we can then identify the $\nu$-anomalous dimension in the RRG approach with the collinear anomaly exponent,

[TABLE]

The comparison between (65) and (73) in addition allows us to express the renormalised soft remainder function as $W^{S}(\tau,\mu)=S(\tau,\mu,\nu_{s}=1/\bar{\tau})$. Interestingly, the latter has a well-defined $\mu$-evolution in the RRG framework that is governed by the RGE

[TABLE]

whereas the soft remainder function does not obey a simple RGE in the collinear anomaly approach (only the product of the soft and collinear remainder functions does so). Moreover, we find that the RGE (76) is not satisfied by our solution for SCET-2 soft functions, and the problem can be traced back to the way we have implemented the rapidity regulator. In other words, the RRG approach intrinsically makes specific assumptions about the form of the rapidity regulator, which in particular must be implemented on the level of connected webs Chiu:2012ir . We are not aware that this difference between the collinear anomaly and the RRG approach has been made so clearly in the literature before.

To summarise, for SCET-2 observables we determine the collinear anomaly exponent $\mathcal{F}(\tau,\mu)$ in (68) or, equivalently, the $\nu$-anomalous dimension $\gamma_{\nu}(\mu_{s})$ in (75) using the relations in (70). As our calculation yields the full bare soft function in (64), it also determines the bare soft remainder function $W_{0}^{S}(\tau)$, which is a useful input in the collinear anomaly approach if the corresponding collinear remainder function $W_{0}^{C}(\tau)$ is known in the same regularisation scheme. Our results for the soft remainder function are, on the other hand, not consistent with the RRG framework since we did not implement the rapidity regulator on the level of connected webs.

5 Revisiting our assumptions

Having established the theoretical framework for the calculation of the correlated-emission contribution to dijet soft functions, we now return to the list of assumptions that we outlined in Section 2. In particular, we can now better understand why these assumptions were made and how some of them could possibly be relaxed in the future. In this section, we in fact already introduce two extensions of our formalism that are valid for multi-differential and Fourier-space soft functions. We now address each of the assumptions from Section 2 in turn.

(A1) Beyond dijet factorisation

The soft functions we consider in this work are defined in terms of two light-like Wilson lines, and they are supposed to be embedded in a dijet factorisation theorem of the form (3). The soft function is, moreover, assumed to have a double-logarithmic evolution in the scales $\mu$ and, possibly, $\nu$. Soft functions that are blind to the jet directions only have a single-logarithmic evolution and thus they cannot be computed directly in the current formalism.101010In some cases it may be possible to write such soft functions as a difference of two dijet soft functions. This can clearly be seen from (10) and (12), where it is not possible to extract a value for the parameter $n$ if the observable is exactly zero in the collinear limit. As the parameter $n$ controls the double-logarithmic terms in the RGE (57), the current formalism cannot be applied to problems with single logarithms per loop order.

We further assumed that the two hard, massless partons are in a back-to-back configuration, i.e. $n\cdot\bar{n}=2$. Although the generalisation to arbitrary kinematical configurations is relatively straight-forward, we plan to relax this assumption only in the context of general $N$-jet soft functions Bell:2018mkk . It would also be interesting to extend the formalism to processes with massive hard partons ($n^{2}\neq 0$), which is relevant for top-quark related processes.

(A2) Loosening the constraints on the $\omega(\{k_{i}\})$ measure

Core to our approach is the structure of the function $\omega(\{k_{i}\})$ appearing in the exponent of the measurement function (4), and whose explicit one-emission and two-emission parametrisations were given in (10) and (12), respectively. While we required that \Re\big{(}\omega(\{k_{i}\})\big{)}>0, we already mentioned in (A2) that the observables are allowed to vanish for configurations with zero weight in the phase-space integrations. We now address this caveat more carefully, and we further elaborate on both the treatment of complex numbers and the regulator dependence of the measurement function.

Integrable divergences

According to the master formulae (30), (36), and (50), the functions $f(y_{k},t_{k})$ and $F_{A,B}(a,b,y,t_{k},t_{l},t_{kl})$ enter our formulae for the numerical integrations – after expansion in the various regulators – in terms of logarithms. It is therefore crucial that these functions do not vanish in the singular limits of the matrix elements, since otherwise the delta function associated with the divergence would put the argument of the logarithm to zero. The functions may, however, vanish for non-singular configurations with zero weight in the phase-space integrals, since the logarithms only constitute integrable divergences in this case. Of course, the logarithmic divergences may still pose a challenge for the numerical integrations, but for all examples we consider in Section 7 the integrable divergence seem to be under control. If, on the other hand, the functions $f(y_{k},t_{k})$ and $F_{A,B}(a,b,y,t_{k},t_{l},t_{kl})$ vanish over wide ranges of phase space (with non-zero weight), the logarithmic divergence is no longer integrable and these situations are therefore excluded by assumption (A2).

Fourier transforms and complex numbers

Whereas our formalism assumes that \Re\big{(}\omega(\{k_{i}\})\big{)}>0, the SoftSERVE implementation requires that \omega(\{k_{i}\})\big{)} is strictly real and non-negative. Whenever the soft function is defined in Fourier rather than Laplace space, the one-emission measurement function will appear as $[i\widetilde{f}(y_{k},t_{k})]^{2\epsilon+\alpha}$ in the NLO master formula (30) with a real-valued function $\widetilde{f}(y_{k},t_{k})$ (and similarly for the NNLO master formulae). In some cases like the one for threshold resummation in Drell-Yan production – see Section 2 – we can absorb the imaginary unit into the definition of the Laplace variable $\tau$, which brings the soft function into the standard form we assume in our framework. In other cases, such as the one for transverse-momentum resummation, this would leave us, however, with a real-valued function $\widetilde{f}(y_{k},t_{k})$ that can take on both positive and negative values, which is not allowed for the SoftSERVE implementation. We can nevertheless use our framework to calculate the real part of such soft functions by carefully treating the imaginary unit in the master formulae as a phase. A detailed description of this method is given in Appendix B, and as an example that exploits the Fourier-space extension we compute the soft function for transverse-momentum resummation in Section 7.

Regulator dependence

According to (A2), the function $\omega(\{k_{i}\})$ is assumed to be independent of the rapidity regulator $\alpha$ and the dimensional regulator $\epsilon$. While this is, of course, always fulfilled for a physical observable, this restricts the set of transformations that one may use to bring the soft function into the form (4). In particular, we are currently aware of a single observable that gives rise to an $\epsilon$-dependent $\omega(\{k_{i}\})$ measure: the $e^{+}e^{-}$ event shape jet broadening. Due to soft recoil effects, the broadening soft function requires a $(d-2)$-dimensional Fourier transform on top of the Laplace transformation to resolve all distributions Becher:2011pf , which brings the soft function into the form (4) with an $\epsilon$-dependent function $\omega(\{k_{i}\};\epsilon)$.

The main effect of the regulator dependence shows up in the expansion of the Laurent series in $\epsilon$, where different orders of the function $f(y_{k},t_{k};\epsilon)$ would contribute at different orders of the Laurent series. As the entire framework hinges on a precise understanding of the behaviour of the input functions in the singular limits, the discussion of the implications of infrared safety would have to be extended to the different orders in the function $f(y_{k},t_{k};\epsilon)$ itself. This could then lead to additional restrictions on the form of the regulator dependence. Such a discussion lies outside the scope of the present paper, and will only be revisited in the future if the specific need arises.

(A3) Mass dimension $\neq$ 1

Relaxing our assumption about the mass dimension of the $\omega(\{k_{i}\})$ measure is probably the easiest generalisation in our list. According to (10), the mass dimension determines the power of the variable $k_{T}$ in the measurement function, which is later integrated out analytically in (29) (similar arguments hold for the two-emission measurement function and the variable $p_{T}$). However, the latter integration is perfectly convergent and well-behaved for any positive, non-zero mass dimension, and the only changes manifest in different numerical constants in some places, like the argument of the Gamma functions or certain exponents.

Switching to different mass dimensions would therefore only result in slightly modified master formulae in our approach. While this could easily be implemented, we have not yet encountered any observable which would require such a tune.

(A4) Broken $n$-$\bar{n}$ symmetry

In its current form we assume that the $\omega(\{k_{i}\})$ measure is symmetric under $n\leftrightarrow\bar{n}$ exchange, which is not necessarily the case for all observables. In Section 7 we consider e.g. the hemisphere soft function, which depends on two invariant masses $M_{L}$ and $M_{R}$ that are not invariant, but rather mapped onto each other, under $n\leftrightarrow\bar{n}$ exchange.

Relaxing the $n$-$\bar{n}$ symmetry is again possible at the expense of doubling the number of input functions that need to be provided by the SoftSERVE user. The required input currently includes the measurement functions $F_{A}$ and $F_{B}$ and the parameter $n$ (the function $f$ is internally determined using the infrared-safety constraints). If the $n$-$\bar{n}$ symmetry is given up, four functions $F_{A}$ through $F_{D}$ and two parameters $n_{(A,B)}$ and $n_{(C,D)}$ would be required instead. The latter reflects the fact that the rapidity scaling of the observable may differ between the two light-like directions if the $n$-$\bar{n}$ symmetry is broken. The extension from two to four input functions can also easily be seen from the symmetry considerations in Section 3.3.

In practice, a broken $n$-$\bar{n}$ symmetry can already be emulated with the current version of SoftSERVE by averaging separate results for the $\{F_{A},F_{B}\}$ and $\{F_{C},F_{D}\}$ sets, similar to a procedure we use for the hemisphere soft function in Section 7 (with special attention required to get the angular dependence correct). The substitutions to generate the relevant input functions are then $y\to 1/y$ to derive $F_{C}$ from $F_{A}$, and $a\to 1/a$ or $b\to 1/b$ to derive both $F_{B}$ from $F_{A}$, and $F_{D}$ from $F_{C}$111111In Figure 4 regions C and D build up the upper layer of four cubes..

(A5) Multi-differential observables

Soft functions for exclusive observables typically depend on more than one kinematic variable, requiring multiple Laplace transformations to resolve all distributions. We can in such cases choose the first Laplace variable $\tau_{1}$ to have dimension 1/mass, and keep the remaining variables $\tau_{i}$ for $i\geq 2$ dimensionless. This generalises our ansatz (4) for the measurement function to

[TABLE]

from which one can derive the one- and two-emission measurement functions $f(y_{k},t_{k};\tau_{2},\ldots)$ and $F(a,b,y,t_{k},t_{l},t_{kl};\tau_{2},\ldots)$ via the usual procedure. In essence, multi-differential soft function can thus be computed by treating the Laplace variables $\tau_{i}$ for $i\geq 2$ as parameters, which need to be sampled over. We demonstrate this strategy in Section 7 by calculating two double-differential soft functions: the hemisphere soft function in $e^{+}e^{-}$ collisions and the soft function for exclusive Drell-Yan production.

(A6) Extended angular dependence

Finally, we note that the angular parametrisation of our integrals, which we presented in detail in Section 3.3.1, is made under the assumption of a back-to-back kinematic setup for the soft Wilson lines. For dijet observables where this is not the case, the angular parametrisation is not sufficient, since more dynamic angles can be resolved by the measurement function. Our master formulae therefore need to be revisited for non-back-to-back dijet observables, and we in fact already implemented such a generalisation in the $N$-jet extension of our formalism Bell:2018mkk .

6 Numerical implementation

The master formulae we derived in the preceding sections are in principle complete as they render the sources of all singularities manifest and allow for a numerical evaluation across a wide field of observables with only a few required properties. In fact, the master formulae can already be used to derive semi-analytic expressions for the anomalous dimensions and collinear anomaly exponents (Bell:2018vaa, ), although the matching corrections seem to be a bit out of reach in such an approach due to their complexity. Solving the equations analytically is nevertheless possible in some isolated cases, in particular in the absence of a non-trivial angular dependence (see the C-parameter in Section 7).

Still, the master formulae are not yet ideally suited for a numerical implementation because of the presence of an overlapping divergence in the variables $a$ and $t_{kl}$. The overlapping divergence could in principle be resolved by multiple sector decomposition steps (it produces three sectors due to the square in $(1-a)^{2}$) before the subtraction and the numerical evaluation. This is precisely what the programs SecDec (Carter:2010hi, ; Borowka:2012yc, ; Borowka:2015mxa, ) and its successor pySecDec (Borowka:2017idc, ) were designed to do, and we can indeed evaluate the master formulae with these programs121212We in fact extensively used (py)SecDec to cross-check our SoftSERVE numbers.. However, there are a few ways to improve on (py)SecDec as regards our purposes, chiefly because there are simplifications possible that (py)SecDec — as a program designed with a larger scope of applications in mind — cannot easily exploit.

The core insight that motivates our tailored numerical approach is the substitution (54), which removes the overlapping divergence in favour of a monomial divergence $u^{-1-2\epsilon}$, at the cost of increasing the complexity of the integrand. Disentangling the divergence between $a$ and $t_{kl}$ means that all divergences are now present in monomial form, which makes a subtraction and expansion procedure trivial. It should be noted that using this substitution can also speed up (py)SecDec runs, as the sector decomposition steps are of course no longer needed in this case either anymore.

The program we subsequently wrote is called SoftSERVE, and it is publicly available at https://softserve.hepforge.org. SoftSERVE mainly implements the master formulae from this paper in C++ syntax, and integrates them using the Cuba library Hahn:2004fe . It is therefore subject to the same assumptions and capabilities of the formalism we developed in this paper, with one additional constraint: It is limited to strictly real measurement functions131313The case of Fourier-space soft functions is special, and will be revisited at the end of this section and in Appendix B.. Below, we will lay out the main reasons for forgoing (py)SecDec and writing a dedicated program, while delegating the technical details of the C++ implementation and the ultimate structure of the program to the SoftSERVE manual, which is provided alongside the program.

The single most important motivation for writing a dedicated program is the problem of rounding errors in conjunction with plus-distributions. The rather technical details of this problem can be found in the SoftSERVE manual although, put succinctly, problems can arise if plus-distributions are integrated against functions that involve large cancellations in the limit originating from the plus-distribution. In such cases the plus-distribution can artificially inflate rounding errors. As this problem arises from rounding errors due to cancellations large enough to exhaust a typical double type variable’s width (i.e. the number of stored digits), the solution is of course to use data types storing more digits. For our approach the problem can only appear in the measurement function, as all other factors appearing are free of large cancellations. The best and most efficient solution is therefore to write a program that can use (slow) multi-precision arithmetic for the calculations related to the measurement functions (and only if told to do so), and that evaluates everything else using (fast) double precision floating-point arithmetic.

Having written a dedicated wrapper for the master formulae, a vanilla run of our program (or a vanilla run of (py)SecDec using the master formulae) reveals problems with numerical convergence that can be traced back to the appearance of square-root and logarithmic divergences. While these are analytically integrable divergences, a numerical approach that relies on sampling the integrand function runs into problems, particularly if the numerical integration involves adaptive variance-reduction techniques. Creative substitutions are again the solution to this problem. In particular, we can remove integrable divergences at both ends of a given integration over a variable $x\in[0,1]$ by substituting

[TABLE]

with suitably chosen parameters $i,j\geq 1$. Values of $i=2,4,\ldots$ then remove logarithmic, square-root, …divergences at $x=0$ while $j$ does so at $x=1$. As an example, consider

[TABLE]

where we used $i=2$ and $j=4$. The expression on the right-hand side is thus no longer divergent at either endpoint of the integration domain after the substitution.

This performance improvement comes at a minor cost: In addition to the functions $F_{A,B}$ and the parameter $n$, a second parameter $m$ is needed that can easily be derived from the functions $F_{A,B}$. As its significance is purely performance-related, we explain this in more detail in the SoftSERVE manual.

Applying these substitutions obscures the original definition of the observable and what ultimately serves as input into the program. We therefore inject a wrapping layer of substitutions between the user input and the numerics, to allow a user to input the relevant data using the physical parametrisation of the measurement function in terms of the variables $a$, $b$, $y$, $t_{k}$, $t_{l}$ and $t_{kl}$, which have physical meaning and allow more thorough checking of the input for consistency. In addition, it is easy to derive the one-loop input from the two-loop input in this form using the infrared-safety constraints (15) and (17), which we also exploit and implement.

An important boost in performance furthermore arises due to the fact that all divergences giving rise to a pole are sourced from monomials like $y^{-1+n\epsilon}$,

[TABLE]

This feature means that the leading pole coefficients need only to be integrated over domains with reduced dimensionality, as the delta distributions render some variables spurious. The full integrands appearing at each order are thus sums over expressions that depend on a reduced set of variables. That sum then again depends on all variables141414As a pedagogical example, consider a function $f(x,y)=g(x)+h(y)$ that depends on two variables, although its summands $g(x)$ and $h(y)$ depend only on one., and in SoftSERVE we simply relabel those variables to reduce the dimensionality, which can speed up the integrations significantly.

Finally, since we know exactly what SoftSERVE will be used for, we can supply it with scripts to facilitate the computations:

• •

execsftsrv runs the set of all colour structures for a given observable,

• •

sftsrvres sums different contributions to the same colour structure,

• •

laprenorm automates the renormalisation procedure outlined in Section 4,

• •

fourierconvert allows us to treat Fourier-space soft functions.

All scripts are supplied in the initial SoftSERVE ​0.9 release in a version without support for explicit calculation of the $C_{F}^{2}$ contribution, i.e. a version designed for observables obeying NAE. These scripts are therefore indexed by a suffix NAE. Following an upcoming publication dealing with uncorrelated emissions BRT and the subsequent release of SoftSERVE ​1.0, unsuffixed scripts that run the calculations for all colour structures will be supplied as well.

It is crucial to note, however, that we are not limited in scope to observables obeying NAE. The $C_{F}T_{F}n_{f}$ and $C_{F}C_{A}$ colour structures for NAE-breaking observables can still be computed using the initial release of SoftSERVE, and the scripts above can still be used for this purpose. The $C_{F}^{2}$ piece must then, however, be calculated by some other means, as the results for $C_{F}^{2}$ presented by the SoftSERVE ​0.9 release (or by using the NAE scripts in any released version) will simply not be correct – they are merely proportional to the square of the one-loop result.

The SoftSERVE package also comes supplied with a range of template observables serving as examples for an inexperienced user. Their salient features are described in the manual, and they include the observables presented in Section 7. Additionally, two sets of integrator settings are supplied, the standard and precision settings. For almost all cases the standard settings will provide sufficient accuracy; they represent the recommended choice of settings. The precision settings are designed to increase the achieved accuracy by roughly one order of magnitude, but operate at the edge of SoftSERVE’s capabilities. They are therefore more prone to unexpected errors and more sensitive to suboptimal coding. The user is of course free to choose his/her own custom settings, as the input files provide access to the relevant flags and options of the Cuba integrators.

As a final comment, the fourierconvert script allows us to calculate real parts of Fourier-space soft functions. As referenced in Section 2, we can treat observables with positive real part, which includes purely imaginary Fourier-space soft functions, where we assume an infinitesimal positive real part to make the analytic integration convergent. The numerical implementation is, on the other hand, limited to strictly real measurement functions, which would exclude the purely imaginary Fourier-space soft functions. The precise nature of these problems and the analytic workaround are delegated to Appendix B, where we find that the real part of a Fourier-space soft function can be calculated from the absolute value of the Fourier-space measurement function, after a reshuffling of the Laurent series. The fourierconvert script performs this reshuffle, taking the correct propagation of uncertainties into account.

7 Results

In this section we collect our results for multiple dijet soft functions that have been calculated using our novel algorithm and its implementation in SoftSERVE. The relevant resummation ingredients for SCET-1 soft functions were defined in Section 4.1, and we present our results in the form

[TABLE]

As all soft functions we consider in this section obey the NAE theorem, we have $\gamma_{1}^{C_{F}}=0$, $c_{2}^{C_{F}}=1/2(c_{1}^{C_{F}})^{2}$, and the missing NNLO coefficients can thus entirely be determined from the correlated-emission contribution.

For SCET-2 soft functions we similarly decompose the collinear anomaly exponent

[TABLE]

and we again have $d_{2}^{C_{F}}=0$ for soft functions that obey the NAE theorem. As explained in Section 4.2, we also obtain results for the bare soft remainder function $W_{0}^{S}(\tau)$, which is however meaningless without its collinear counterpart $W_{0}^{C}(\tau)$ in the collinear anomaly framework. We therefore restrict our attention to the collinear anomaly exponent for SCET-2 observables in the following.

The numbers we present in this section can be reproduced using the template files that are provided in the SoftSERVE distribution (bearing in mind that Monte Carlo integrations yield statistical predictions). Whenever we quote numbers, we used the precision setting for the numerical integrations unless mentioned otherwise. For the plots we used the standard setting, since the uncertainties are anyway not visible on the scale of the plots. The standard setting yields runtimes of a few minutes on a single 8-core machine, and it usually produces results with percent accuracy. Calculations using the precision setting, on the other hand, typically run for a few hours and they yield numbers with per mille accuracy or better. The Cuba library supports parallelisation, and so these run-time estimates are highly dependent on the number of available cores. As the SoftSERVE manual contains instructions about how to set up the corresponding input files, we only present the relevant input functions along with the numerical results in this section.

For many of the observables we consider in this section, the soft function has been determined to NNLO before, either by an explicit calculation or via a fit to a fixed-order code. The available results provide useful cross checks for our code, and they allow us to assess its numerical performance. For those results that were not available in the literature, we checked our SoftSERVE numbers with independent (py)SecDec runs. In some cases the uncertainty estimate provided by the integrator is smaller than the rounding error introduced by the truncation to six digits, and in these cases we simply give an uncertainty for the last digit (i.e. an error estimate of $10^{-6}$ for a result of $9.8696$ will show up as $9.8696(1)$).

7.1 $e^{+}e^{-}$ event shapes

The formalism we have developed in this paper applies to soft functions that are defined in terms of two light-like directions, and $e^{+}e^{-}$ event-shape variables that obey a hard-jet-soft factorisation theorem in the dijet limit are primary examples that fall into this class. In the following we present results for popular event shapes like C-parameter, thrust, and total jet broadening, but we also study the less known event shape angularities and a particular example of a double-differential observable.

C-parameter

The C-parameter was one of the observables that we already introduced in Section 2. Starting from the definition

[TABLE]

the relevant input functions were given in Table 2, with $n=1$, $f(y_{k},t_{k})=1/(1+y_{k})$, and151515Throughout this section we suppress the angular variables in the arguments of the two-emission measurement function if the observable does not depend on any of these angles. We also remind the reader that the expression for the measurement function in region B is in general not unique, due to the freedom in the definition (53).

[TABLE]

We then obtain for the C-parameter soft function

[TABLE]

where the available results from Hoang:2014wka are listed in the square brackets. While the first four numbers are known analytically, the last two numbers were obtained in Hoang:2014wka via a fit to the EVENT2 generator. Our numbers confirm these fit values, but they are significantly more accurate. For the C-parameter, we could actually derive analytic results following the strategy that we used for the derivation of anomalous dimensions in Bell:2018vaa , yielding

[TABLE]

which confirms our SoftSERVE numbers at the $2\sigma$ level.

Thrust

Thrust is the canonical event shape, which on the level of the soft function is defined as

[TABLE]

Thrust is again a SCET-1 observable with $n=1$, $f(y_{k},t_{k})=1$, and

[TABLE]

For the thrust soft function we find

[TABLE]

which is in agreement with the analytic results from Kelley:2011ng ; Monni:2011gb .

From the numerical perspective, observables like thrust are more problematic than the C-parameter because of the $1/y$ divergence in the second term of (88), which is formally cut off by the step function. As the limit $y\rightarrow 0$ is a singular limit of the matrix element, the contribution from this phase-space region is enhanced and it therefore sources the numerical instability. As a consequence we find that observables like thrust are more prone to producing ill-defined results, and their uncertainty estimates are often less trustworthy. The presence of a plus-distribution in the subtraction also means that this problem cannot be solved, but merely mitigated, with some ideas laid out in the SoftSERVE manual.

Jet broadening

The first SCET-2 observable in our list is total jet broadening, which has been studied within SCET in Becher:2011pf ; Chiu:2012ir ; Becher:2012qc . When the broadening is measured with respect to the thrust axis, it is well known that soft recoil effects complicate the resummation, which could be circumvented by choosing a different reference axis Larkoski:2014uqa . We will come back to a recoil-free definition of jet broadening later, but for the moment we simply switch off the recoil effects in the original thrust-axis definition since this allows us to compare our results to an existing NNLO calculation. We thus start from the definition

[TABLE]

which yields $n=0$ as required for a SCET-2 observable, along with $f(y_{k},t_{k})=1/2$ and

[TABLE]

Using SoftSERVE we then find for the collinear anomaly exponent

[TABLE]

which is in excellent agreement with the analytic results from Becher:2012qc .

Angularities

Angularities represent a generalisation of the thrust and broadening event shapes that depend on a continuous parameter $A$. According to their standard definition, the angularities are measured with respect to the thrust axis, and they are defined as

[TABLE]

which reduces to thrust for $A=0$, and is proportional to total broadening for $A=1$. For values of $A<1$ considered here, the angularities fall into the SCET-1 class with $n=1-A$, $f(y_{k},t_{k})=1$, and

[TABLE]

One easily verifies that these expressions reduce to those of thrust in the corresponding limit, $A\rightarrow 0$. The angularity soft function has been computed to NLO in Hornig:2009vb , where it was found that $\gamma_{0}^{S}(A)=0$ and $c_{1}^{C_{F}}(A)=-\pi^{2}/(1-A)$. Our SoftSERVE numbers indeed confirm these expressions and our results for the two-loop soft anomalous dimension and the finite term of the renormalised soft function are displayed as a function of the parameter $A$ in Figure 5. While our results for the constants $c_{2}^{C_{A}}$ and $c_{2}^{n_{f}}$ are new (preliminary results were reported in Bell:2015lsf ), we already derived the two-loop soft anomalous dimension in Bell:2018vaa . Our results have in fact already been utilised to resum the angularities distribution to NNLL accuracy in Procura:2018zpn and to NNLL*′* accuracy in (Bell:2018gce, ), and they thus represent the first phenomenological application of SoftSERVE.

As mentioned in the previous paragraph, the angularities share a numerical instability with thrust. We therefore find that some SoftSERVE runs using the precision setting produce ill-defined results, and the error bars are in general also less robust than for other observables.

Broadening-axis angularities

Similar to the total broadening event shape, the thrust-axis angularities are subject to soft recoil effects which become increasingly more important in the limit $A\to 1$. One way of avoiding these complications consists in defining the angularities with respect to the broadening axis rather than the thrust axis. Following Larkoski:2014uqa , we then start from

[TABLE]

where we omitted a factor $2^{B-1}$, which makes the expression for the one-loop soft anomalous dimension easier.161616According to (4), a change in the normalisation of $\omega(\{k_{i}\})$ always results in a rescaling of the Laplace parameter $\tau$. This in turn modifies the numerical value of the non-cusp anomalous dimension in the RGE (57) and the $d_{i}$ coefficients in (68). We are, however, free to turn the argument around and choose the normalisation of $\omega(\{k_{i}\})$ such that the expression of the one-loop anomalous dimension becomes simple. Thrust corresponds to $B=2$ in this notation, and we are interested in values $B\leq 2$ that extend beyond the broadening limit $B=1$. The broadening-axis angularities are thus characterised by $n=B-1$, $f(y_{k},t_{k})=(1+y_{k})^{1-B/2}$, and

[TABLE]

where we defined

[TABLE]

The considered soft function is known to NLO from Larkoski:2014uqa , and our results confirm these findings with $\gamma_{0}^{S}(B)=0$ and $c_{1}^{C_{F}}(B)=(B^{2}-3B-1)/(B-1)\,\pi^{2}/3$. Our NNLO results, on the other hand, are new and they are shown in Figure 6 as a function of the parameter $B$. While we cannot run SoftSERVE in the SCET-1 mode for $B=1$, the plots illustrate that we can compute the anomalous dimension and the matching correction sufficiently close to the SCET-2 limit – done here using $B=0.99$ and $B=1.01$ – and we can even interpolate between those points to convert the soft anomalous dimension into the SCET-2 anomaly exponent using the formulae provided in Bell:2018vaa .

Hemisphere masses

We finally consider the hemisphere soft function Hornig:2011iu ; Kelley:2011ng , which represents our first example of a double-differential observable. In this case, we take two Laplace transformations with respect to the hemisphere masses $M_{L}$ and $M_{R}$, and we denote the respective Laplace variables by $\tau_{L}$ and $\tau_{R}$. We furthermore introduce the variables

[TABLE]

such that the dependence on the dimensionful variable $\tau_{1}$ factorises by construction. The observable, on the other hand, now depends on the second Laplace variable $\tau_{2}$ via

[TABLE]

with $\bar{\tau}_{2}=1-\tau_{2}$, and we have furthermore adjusted the normalisation for later convenience – see footnote 16. The observable is thus described by $n=1$, $f(y_{k},t_{k};\tau_{2})=\sqrt{\tau_{2}/\bar{\tau}_{2}}$, and

[TABLE]

from which we recover the thrust expressions for $\tau_{2}=1/2$. The analysis of the hemisphere soft function is, moreover, complicated by the fact that the definition in (99) is not symmetric under $n$-$\bar{n}$ exchange, since the roles of $M_{L}$ and $M_{R}$ – and hence $\tau_{2}$ and $\bar{\tau}_{2}$ – are interchanged under this symmetry. This seems to be in conflict with assumption (A4), but we can simply restore this symmetry by substituting

[TABLE]

on the level of the bare soft function171717This procedure illustrates the comment made in Section 5 concerning the emulation of a broken $n$-$\bar{n}$-symmetry, with $S_{0}(\tau_{1},\bar{\tau}_{2})$ effectively representing the contribution from regions C and D.. As the soft anomalous dimension in (62) and the matching correction in (63) depend only linearly on the $C_{F}C_{A}$ and $C_{F}T_{F}n_{f}$ pole coefficients, the average in (101) can actually be directly implemented on the level of these quantities as well.

Due to the particular normalisation in (99), the anomalous dimension of the hemisphere soft function matches the one for thrust, and we display the finite terms of the two-loop soft function as a function of the Laplace variable $\tau_{2}$ in Figure 7 (the interval $\tau_{2}\in[0.5,1]$ is just a mirror image of $\tau_{2}\in[0,0.5]$ due to the symmetrisation procedure). In addition, we have Laplace-transformed the analytic expressions from Kelley:2011ng , and the result is shown in Figure 7 by the solid lines. As can be seen from the plots, we observe perfect agreement for this double-differential observable, and we of course also confirm the corresponding one-loop expression with $c_{1}^{C_{F}}(\tau_{2})=-8\,\mathrm{artanh}^{2}(1-2\tau_{2})-\pi^{2}$.

As a generalisation of thrust, the hemisphere soft function is prone to numerical instabilities, and for two out of twenty runs that were needed to produce the plots in Figure 7, two resulted in ill-defined expressions for the default setting of the numerical integrator. This problem could, however, be remedied by increasing the value of the Divonne border parameter from $10^{-8}$ to $10^{-6}$. While doing so, we verified that the variation of this parameter does not introduce systematic uncertainties that are relevant at the quoted accuracy.

7.2 Hadron-collider observables

Another important class of soft functions that fall into the considered ‘dijet’ category are hadron-collider observables for which central jets are vetoed, either via a kinematical restriction or an explicit jet algorithm. As the latter are usually not consistent with NAE, we postpone their discussion to a future study BRT , and we focus on soft functions that are relevant for soft-gluon resummation, transverse-momentum resummation, and hadronic event shapes in the following. For hadron-collider soft functions, only $q\bar{q}$–initiated processes are of the type (1), whereas the Wilson lines for other channels refer to different colour representations. The soft function does, however, obey Casimir scaling to the considered order, and the results can therefore easily be translated to other channels by rescaling the expressions in (81) and (82) with $C_{i}/C_{F}$, where $C_{i}$ is the Casimir operator of the said channel.

Threshold Drell-Yan production

The production of a lepton pair $pp\rightarrow l_{1}l_{2}X$ at threshold represents the simplest hadron-collider soft function we can consider. We already encountered its definition in Section 2, where we found that

[TABLE]

which yields $n=-1$, $f(y_{k},t_{k})=1+y_{k}$, and

[TABLE]

Using SoftSERVE, we then find

[TABLE]

which nicely agrees with the analytic NNLO results from (Belitsky:1998tc, ).

W-production at large transverse momentum

We next consider the soft function for $W$-production at large transverse momentum from Becher:2012za . Although this is a soft function that involves more than two Wilson lines, it represents de facto a dijet soft function, since the gluon attachments to the Wilson line $S_{n_{J}}$ vanish up to NNLO and we are furthermore free to choose $n_{1}\cdot n_{2}=2$ along with $n_{1}\cdot n_{J}=n_{2}\cdot n_{J}=2$ due to rescaling invariance of the Wilson lines Becher:2012za . The vector $n_{J}^{\mu}$ now introduces a non-trivial angular dependence with

[TABLE]

where $\theta_{i}=\sphericalangle(\vec{n}_{J}^{\perp},\vec{k}_{i}^{\perp})$. Up to NNLO the soft function is thus characterised by $n=-1$, $f(y_{k},t_{k})=1+y_{k}-2\sqrt{y_{k}}(1-2t_{k})$, and

[TABLE]

and we obtain

[TABLE]

which is again in agreement with the analytic results from Becher:2012za .

Exclusive Drell-Yan production

The soft function for exclusive Drell-Yan production Li:2011zp is another example of a double-differential observable. Due to rescaling invariance of the Wilson lines, the position-space soft function can in this case only depend on $\tau_{1}=i/2\,\sqrt{x_{+}x_{-}}$ and $\tau_{2}=\sqrt{x_{T}^{2}/x_{+}x_{-}}$. Whereas the dependence on the dimensionful variable $\tau_{1}$ factorises, the observable then becomes a non-trivial function of the variable $\tau_{2}$ with

[TABLE]

The soft function thus generalises the two preceding examples, and we recover the threshold Drell-Yan soft function for $\tau_{2}=0$, while the one for $W$-production at large transverse momentum corresponds to $\tau_{2}=1$. In terms of our parametrisations, we find $n=-1$, $f(y_{k},t_{k};\tau_{2})=1+y_{k}-2\tau_{2}\sqrt{y_{k}}(1-2t_{k})$, and

[TABLE]

It turns out that the respective anomalous dimension is independent of $\tau_{2}$, and it can therefore be read off from the two previous examples. The finite term of the renormalised NNLO soft function is, on the other hand, shown in Figure 8 together with the numbers from Li:2011zp . We again find perfect agreement with the known analytic results, which is of course also true for the one-loop constant $c_{1}^{C_{F}}(\tau_{2})=4\,\mathrm{Li}_{2}(\tau_{2}^{2})+\pi^{2}/3$.

Transverse-momentum resummation

The soft function for transverse-momentum resummation in Drell-Yan production was the third example we considered in Section 2. There we saw that

[TABLE]

where the factor of $i$ reflects a Fourier transformation. The measurement function is thus purely imaginary in this case and, unlike in the preceding example, we cannot simply factor out an imaginary unit along with the Laplace variable $\tau$, since $\omega^{p_{T}}(\{k_{i}\})$ would in this case not be bounded to be positive, which would contradict assumption (A2). Nevertheless we argued in Section 5 that we can compute Fourier-space soft functions by using the absolute value of the naive measurement function. The required input functions for the transverse-momentum soft function are therefore given by $n=0$, $f(y_{k},t_{k})=2|1-2t_{k}|$, and

[TABLE]

Further instructions for the computation of Fourier-space soft functions can be found in Appendix B. Applying the fourierconvert script before renormalisation then leads to

[TABLE]

The slightly reduced accuracy is due to the appearance of integrable logarithmic divergences in the bulk of the integration region, since the absolute value can vanish for non-trivial angular configurations. Nevertheless, the agreement with the known results from Becher:2010tm ; Gehrmann:2014yya is satisfactory.

Transverse thrust

We finally consider the soft function for the hadronic event shape transverse thrust Becher:2015gsa . As this example involves four light-like Wilson lines, the computation of the full NNLO soft function clearly falls outside the scope of the present paper. It was shown, however, in Becher:2015lmy that the underlying anomalous dimensions can be reconstructed from the information on dijet soft functions, and we therefore concentrate here on the computation of the anomalous dimensions for transverse thrust.

To this end, we split the $2\rightarrow 2$ process into two toy processes for $e^{+}e^{-}\to q\bar{q}$ and $q\bar{q}\to e^{+}e^{-}$ scattering. We first consider the latter, which according to Becher:2015gsa gives rise to

[TABLE]

where $c_{0}=e^{4G/\pi}$ depends on Catalan’s constant $G\simeq 0.915966$, and the vector $\vec{n}_{\perp}$ singles out a direction in the plane transverse to the beam direction. We thus find that the toy observable is of SCET-2 type with $n=0$, $f(y_{k},t_{k})=2c_{0}(1-|1-2t_{k}|)$, and

[TABLE]

Using SoftSERVE we then determine the corresponding collinear anomaly exponent, which yields

[TABLE]

where the numerical two-loop results in the square brackets are taken from Becher:2015lmy . Our results indeed confirm these numbers, but they are again significantly more precise.

The soft function for the second toy process starts from a similar definition Becher:2015gsa ,

[TABLE]

where $s=\sin\theta_{B}$ depends on the angle between the beam and the jet axis. Its decomposition into light-cone coordinates is, however, significantly more involved, since the components $k_{\top}^{\mu}$ transverse to the thrust axis differ from the components $k_{\perp}^{\mu}$ in (116) that are transverse to the beam axis – see Becher:2015gsa . We then find that this soft function is described by $n=1$ and

[TABLE]

where $c=\cos\theta_{B}$ and $c_{k}=\cos\theta_{k}=1-2t_{k}$. The expressions for the two-emission measurement functions are rather lengthy with

[TABLE]

and

[TABLE]

where we remind the reader that the last expression is not unique, due to the freedom in the definition (53).

The complicated structure of the measurement functions clearly exhibits that an analytic calculation is very challenging – if not out of reach – for this class of soft functions. Yet our numerical approach also suffers from large cancellations between the square roots and the absolute values for small values of the parameter $y$. This causes major problems for the numerical evaluation due to rounding errors in floating point computations. As this is mostly an implementation issue, we refer the user to the manual of SoftSERVEfor a detailed discussion. The consequence is that SoftSERVE must be run using multi-precision variables provided by the boost boost or GMP/MPFR GMP/MPFR libraries for this observable. This significantly slows down the program, and hence our transverse thrust runs used reduced accuracy settings, as specified in more detail in the SoftSERVE manual.

Although the measurement functions for this observable depend on the angle $\theta_{B}$ between the beam and jet axes, it turns out that the anomalous dimension is independent of this angle, and we obtain

[TABLE]

where the results of Becher:2015gsa were obtained via a fit to the EVENT2 generator.

As explained at the beginning of this section, the finite terms of the corresponding soft functions are not useful, as the split into $e^{+}e^{-}\to q\bar{q}$ and $q\bar{q}\to e^{+}e^{-}$ observables is only justified on the level of the anomalous dimension and the collinear anomaly exponent. So while the SoftSERVE runs of course produce numbers for the finite terms as well, they cannot be used for any meaningful interpretation in this case.

8 Conclusions

We have presented a novel formalism that allows for an automated calculation of soft functions that are defined in terms of two light-like Wilson lines. Our method is based on a universal phase-space parametrisation which we use to isolate the implicit divergences of the phase-space integrals, and which allows us to perform the expansion in the various regulators directly on the integrand level. The remaining integrations can then be performed numerically, and we developed a C++ package called SoftSERVE which precisely performs these integrations using the Divonne integrator of the Cuba library.

Our method is currently restricted to the correlated-emission contribution, which is sufficient for the computation of NNLO soft functions that obey the NAE theorem. As exemplified by the large number of results in Section 7, our method is fairly general and it relies on a handful of assumptions that we discussed in detail in Sections 2 and 5. Further extensions of our formalism that deal with the uncorrelated-emission contribution Bell:2018jvf and with more than two light-like directions including non-back-to-back kinematics Bell:2018mkk are currently in progress.

With the publication of the present paper, we release SoftSERVE ​0.9, which is publicly available at https://softserve.hepforge.org/. The current version of SoftSERVE not only allows for the calculation of bare dijet soft functions, but it also provides scripts for their renormalisation according to the conventions we introduced in Section 4. In this paper we refrained from presenting the technical aspects of the numerical implementation, which are explained in detail in the SoftSERVE manual.

Among the plethora of results we obtained in Section 7, our NNLO numbers for the thrust-axis and broadening-axis angularities are new, and we for the first time obtained the C-parameter soft function in an analytic form. More generally, our results for the two-loop anomalous dimensions and the two-loop finite terms are required for NNLL and NNLL*′* resummations within SCET, respectively. We believe that SoftSERVE can open the path for high-precision SCET resummations in the future, as has recently been illustrated for the $e^{+}e^{-}$ event shape angularities in Bell:2018gce .

Yet the process of automating resummation within SCET requires many further ingredients beyond the automated calculation of soft functions. At present there exists an automated NLL resummation code that is based on the coherent-branching formalism, CAESAR Banfi:2004yd , which has been extended to NNLL accuracy for $e^{+}e^{-}$ event shapes in Banfi:2014sua . While a combination of the CAESAR approach with methods from SCET is currently under investigation Bauer:2018svx , we believe that an alternative approach that is purely based on effective-field-theory techniques could be a valuable addition. We think that SoftSERVE could provide one essential pillar for such an automated resummation code in SCET.

Acknowledgements.

We would like to sincerely thank the developers of SecDec for their support throughout the duration of this work, and in particular for the dedicated development of the multi-regulator functionality in pySecDec, which facilitated crosschecks for SCET-2 observables. G.B. is supported by the Deutsche Forschungsgemeinschaft (DFG) within Research Unit FOR 1873, and R.R. is supported by the Swiss National Science Foundation (SNF) under grants CRSII2$\_$160814 and 200020$\_$182038. J.T. and R.R. acknowledge travel support from the University of Oxford Department of Physics and Universität Siegen. J.T. also acknowledges research and travel support from the Senior Scholarship Trust of Hertford College (Oxford) and DESY Hamburg.

Appendix A Equality of one-emission and two-emission values of $n$

As mentioned in Section 2, the leading scaling in the variables $y_{k}$ and $y$ is the same between the one-emission and two-emission measurement functions, see (10) and (12). As we will show in this appendix, this is a non-trivial finding for generic observables, in particular for those that violate the NAE theorem.

To this end, we first rewrite the two-emission measurement function in the form

[TABLE]

where $\omega_{c}(\{k,l\})$ encodes the effects that break NAE. Whenever the emission with momentum $k^{\mu}$ becomes soft, infrared safety implies that the measurement function reduces to $\omega(\{l\})$, see (2.2), and hence the NAE-violating term must satisfy $\omega_{c}(\{0,l\})=0$. Similarly, we obtain $\omega_{c}(\{k,0\})=0$ when $l^{\mu}$ becomes soft, and $\omega_{c}(\{k,\alpha k\})=0$ with $\alpha>0$, which arises when the two emissions become collinear to each other.

Since both $y_{k}=k_{+}/k_{-}\sim y$ and $y_{l}=l_{+}/l_{-}\sim y$ in the parametrisation (11), the leading scaling in the variable $y$ in the first two terms of (121) is equal to the one-emission case. The question therefore becomes whether the NAE-breaking term $\omega_{c}(\{k,l\})$ can enforce a different scaling of the observable.

The answer lies in the factorisation theorem (3), and the fact that the poles must cancel between the bare hard, jet and soft functions. As the hard function is the same for all dijet observables, the crucial cancellation is the one between the collinear and soft functions. The role of the parameter $n$ is two fold in this context: First, it determines the coefficient of the leading pole in the soft function, as the monomial $y^{-1+n\epsilon}$ gives rise to a term $\delta(y)/n\epsilon$. And second, the analysis in Section 2.3 showed that the parameter $n$ is related to the power counting of the momentum modes in the effective theory. Following this analysis, we thus know what the scaling of the three sectors — soft, collinear and anti-collinear — is for two emissions, as the form of the factorisation theorem implies that no new modes appear at higher orders. As the power counting of the momentum components is uniform in the soft sector, the form of the NAE-breaking term $\omega_{c}(\{k,l\})$ is unconstrained, since consistency demands that its mass dimension is the same as that of the one-emission terms, which immediately puts the three terms in (121) at the same order in the power counting. This does not hold, however, in the collinear sectors.

We therefore consider the collinear sector with scaling $p_{c}^{\mu}=Q(1,\lambda^{2/(n+1)},\lambda^{1/(n+1)})$ more closely. As we have seen in (18), the measurement function for one collinear emission scales in the soft limit $z_{k}\to 0$ as

[TABLE]

Based on power-counting and dimensional arguments, the collinear analog of the NAE-breaking term then similarly scales in the soft limit as

[TABLE]

where $q$ can stand in for both $k$ or $l$ (i.e. $\sqrt{k_{+}l_{-}}$ and $\sqrt{k_{+}k_{-}}$ both count as $\sqrt{q_{+}q_{-}}$ here). The question then is whether $n^{\prime}$ can differ from $n$. There are three possibilities:

1. $n=n^{\prime}$:

The three terms in (121) contribute in the soft limit to the jet function, and its leading poles are determined by the one-emission value of $n$. To ensure the cancellation of the leading poles, the soft function must follow suit and its two-emission measurement function must therefore scale as $y^{n/2}$. 2. $n<n^{\prime}$:

In this case the NAE-breaking term is suppressed in the jet function in the soft limit. As the jet function does not see this term in the soft region, its leading poles match those of an auxiliary observable (the NAE version of the observable), for which the leading poles are again determined by the one-emission value of $n$ as described above. 3. $n>n^{\prime}$:

In this case the NAE-breaking term dominates over the one-emission terms in the soft limit. But the observable would then not be infrared safe, which violates one of our underlying assumptions.

We therefore conclude that the factorisation theorem (3), in particular the absence of additional momentum modes, enforces the equality of the one-emission and two-emission values of $n$.

Appendix B Fourier-space soft functions

In Section 5 we explained that the numerical evaluation of Fourier-space soft functions using SoftSERVE may require a workaround, since SoftSERVE assumes that the measurement function is strictly real and non-negative. In this appendix we describe the details of this workaround.

To begin with, we consider a one-emission measurement function of the form

[TABLE]

where the factor of $i$ arises because of the Fourier transformation, $\tau$ is the respective Fourier variable and the function $\widetilde{f}(y_{k},t_{k})$ is assumed to be real-valued. If the latter is non-negative, the factor of $i$ can simply be factorised alongside $\tau$ outside the integral, and the SoftSERVE implementation proceeds in the standard form (threshold resummation for Drell-Yan production with $\tau=ix^{0}/2$ is of this type). In general, however, the function $\widetilde{f}(y_{k},t_{k})$ may well attain also negative values, as is the case e.g. for transverse-momentum resummation with $\widetilde{f}(y_{k},t_{k})=-2(1-2t_{k})$. As SoftSERVE only accepts non-negative measurement functions, it seems as if we cannot evaluate such functions with the existing code. There exists, however, a workaround as long as one is interested only in the real part of the soft function (which is usually the case).

To this end, we review the analytic steps of our derivation, and we split the calculation into two parts according to the two possible signs of the measurement function $\widetilde{f}(y_{k},t_{k})$. In each of these regions, we then factor out either $(+i)$ or $(-i)$ such that the remnant measurement function is non-negative. As the measurement function enters the master formula (30) as $f(y_{k},t_{k})^{2\epsilon+\alpha}$, we can write the soft function in the form

[TABLE]

where we indicated that S_{R}\big{[}f\big{]} in (30) depends on the function $f(y_{k},t_{k})=i\widetilde{f}(y_{k},t_{k})$, while S_{R}^{\pm}\big{[}f\big{]} depend on the magnitude of $\widetilde{f}(y_{k},t_{k})$. The latter are, moreover, restricted to an integration domain in which $\widetilde{f}(y_{k},t_{k})$ is either positive (S_{R}^{+}\big{[}|\widetilde{f}|\big{]}) or negative (S_{R}^{-}\big{[}|\widetilde{f}|\big{]}), and they are not accessible in our formalism, since we always assume that the integrations are performed over the full domain. But if we decompose the soft function into its real and imaginary parts, we obtain

[TABLE]

where we used the fact that S_{R}^{\pm}\big{[}|\widetilde{f}|\big{]} are real-valued, and that S_{R}^{+}\big{[}|\widetilde{f}|\big{]}+S_{R}^{-}\big{[}|\widetilde{f}|\big{]}=S_{R}\big{[}|\widetilde{f}|\big{]}, i.e. the sum over two complementary integration domains gives back the full integral. We thus see that the observable-dependent phase-space restrictions drop out in the real part of the soft function. In other words we simply have to evaluate the soft function with $\widetilde{f}(y_{k},t_{k})$ replaced by its absolute value, which is yet to be multiplied by a regulator-depend factor.

The same arguments apply to the mixed real-virtual interference S_{RV}\big{[}f\big{]} and the double real-emission contribution S_{RR}\big{[}F\big{]}, which depend on different powers of the measurement functions $f(y_{k},t_{k})^{4\epsilon+\alpha}$ and $F(a,b,y,t_{k}^{+},t_{l},t_{kl})^{4\epsilon+2\alpha}$, respectively, see (36) and (50). The real part of a bare Fourier-space soft function then becomes

[TABLE]

If one is interested in calculating the real part of a Fourier-space soft function, one should thus run SoftSERVE using the absolute value of the measurement function (which is by construction real and non-negative). The individual pieces of this calculation then need to be multiplied with different regulator-dependent factors, which reshuffle the coefficients in the Laurent expansion. The SoftSERVE package contains a script to perform this Laurent series reshuffle keeping track of the error bars — its usage is explained in the manual. As an application of this technique we consider the soft function for transverse-momentum resummation in Section 7.