Two-loop $N$-jettiness soft function for $pp \to 2j$ production

TL;DR

This paper computes the two-loop soft function for various jet production processes, enabling precise NNLO predictions using the N-jettiness subtraction scheme.

Contribution

It provides the missing two-loop soft function calculation crucial for NNLO accuracy in N-jettiness based predictions.

Findings

Two-loop soft function calculated for pp → jj, e+e- → 4j, ep → 3j.

Enables full NNLO predictions for these processes.

Supports N-jettiness subtraction scheme implementation.

Abstract

We calculate the two-loop soft function for $pp \to jj$, $e^+e^- \to 4j$ or $e p \to 3j$ when the $N$-jettiness observable is measured. The result presented here makes up the necessary piece for realizing the full next-to-next-to-leading order predictions for these processes using the $N$-jettiness subtraction scheme.

Figures

Equations

$$\mathrm{d}\sigma ({\tau }_{\mathrm{cut}})=BB\mathrm{Tr}[S\cdot H]\prod _{i=1}{J}_i+\mathcal{O}\left({\tau }_{\mathrm{cut}}\right),$$

$$\tau =\sum _i\min \left(n_1\cdot q_i,n_2\cdot q_i,\dots ,n_N\cdot q_i\right),$$

$$S(\tau )=\prod _i^n\int \frac{{\mathrm{d}}^d{q}_i}{(2\pi {)}^{d-1}}{\delta }^{+}(q_i^2)\mathcal{S}(\{q_i\}){\Theta }_{N,i}(\{q_i\}),$$

$${\Theta }_{N,i}(\{q_i\})=\sum _{j=1}^N\delta (\tau -\sum _i^n{n}_{j_i}\cdot q_i)\prod _{k_i\ne j_i}^{N-1}\theta (n_{k_i}\cdot q_i-n_{j_i}\cdot q_i),$$

$${\int }_0^1\prod _i\mathrm{d}{x}_i{x}_i^{-1-a_iϵ}F(\{x_i\})={\int }_0^1\prod _i\mathrm{d}{x}_i\left[-\frac{1}{a_iϵ}\delta (x_i)+{\left(\frac{1}{x_i}\right)}_{+}+\dots \right]F(\{x_i\}),$$

$$k^{\mu }=\frac{n_j^{\mu }}{n_{ij}}{k}^{+}+\frac{n_i^{\mu }}{n_{ij}}{k}^{-}+k_{\perp }^{\mu },\mathrm{denoted}\mathrm{as}{k}^{\mu }=(k^{+},k^{-},k_{\perp }^i),$$

$$n_k^{\mu }=(n_k^{+},n_k^{-},0,n_k^{\perp };0),$$

$$n_l^{\mu }=(n_l^{+},n_l^{-},n_l^{\perp }{s}_{{\varphi }_l},n_l^{\perp }{c}_{{\varphi }_l};0),$$

$$q_1^{\mu }=(q_1^{+},q_1^{-},q_1^{\perp }{s}_{{\varphi }_1}{c}_{{\alpha }_1},q_1^{\perp }{c}_{{\varphi }_1};q_1^{\perp }{s}_{{\varphi }_1}{s}_{{\alpha }_1},0),$$

$$c_{{\varphi }_l}=\frac{n_{ik}{n}_{jl}+n_{il}{n}_{jk}-n_{ij}{n}_{kl}}{2\sqrt{n_{ik}{n}_{il}{n}_{jk}{n}_{jl}}}.$$

$$\int \frac{{\mathrm{d}}^d{q}_1}{(2\pi {)}^{d-1}}{\delta }^{+}(q_1^2)F(q_1)\delta (\tau -n_i\cdot q_1)\theta (n_j\cdot q_1-n_i\cdot q_1)$$

$$q_1^{-}=\frac{\tau }{x_1},c_{{\varphi }_1}=\pi x_2,c_{{\alpha }_1}=1-2x_3.$$

$$q_1^{\mu }=(q_1^{+},q_1^{-},0,q_1^{\perp };0),$$

$$n_k^{\mu }=(n_k^{+},n_k^{-},-n_k^{\perp }{s}_{{\varphi }_1},n_k^{\perp }{c}_{{\varphi }_1};0),$$

$$n_l^{\mu }=(n_l^{+},n_l^{-},n_l^{\perp }{s}_{{\varphi }_l-{\varphi }_1}^{{\alpha }_1},n_l^{\perp }{c}_{{\varphi }_l-{\varphi }_1}^{{\alpha }_1};-n_l^{\perp }{s}_{{\varphi }_l}{s}_{{\alpha }_1},0),$$

$$q_2^{\mu }=(q_2^{+},q_2^{-},q_2^{\perp }{s}_{{\varphi }_2}{c}_{{\alpha }_2},q_2^{\perp }{c}_{{\varphi }_2};q_2^{\perp }{s}_{{\varphi }_2}{s}_{{\alpha }_2}{c}_{\beta },q_2^{\perp }{s}_{{\varphi }_2}{s}_{{\alpha }_2}{s}_{\beta },0),$$

$$q_1\cdot q_2=\frac{1}{n_{ij}}\left[q_1^{+}{q}_2^{-}+q_1^{-}{q}_2^{+}-2\sqrt{q_1^{+}{q}_2^{-}{q}_1^{-}{q}_2^{+}}{c}_{{\varphi }_2}\right],$$

$$\eta =\frac{{\left(\sqrt{q_1^{+}{q}_2^{-}}-\sqrt{q_1^{-}{q}_2^{+}}\right)}^2(1-x^{\prime })}{{\left(\sqrt{q_1^{+}{q}_2^{-}}-\sqrt{q_1^{-}{q}_2^{+}}\right)}^2+4\sqrt{q_1^{+}{q}_2^{-}{q}_1^{-}{q}_2^{+}}{x}^{\prime }},x^{\prime }={\sin }^2\left(\frac{\pi }{2}x\right),$$

$$\eta =x^{\prime },x^{\prime }={\sin }^2\left(\frac{\pi }{2}x\right),$$

$$S^{(1)}=\frac{{\alpha }_s}{4\pi }\sum _{i\ne j}{T}_i\cdot T_j\left(4C_{ij}^{(1,-1)}{\tilde{L}}^2+4C_{ij}^{(1,0)}\tilde{L}+\left[2C_{ij}^{(1,1)}+\frac{2}{3}{\pi }^2{C}_{ij}^{(1,-1)}\right]\right)\equiv \sum _{i\ne j}{T}_i\cdot T_j{S}_{ij}^{(1)},$$

$$S^{(2)}$$

$${\mathcal{S}}^{(1)}=-(4\pi {\alpha }_s{\mu }^{2ϵ})\times 2\times \sum _{i\ne j}{T}_i\cdot T_j{\mathcal{S}}_{ij}^{(1)}(q),$$

$${\mathcal{S}}_{ij}^{(1)}(q)=\frac{n_{ij}}{2n_i\cdot q{n}_j\cdot q}.$$

$${\mathcal{S}}_{RV,dipole}^{(2)}=\frac{(4\pi {\alpha }_s{\mu }^{2ϵ}{)}^2}{4{\pi }^2}\times \frac{(4\pi {)}^{ϵ}}{{ϵ}^2}\frac{{\Gamma }^4(1-ϵ){\Gamma }^3(1+ϵ)}{{\Gamma }^2(1-2ϵ)\Gamma (1+2ϵ)}\times C_A\sum _{i\ne j}{T}_i\cdot T_j{\left[{\mathcal{S}}_{ij}^{(1)}(q)\right]}^{1+ϵ},$$

$${\mathcal{S}}_{RV,tri.}^{(2)}$$

$${\mathcal{S}}_{q\bar{q}}^{(2)}=-\frac{1}{2}(4\pi {\alpha }_s{\mu }^{2ϵ}{)}^2{T}_R{N}_F\sum _{i\ne j}{T}_i\cdot T_j{\mathcal{J}}_{ij},$$

$${\mathcal{J}}_{ij}={\mathcal{I}}_{ii}+{\mathcal{I}}_{jj}-2{\mathcal{I}}_{ij},$$

$${\mathcal{I}}_{ij}=-\frac{2n_{ij}{q}_1\cdot q_2+n_i\cdot (q_1-q_2)n_j\cdot (q_1-q_2)}{2(q_1\cdot q_2{)}^2{n}_i\cdot (q_1+q_2)n_j\cdot (q_1+q_2)},$$

$${\mathcal{S}}_{gg,non-abe.}^{(2)}=-\frac{1}{2}{\left(4\pi {\alpha }_s{\mu }^{2ϵ}\right)}^2{C}_A\sum _{i\ne j}{T}_i\cdot T_j{\mathcal{T}}_{ij},$$

$${\mathcal{T}}_{ij}=2{\mathcal{S}}_{ij}-{\mathcal{S}}_{ii}-{\mathcal{S}}_{jj},$$

$${\mathcal{S}}_{ij}={\mathcal{S}}_{ij}^{s.o.}-\frac{1}{2}\frac{n_i\cdot q_1{n}_j\cdot q_2+n_i\cdot q_2{n}_j\cdot q_1}{n_i\cdot (q_1+q_2)n_j\cdot (q_1+q_2)}{\mathcal{S}}_{ij}^{s.o.}+(1-ϵ){\mathcal{I}}_{ij}^{\mathrm{I}}+2{\mathcal{I}}_{ij}^{\mathrm{II}},$$

$${\mathcal{S}}_{ij}^{s.o.}$$

Full text

I Introduction

The accumulation of the data set and the continuously increasing precision of the experimental analyses bring the physics program at the Large Hadron Collider (LHC) and other experiments the ability to search for new physics from the small deviations of the data from the Standard Model predictions. In order to fulfill such searches, the precision of the theoretical predictions from Quantum Chromodynamics (QCD) has to match the small experimental errors to reliably interpret the collider data. Currently, the only theoretical tool, with no known substitutes, enables us to make predictions for collider phenomenologies out of the first QCD principle is the perturbative calculation based on the systematic expansion of the strong coupling constant $\alpha_{s}$ or $\alpha_{s}L$ if large logarithmic correction $L$ is present. In recent years, tremendous efforts have been made in realizing the theoretical predictions for the physical processes at the colliders with a full control of the final state kinematics beyond the next-to-leading order (NLO) accuracy. On the fixed order side, significant progress has been made in the past few years, which includes achieving the fully differential gluon-gluon fusion Higgs ($H$) production Cieri:2018oms ; Dulat:2018bfe and the single jet production in deep-inelastic scattering Currie:2018fgr at next-to-next-to-next-to-leading order (N3LO) in $\alpha_{s}$ and the next-to-next-to-leading order (NNLO) calculations of the single inclusive jet or dijet productions Currie:2017eqf , $V/H+j$ processes Boughezal:2015dva ; Ridder:2015dxa ; Boughezal:2015ded ; Gehrmann-DeRidder:2017mvr ; Boughezal:2015dra ; Chen:2014gva ; Boughezal:2015aha , Higgs production through the vector boson fusion Cacciari:2015jma ; Cruz-Martinez:2018rod , single top Brucherseifer:2014ama ; Berger:2016oht and $t{\bar{t}}$ productions Czakon:2015owf ; Catani:2019iny at the LHC, the jet productions in deep-inelastic scattering Abelof:2016pby ; Currie:2017tpe , and the heavy flavor production in deep-inelastic neutrino scattering Berger:2016inr .

The rapid growth in the higher order predictions benefit from either the new higher loop calculations achieved or the theoretical schemes developed for dealing with the singularities in the real emissions. Recent achievements of the analytic two-loop five-point amplitude Abreu:2018aqd ; Chicherin:2018yne ; Chicherin:2018old highlight the current status of the loop calculations for the collider processes, and meanwhile new ideas are keeping on popping out Borowka:2018dsa ; Liu:2017jxz ; Liu:2018dmc ; Xu:2018eos . To handle the infrared (IR) singularities in the real corrections, various schemes have been proposed and successfully applied to the LHC phenomenology, such as the local subtraction schemes including the antenna subtraction GehrmannDeRidder:2005cm , the sector improved residue subtraction Czakon:2010td and the nested soft-collinear subtraction Caola:2017dug , as well as the physical observable based global schemes like the $q_{T}$-subtraction Catani:2007vq , inclusive jet mass subtraction Gao:2012ja and the $N$-jettiness subtraction Boughezal:2015dva ; Gaunt:2015pea . Although the global schemes are always criticized for the numerical instability, they have realized many interesting higher order predictions for the physical processes at the hadron colliders including the numerically most challenging Higgs production at N3LO Cieri:2018oms .

Among the local schemes, $N$-jettiness scheme is the known scheme being able to handle the IR singularities for generic jet processes at the hadron colliders. The scheme utilizes a threshold cut-off $\tau_{\rm cut}$ in the $N$-jettiness event shape $\tau$ Stewart:2010tn to distinguish between the NNLO $N$ jet (when $\tau<\tau_{\rm cut}$) and the NLO $N+1$ jet configurations ($\tau>\tau_{\rm cut}$). Below the $\tau_{\rm cut}$, a Soft-Collinear-Effective-Theory Bauer:2000ew ; Bauer:2000yr ; Bauer:2001ct ; Bauer:2001yt ; Bauer:2002nz approach based on the factorization theorem Stewart:2010tn

[TABLE]

is used to approximate the full QCD contribution. Here $B$ stands for the beam function due to the initial state energetic radiations while $J$’s are the jet functions for the final state collinear radiations. The hard function $H$ encodes the hard loop information. $S$ is the soft function for soft radiations. The trace is over the relevant color space. All the ingredients in the factorization theorem are universal except for the hard function $H$ which depends on the specific process under consideration. The difference between the factorization theorem and the full QCD is suppressed by power corrections in $\tau_{\rm cut}$, which have been studied extensively recently for the [math]-jettiness case Moult:2016fqy ; Boughezal:2016zws ; Moult:2017jsg ; Boughezal:2018mvf ; Ebert:2018lzn . The universal beam and the jet functions are all known analytically to two-loops Gaunt:2014xga ; Gaunt:2014cfa ; Becher:2006qw ; Becher:2010pd and recently some of them are even calculated at three-loop order Bruser:2018rad ; Banerjee:2018ozf ; Melnikov:2018jxb . As for the quark beam function, its NNLO longitudinal polarized counterpart is also known in the literature Boughezal:2017tdd . The calculation of the soft function is complicated by its explicit dependence on the full $N$-jettiness measure, therefore its analytic NNLO calculation is by far out of reach and people turn to numerical solutions. A generic framework was proposed in Boughezal:2015eha for calculating such soft functions with arbitrary numbers of jets and the $1$-jettiness soft functions in $pp$-collision were calculated using this scheme or similar approach for both the massless jet Boughezal:2015eha ; Campbell:2017hsw and the massive heavy flavor productions Li:2016tvb ; Li:2018tsq . However, at NNLO, the most general calculations can only be covered when one considers the $N$-jettiness soft function with four external hard legs, which is missing in the original calculations. In this manuscript, we present such calculation by extending the previously developed numerical method and show sufficient calculating details. We noticed a recent conference proceeding Bell:2018mkk on the $2$-jettiness soft function using pySecDec. Our work uses a different and independent approach which can serve as a cross check of the calculations.

The rest of this work is organized as follows. We discuss the numerical setups and the phase space parameterization in detail in Section II. In Section III, we present our numerical results for the $2$-jettiness observable in hardon-hadron collisions. Last we conclude in Section IV

II Computational Setup

In this section, we summarize the setups used in our calculation. We follow the general numerical framework proposed in Boughezal:2015eha , and show its capability in dealing with $N$-jettiness with four hard reference vectors in $pp\to jj$, $ep\to 3j$ and $0\to 4j$. We note that the calculations with four reference vectors go through all possible configurations in calculating any $N$-jettiness soft functions at NNLO. We clarify some of the computational details of this framework.

For simplicity, we stick to the hadronic $N$-jettiness definition which reads

[TABLE]

where $n_{i}$’s are the $N$ light-like hard reference vectors along the initial state hadronic beams or the final state jet axes and $q_{i}$’s are the four momentum of the soft radiations. We note that here all the hard vectors $n_{i}$’s are physical objects and are defined in $4$-dimension while the soft momentum $q_{i}$’s can potentially be un-resolvable and are therefore $d$-dimensional quantities, with $d=4-2\epsilon$ in the dimensional regularization. The soft function with $n$ real soft emissions then can be calculated by performing the phase space integration

[TABLE]

where we have introduced the $N$-jettiness measurement function

[TABLE]

and ${\cal S}$ is the soft current matrix element for the soft radiations which, up to NNLO, has been known for a long time Catani:1999ss ; Catani:2000pi . For $N$-jettiness with $3$ reference vectors , to NNLO, the soft current can be expressed in terms of a sum of dipole contributions, while starting from four reference vectors, like the case we study here, new triple-pole structure involving 3 eikonal lines arises which is induced by the one-loop soft current. The triple-pole contribution, which will be calculated later in this work, is the only remaining piece for calculating the $N$-jettiness soft function with arbitrary external eikonal lines at ${\cal O}(\alpha_{s}^{2})$. For completeness, we summarize the explicit form of the soft current ${\cal S}$ to ${\cal O}(\alpha_{s}^{2})$ in the Appendix V.1.

In order to compute the soft function numerically, the key step is to find a suitable momentum parameterization to allow us to write the phase space integration in the form of

[TABLE]

which is ready for the Laurent expansion so as for us to identify all the $\epsilon$-poles and evaluate the coefficients of the poles numerically. We thus now turn to the discussion of the phase space parameterization for single and double real emissions with four reference directions, separately. We note that though we are focusing on the four reference vector case, the parameterization is general enough to be applicable directly to any $N$-jettiness computation at the NNLO. In all cases below, we utilize the Sudakov decomposition to write an arbitrary light-like four vector $k^{\mu}$ with respect to two of the reference vectors $n_{i}$ and $n_{j}$ as

[TABLE]

with $k^{+}=n_{i}\cdot k$, $k^{-}=n_{j}\cdot k$ and $k_{\perp}\cdot n_{i}=k_{\perp}\cdot n_{j}=0$ and $|k_{\perp}|^{2}=\frac{2}{n_{ij}}k^{+}k^{-}$. Here we have introduced the notation $n_{ij}\equiv n_{i}\cdot n_{j}$.

II.1 single real emission phase space parameterization

With 4 reference vectors, $n_{i}$, $n_{j}$, $n_{k}$ and $n_{l}$, at hand, for one real emission with momentum $q_{1}^{\mu}$ for computing the NLO soft function or the NNLO real-virtual contribution, we can parameterize all the involved vectors with respect to $n_{i}$ and $n_{j}$, which are

[TABLE]

where we have used the freedom to choose the azimuthal angle of the vector $n_{k}$ and the additional $\epsilon$-angles other than $\alpha_{1}$ to be [math]. Here $\alpha_{1}\,,\phi_{l}\,,\phi_{1}\in[0,\pi]$ with

[TABLE]

Here the non-zero $\phi_{l}$ makes $n_{k}$ and $n_{l}$ span the azimuthal plain and therefore $\alpha_{1}$ is necessary for parameterizing the $d$-dimensional momentum $q_{1}$, since $q_{1}$ in general lies outside the pain $n_{k}$ and $n_{l}$ belong to. We note that since the reference vectors $n_{i}$’s are physical, for arbitrary numbers of the reference vectors, their “$\perp$”-components will all be lying in one azimuthal plain, and thus Eq. (II.1) is the most general parameterization for computing the $N$-jettiness function with arbitrary $N$. For the processes similar to what we are considering here, the transverse components $n_{k}^{\perp}$ and $n_{l}^{\perp}$ have to align with each other ($\phi_{l}=0$ or $\pi$) due to the momentum conservation, therefore one can simplify the parameterization by setting $\alpha_{1}=0$ which reduces to the parameterization used in Boughezal:2015eha ; Campbell:2017hsw for computing $1$-jettiness.

Although in general $n_{i}$ and $n_{j}$ can be chosen arbitrarily for performing the Sudakov decomposition, in order to correctly isolate the $\epsilon$-poles, it is more useful to choose $n_{i}$ so that $n_{i}\cdot q$ contributes to $\tau$ as set by $\Theta_{N,i}$ while $n_{j}\cdot q$ appears as one of the singular poles in the denominator of the soft current matrix ${\cal S}$. With this parameterization, for non-zero $\alpha$, the phase space integration for single emission can be modified to

[TABLE]

where we have made the variable changes

[TABLE]

It can also be shown that for $\alpha_{1}=0$, one can further modify the phase space integration to make $\phi_{1}\in[0,2\pi]$ when taking into account the properties of the soft current and the $N$-jettiness observable, as has been done in Boughezal:2015eha .

Eq. (9) is our starting point for calculating the single emission contributions to the soft function numerically and we list in the Appendix V.2 all the final forms of the phase space integrations suitable for direct numerical evaluation of the one emission contributions.

II.2 double real emission phase space parameterization

Now we sketch the parameterization for the double real emission with momenta $q_{1}$ and $q_{2}$ at NNLO. We want to add on top of Eq. (II.1) the momentum $q_{2}$ parameterization, meanwhile require that the dot product $q_{1}^{\perp}\cdot q_{2}^{\perp}$ only relies on the azimuthal angle between these two vectors, so as to use the non-linear transformation Boughezal:2015eha in a later stage for the purpose of pole isolation. To find the particular parameterization of $q_{2}$, we first rotate Eq. (II.1) around the $n_{k^{\perp}}$ axis by $\alpha_{1}$ to eliminate the $\alpha_{1}$ dependence in $q_{1}^{\mu}$, then we rotate in the azimuthal plain by $\phi_{1}$ to align the $x$-axis with $q_{1}^{\perp}$ to further remove the dependence of $\phi_{1}$ in the parameterization of $q_{1}^{\mu}$. These rotational operations lead to the momentum representation

[TABLE]

where we introduced the notation $c^{\alpha_{1}}_{\phi_{l}-\phi_{1}}\equiv c_{\phi_{l}}c_{\phi_{1}}+s_{\phi_{l}}s_{\phi_{1}}c_{\alpha_{1}}$. It is now straightforward to see that to satisfy our requirements, $q_{2}^{\mu}$ can be parameterized (in this frame after rotation) as

[TABLE]

so that $q^{\perp}_{1}\cdot q^{\perp}_{2}=-|q^{\perp}_{1}||q_{2}^{\perp}|\,c_{\phi_{2}}$. The $\beta$ angle here is necessary to guarantee that the $\epsilon$-components of $q_{1}$ and $q_{2}$ are not always aligned with each other. Due to the Lorenz invariance of the individual emission phase space measure, the rotational operations used to derive the parameterization of $q_{2}$ do not introduce any complexity to the Jacobian of the phase space integration. We also note that the parameterization derived here is the most general for computing the double real contributions to any $N$-jettiness soft functions. As for the processes we are interested in in this work, using the same momentum conservation for the single real emission case, one can set $\alpha_{1}$ and $\beta$ to [math] to simplify the phase space integration.

Within this parameterization, the dot product $q_{1}\cdot q_{2}$ takes the form

[TABLE]

which generates a linear singularity in the double real soft current if $q_{1}$ be collinear to $q_{2}$ is possible, which prevents us from extracting the $\epsilon$-poles using the Laurent expansion in Eq. (5). To extract the poles, the angle $c_{\phi_{2}}$ is further parameterized as $c_{\phi_{2}}=1-2\eta$, with the non-linear parameterization

[TABLE]

if the soft momentum $q_{1}$ and $q_{2}$ can be collinear to each other, and

[TABLE]

if $q_{1}$ and $q_{2}$ can no way be collinear as constrained by the jettiness measurement, for instance when $\tau=n_{i}\cdot q_{1}+n_{j}\cdot q_{2}$.

Last we comment on the choice of $n_{i}$ and $n_{j}$ in the Sudakov decomposition. The same logic follows the single emission case. The reference vectors $n_{i}$ and $n_{j}$ are chosen if $n_{i}\cdot q_{1}$ and $n_{j}\cdot q_{2}$ contribute to $\tau$, i.e. $\tau=n_{i}\cdot q_{1}+n_{j}\cdot q_{2}$. Otherwise, if both $n_{i}\cdot q_{1}$ and $n_{i}\cdot q_{2}$ contribute to the jettiness observable, then $n_{i}$ will be naturally picked as one reference axis for the decomposition, and for convenience, $n_{j}$ is chosen in the way that the double real soft current is singular as $n_{j}\cdot q_{1}\to 0$ and/or $n_{j}\cdot q_{2}\to 0$. Additional variable changes are needed to map all variables onto the regime $[0,1]$ suitable for doing the Laurent expansion. The mapping is standard in the sector decomposition calculations and has been detailed in Campbell:2017hsw ; Li:2018tsq . After the mapping, in the cases where one encounters the fractional power of the form $s^{-\frac{1}{2}+\epsilon}$ with $s\in[0,1]$, a further variable transform $s=\sin^{2}\left(\frac{\pi}{2}x\right)$, with $x\in[0,1]$, is then used.

III Numerical Results

In this section, we present our main results. We highlight our results using the $2$-jettiness soft function to ${\cal O}(\alpha_{s}^{2})$ in $pp\to jj$. Other processes such as $e^{+}e^{-}\to 4j$ can be easily obtained by changing properly the values of $\lambda_{ij}$ in Eq. (21), without doing any additional integrations. The $2$-jettiness soft function depends on four hard directions, which we denote as $n_{1}$, $n_{2}$, $n_{3}$ and $n_{4}$. We align $n_{1}$ and $n_{2}$ with the incoming beam axes in the $\pm z$-directions, and let $n_{3}$ and $n_{4}$ lie along the outgoing jet directions. For the purpose of presenting, though not at all required in our calculation, we let the incoming partons carry exactly the same energy in the collision and show the $2$-jettiness soft function as a function of $n_{13}$ and we present the results in the Laplacian space $\int_{0}^{\infty}\mathrm{d}\tau\,e^{-\tau\,z}\sigma(\tau)$.

We start from the NLO soft function. At NLO, the renormalized $2$-jettiness soft function can be written as a sum of dipole terms

[TABLE]

where $C_{ij}^{(1,i)}$ are the coefficients of the $\epsilon^{i}$-poles returned by our calculation using the method for single emission as sketched before. The coefficients of the logarithmic terms $\tilde{L}^{2}$ and $\tilde{L}$ can be predicted by solving the renormalization group equation (RGE) of the soft function and therefore serve as checks of our numerical approach.

In fig. 1, we show the NLO coefficients of the $\tilde{L}^{n}$ ($n=2,1,0$) terms for dipole $S_{12}^{(1)}$ and $S_{13}^{(1)}$ as functions of $n_{13}$. These are the only independent terms in the special kinematics we have chosen.

In both cases, we have normalized the results to $\frac{\alpha_{s}}{4\pi}$. The red, blue and black dots represent the coefficients for $\tilde{L}^{2}$, $\tilde{L}^{1}$ and $\tilde{L}^{0}$ from our numerical set-ups, respectively. The solid lines are obtained by solving the RGE analytically. We see perfect agreements between the analytic results and our numerical predictions, which implies the validity of the computational framework at this order.

At NNLO, other than the abelian terms which can be obtained by exponentiate the NLO results, the soft function receives non-abelian contributions from the $q{\bar{q}}$-emission, the $gg$-emission, the real-virtual corrections and the renormalization. The real-virtual corrections can be further broke down to the dipole and the triple-pole contributions, and the rest are of the dipole forms. The summation of all the dipole contributions can be organized by the color factors $C_{A}$ and $N_{F}T_{R}$, while the triple-pole real-virtual corrections can be massaged into one single color term $f_{abc}T^{a}_{1}T^{b}_{2}T^{c}_{3}$ with the help of the color charge conservation $\sum T_{i}=0$. The final results of the NNLO soft function can therefore be written as

[TABLE]

where all the NNLO dipoles $S_{ij}^{(2)}$ and triple-pole $S_{123}^{(2)}$ in the non-abelian terms can be written as a series in terms of the logarithms $\tilde{L}$. And once again, all the coefficients of those logarithmic terms $\tilde{L}^{n}$ with $n=3,2,1$ can be predicted by RGE with the knowledge of the full NLO results obtained before.

The coefficients of the $N_{F}T_{R}$ terms for the dipoles $S_{12}^{(2)}$ and $S_{13}^{(2)}$ are shown in fig. 2, which include the contributions from the $q{\bar{q}}$ emission and the $N_{F}$ term in the renormalization.

Perfect agreements are observed between the numerical calculations of the $\tilde{L}^{n}$ coefficients for $n=3,2,1$ as shown in colored dots and the RGE predictions in solid lines, which validates not only the NNLO computations of the $N_{F}$ contributions but also the correctness of the non-logarithmic term in the NLO calculations. The NNLO non-logarithmic prediction (in black dots) displayed here is not presented in the original NNLO $N$-jettiness soft paper Boughezal:2015eha and is one of our main results of this work. We note that the NNLO results here and below are also presented in a conference proceeding Bell:2018mkk in a different way with a different approach.

Same situation happens to the $C_{A}$ term, as displayed in fig. 3, in which we also compared $\tilde{L}^{n}$ coefficients from the numerical calculation against the results from the RGE to find no difference. The non-logarithmic coefficients are represented by the black dots.

Last we give the prediction for the triple-pole term $S_{123}^{(2)}$, which arises firstly for four hard reference directions. Just like the dipole contributions, the logarithms in the triple-pole are still predictable from the RGE. The analytic predictions are once again found to agree with our numerical results, as clearly seen in the left panel of fig. 4. Other than that, we also compare the $\epsilon$-pole terms without the $\tilde{L}$’s from our numerical calculations with the ones predicted by the RGEs and find satisfactory agreements. The result for the non-logarithmic term (left panel in fig. 4) as well as the term entirely unpredictable from RGE (right panel, normalized further to $-2\pi$) are indicated in black dots, which finalizes all the NNLO contributions to the $2$-jettiness for the $pp\to jj$ process.

IV Conclusions

In this work, we present the calculation of the $N$-jettiness soft function with four hard reference directions using the computational framework proposed in Boughezal:2015eha using the approach of sector decomposition. We generalize the parameterizations used for the $1$-jettiness calculation to the arbitrary $N$-jettiness case. We specifically calculated the NNLO $N$-jettiness soft function with four external hard legs, where the last new piece at this order, the triple-pole configuration, first arises. We managed to isolate all the singularities from the soft function integrands and reduce the computation to a set of numerical integrals which are ready to evaluate. We check the numerically computed logarithmic terms in the soft function against the predictions from the RGEs and found perfect agreements. The non-logarithmic terms which can not be obtained through the known RGEs are calculated directly in this work. We expect the result obtained in this manuscript to find its quick applications in both the fixed NNLO calculations and achieving the parton shower matching Alioli:2013hqa for the relevant processes.

V Acknowledgments

This work is supported by the National Natural Science Foundation of China under Grant No. 11775023 and the Fundamental Research Funds for the Central Universities. S.J. is supported by the BNU JingShi Undergraduate Fellowship.

Appendix

V.1 soft current

For completeness, we list all the soft current matrix elements here, which were first derived in the seminal works Catani:1999ss and Catani:2000pi by Catani and Grazzini.

The NLO soft current with one soft emission $q$ is given by

[TABLE]

with the NLO dipole given by

[TABLE]

The NNLO real-virtual correction can be written as the sum of the dipole ${\cal S}^{(2)}_{RV,dipole}$ and triple-pole ${\cal S}^{(2)}_{RV,tri.}$ contributions, with the dipole contribution be

[TABLE]

and the triple-pole term

[TABLE]

where $\lambda_{ij}=1$ if $i$ and $j$ are both incoming or outgoing, otherwise $\lambda_{ij}=0$.

The NNLO real-real emission is given by

[TABLE]

from the $q{\bar{q}}$ emission. To get this form we have used the color charge conservation $\sum_{i}T_{i}=0$. Here

[TABLE]

and

[TABLE]

where $q_{1}$ and $q_{2}$ are the momentum carried by the real emissions. We denote the term proportional to $2\,n_{ij}\,q_{1}\cdot q_{2}$ as ${\cal I}_{ij}^{\rm II}$ and the rest ${\cal I}_{ij}^{\rm I}$.

And the $gg$ emission gives the non-abelian contribution

[TABLE]

with

[TABLE]

where

[TABLE]

and the current in the strongly ordering limit

[TABLE]

The abelian piece for the $gg$ double real emission is

[TABLE]

V.2 numerical integrations for single emission

Here we list all final numerical integrations for evaluate single emission contributions. We note again that though we only show the $N$-jettiness with four hard directions $n_{i}$, $n_{j}$, $n_{k}$ and $n_{l}$, the parameterization here presented here is general enough for any $N$-jettiness soft functions. To extend to the case with more reference vectors, one simply insert more $\theta$ functions originated from the $N$-jettiness measurements. We use the notation $q||n_{i}$ to denote the separation of $q$ and $n_{i}$ is the smallest and thus $\tau=n_{i}\cdot q$.

V.2.1 final forms for numerical integration at NLO

At NLO, the contribution can be written as a sum of dipoles. We assume that the soft momentum $q$ is emitted from dipole $ij$. For the integrations shown below, we have extracted out an overall factor $-16\pi^{2}\left(\frac{\alpha_{s}}{2\pi}\right)\left(\frac{e^{\gamma_{E}}}{4\pi}\right)^{\epsilon}T_{i}\cdot T_{j}$.

• •

case 1: $q_{1}\parallel n_{i}$.

For $q_{1}\parallel n_{i}$, the final form suitable for numerical integration is

[TABLE]

where we have defined

[TABLE]

and

[TABLE]

As for the case in which $q_{1}\parallel n_{j}$, the result can be obtained by switching $i$ and $j$ in Eq. (30).

• •

case 2: $q_{1}\parallel n_{k}$.

While for $q_{1}\parallel n_{k}$, we choose $n_{k}$ and $n_{i}$ to be the reference vectors, to have $q_{1}^{+}=n_{k}\cdot q_{1},q_{1}^{-}=n_{i}\cdot q_{1}$. Following the same variable changes shown in Eq. (10), we get the final form for the numerical evaluation

[TABLE]

where

[TABLE]

Again, we can switch $k$ and $l$ to obtain the $q_{1}\parallel n_{l}$ contribution.

V.2.2 final forms for $2$-parton correlated real-virtual contribution

For the $2$-parton correlated real-virtual contribution, we have for our numerical evaluation

• •

case 1: $q_{1}\parallel n_{i}$.

[TABLE]

• •

case 2: $q_{1}\parallel n_{k}$.

[TABLE]

Here we have normalized to

[TABLE]

The $q||n_{j}$ and $q||n_{l}$ cases can be easily obtained by switching properly the indices.

V.2.3 final forms for 3-parton correlated real-virtual contribution

Now we turn to the triple-pole case in the real-virtual correction, in which $3$-parton correlated emission contributes. This configuration first arises in this case with four external legs. We assume the soft emission is from triple-pole $ijk$. And we normalize our results to

[TABLE]

• •

case 1: $q_{1}\parallel n_{i}$.

[TABLE]

where $F^{ij}_{i}$ has been defined before.

• •

case 2: $q_{1}\parallel n_{j}$

We let $q_{1}^{+}=n_{j}\cdot q_{1}$, $q_{1}^{-}=n_{i}\cdot q_{1}$ to find

[TABLE]

• •

case 3: $q_{1}\parallel n_{k}$

We let $q_{1}^{+}=n_{k}\cdot q_{1}$ and $q_{1}^{-}=n_{l}\cdot q_{1}$ to get

[TABLE]

Here we redefine $F^{ij}_{k}(x_{3})$ as

[TABLE]

• •

case 4: $q_{1}\parallel n_{l}$

We let $q_{1}^{+}=n_{l}\cdot q_{1}$ and $q_{1}^{-}=n_{k}\cdot q_{1}$ to have

[TABLE]

where

[TABLE]