Next-to-leading order QCD corrections to di-photon production in association with up to three jets at the Large Hadron Collider

We present the computation of next-to-leading order (NLO) QCD corrections to di-photon production in association with two or three hard jets in pp collisions at a center-of-mass energy of 8 TeV. The inclusion of NLO corrections is shown to substantially reduce the theoretical uncertainties estimated from scale variations on total cross section predictions. We study a range of differential distributions relevant for phenomenological studies of photon pair production in association with jets at the LHC. Using an efficient computational set-up we performed a detailed study of uncertainties due to parton distribution functions. The computation of the virtual corrections is performed using new features of the C++ library NJET.


Introduction
There has been a lot of recent interest in the study of pp → γγ + jets processes as an important background to pp → H → γγ which is one of the cleanest decay channels for studies of the Higgs properties [1,2]. These processes also have importance outside the realm of Higgs physics testing our ability to model isolated photon radiation in association with strong interactions, see for example Refs. [3,4] for recent experimental studies of the closely related process of photon production in association with hard jets. From a theoretical point of view di-photon production is under good control with corrections known up to NNLO in QCD [5]. NLO QCD corrections to pp → γγ + j have been available for some time [6] and have recently been re-explored [7] to investigate the impact of using different photon isolation criteria comparing the smooth cone isolation [8] with the fixed cone isolation favoured in experimental studies. The computation of pp → γγ+2j including NLO QCD corrections has been completed quite recently [9,10].
In this paper we extend the available range of predictions for pp → γγ + jets to include up to three hard jets with NLO accuracy in QCD. The results show a significant reduction in the uncertainty of the theoretical predictions and have noticeable corrections to the shapes of the distributions when going from LO to NLO. As well as providing new phenomenological studies relevant for the current measurements at ATLAS and CMS we also present new developments to the NJet C++ code enabling more efficient computations of high multiplicity processes at NLO.
This article is organized as follows: we begin by outlining our computational set-up in Section 2, where for completeness we review the well known decomposition of nextto-leading order differential cross sections and describe the interface of NJet with the Sherpa Monte-Carlo, which we used for the computation of the unresolved real radiation contributions and phase-space integration. In Section 3 we provide results for the LHC at centre-of-mass energy of 8 TeV for both pp → γγ + 2j and pp → γγ + 3j. We present differential distributions for some important observables, particularly those used in Higgs productions studies with vector-boson fusion phase space in the case of pp → γγ + 2j. We present a study of the dependence on the renormalization scale of the NLO predictions and investigate the uncertainty due to the choice of Parton Distribution Functions (PDFs) on total cross sections and distributions. In Section 4 we present our conclusions.

Computational set-up
The computation is performed in the five-flavour scheme with massless b-quark included in the initial state. The basic partonic sub-processes considered are: for pp → γγ+3j, from which all relevant channels can be derived using crossing symmetries. Channels with like-flavour fermions are obtained from the above using the appropriate (anti-)symmetrization relations. We do not include loop-induced and formally higher order sub-processes 0 → γγ + 4g and 0 → γγ + 5g. We can schematically write down the NLO partonic cross section as a sum of four finite contributions which can be integrated separately over their respective phase spaces, where dσ B n denotes the leading order contribution, dσ I N contains the integrated dipole subtraction terms, including factorization contributions from initial state singularities, dσ V n the one-loop virtual corrections and dσ RS n+1 the infra-red finite contributions from realradiation and dipole subtraction terms. This hard partonic cross section is then convoluted with the parton distribution functions to obtain the cross sections for hadronic collisions.
The computation of the Born and real-emission matrix elements is performed using the colour dressed recursive Berends-Giele formulation [26] implemented in the Comix amplitude generator [27]. The subtraction of infra-red singularities is performed using the Catani-Seymour [28] dipole method. The evaluation of these contributions is performed using the Sherpa [29,30] package, which is also used for the determination and organization of the partonic subprocesses and the integration over the phase space. The one-loop virtual amplitudes are provided using the automated generalized unitarity framework implemented in the latest version of the NJet library 1 . NJet 2.0 is an updated code based on NJet [31] and NGluon [32]. The primitive kinematic objects are constructed using a numerical generalized unitarity algorithm [11-13, 20, 33-36] with Berends-Giele recursion relations for the tree-level input [26]. The extended code can compute arbitrary primitive amplitudes for vector bosons (W, Z and γ) and massless QCD partons. Full colour sums are implemented for pure QCD with up to five jets, vector bosons with up to five jets and di-photon production with up to four jets. We interface this code to the Sherpa Monte-Carlo via the updated Binoth Les Houches Accord [37,38] to obtain the virtual events. In addition to the standard internal checks we have managed to check individual phase-space points for the processes used in this paper against those obtained with GoSam [23] and MadLoop [21]. For both processes we have neglected the small effect of top quark loops in the virtual amplitudes. In the case of pp → γγ + 3j we also neglect the contribution from vector loops where the photons couple directly to a virtual fermion loop. These corrections have been included in pp → γγ + 2j contributing less than 0.5% of the total cross section, therefore they are expected to be negligible for pp → γγ + 3j.
For the current calculation we made use of some new features to optimise the computation time needed for the virtual corrections. Firstly we use a C++ library Vc [39] to utilize vector capabilities of modern CPUs and gain a factor of ∼ 2 in the computation speed. We also separate leading and sub-leading contributions in colour such that the simpler, dominant contributions can be sampled more often than the sub-leading terms. The definition of our leading terms include all multi-quark processes in the large N c limit and processes with two or more gluons in the final state using the de-symmetrized form of the colour sum that efficiently exploits the Bose symmetry of the phase space [31,40]. The de-symmetrized sums give full colour information and are faster than leading colour when we have many final state gluons. In Figure 1 we show virtual corrections to the 3rd leading jet transverse momentum in pp → γγ + 3j. The corrections from the sub-leading part are around 10% on average with a slight rise to around 20% at large p T . In the case of pp → γγ + 3j the virtual cross sections are about 1/3 of the size of the total cross section. For the current implementation of pp → γγ + 3j, the leading virtual events are generated approximately 7 times faster than the sub-leading events.
In order to maximise the phenomenological predictions that we can extract from processes with such complicated final states, we make use of the ROOT Ntuple format provided by Sherpa [41]. During the course of event generation, additional weights from the poles in the loop process and from the subtraction terms are stored along with a full list of kinematic variables and couplings. This allows the events to be re-weighted to different scale and PDF choices during analysis. Even within this approach a full study of PDF uncertainties and scale variations can be computationally intensive. During this work we have developed an interface to the APPLgrid library [42] to allow for specific observables to be efficiently re-weighted changing the scales and the specific PDF set used.

Numerical results
All the results presented in this section are for pp collisions with a centre-of-mass energy of 8 TeV. We consider the following kinematic cuts on the external momenta, which are inspired by typical experimental cuts used in the analyses at LHC where the photon transverse momenta have been ordered by size. The jets are defined using the anti-k T algorithm [43] with cone size R = 0.5 as implemented in FastJet [44].
Photons are selected using the Frixione smooth cone isolation criterion [8]. A photon is considered isolated if the total hadronic energy inside all cones of radius r γ < R with = 0.05, R = 0.4 and n = 1. We use the NLO CT10 PDF set [45] for our central predictions with the strong coupling running from α s (M Z ) = 0.118, and the electromagnetic coupling fixed at α = 1/137.036. In particular we use the same (NLO) PDF set and definition of the strong coupling constant both for LO and NLO predictions. Using a NLO PDF set for the LO computation includes higher order terms that go beyond a fixed order prediction, nevertheless such a set-up allows us to separate NLO effects coming from the running of the strong coupling and from PDFs and to highlight the impact of corrections coming from the NLO matrix elements. We choose a dynamical value for the factorization and renormalization scales which are kept equal, µ R = µ F , when performing scale variations. We have investigated the dependence on a number of different functional forms which we will denote as: where m γγ = p 2 T,γ 1 + p 2 T,γ 2 . The quantities H T and Σ 2 are constructed after the clustering of final state partons into jets. Notice that partonic and jet-level scales will only differ at NLO, where the additional unresolved radiation enters in the clustering algorithm.

Results for pp → γγ + 2j
We first consider the production of a photon pair in association with two jets. In this case we compare our predictions with the recent results of reference [9] and present additional studies of PDF variations and dynamical scale choices. For the latter we rely on the possibility to substantially reduce the computational cost by the use of our APPLgrid set-up.
In Figure 2 we show the dependence of the total inclusive cross section upon variation of the renormalization and factorization scales with µ R = µ F . We consider the five different dynamical scales defined in Eq. (3.6). Though the scale choices are closely related to each other we find significant deviations for the value of the total cross section. Nevertheless the NLO predictions show a significant reduction in the dependence on the scale variation compared to the LO ones. Taking the envelope of all the scales considered we see that the LO predictions vary in the interval 1.64−3.04 whereas NLO ones lie within 2.46−3.58 when the scales are varied over the range x ∈ [0.5, 2] around the central choice. This represents a reduction in the scale variation uncertainty from ∼ 30% at LO to ∼ 20% at NLO. NJet + Sherpa pp → γγ + 2j @ 8 TeV We notice that choosing µ R = Σ 2 as a scale leads to significantly larger predictions for the total cross section than with the other scales considered here. Moreover the scale variation profile closely resembles the one of a leading order prediction. We therefore consider this choice of scale disfavoured with respect to the others, which provide results that are in better agreement with each other. In general, we find that the H T -based scales lead to a broader profile of scale variations and on average favour harder events than the Σ-based scales. On the other hand we find the shapes of the distributions to be quite stable with respect to the scale choice.
The results for the total cross section and distributions using µ R = √ Σ 2 /2 are in good agreement with those obtained previously by Gehrmann, Greiner and Heinrich [9]. In the following we opt for the scale H T for computing our predictions for the total cross section and the differential distributions presented here. This scale has been widely used in studies of W + jets (see for example [18]).
The values for the total cross sections at both LO and NLO computed our default choices for scale, PDF set and physical parameters are found to be: where the sub-scripts(super-scripts) show the maximum deviation from the central value in the negative(positive) direction over the range x ∈ [0.5, 2] for µ R = x H T /2. Monte-Carlo integration errors are shown in brackets. Figure 3 shows the differential distributions for the ordered jet transverse momenta. The results for the scale H T /2 are consistent with those obtained at the scale √ Σ 2 with a significant reduction in scale variation from around 20% at LO to ∼ 10% at NLO in both cases. The K-factor is fairly constant at around 1.1 for p T higher than 200 GeV rising to 1.4 as we approach the p T cut. This larger K-factor in the low p T region could be the indication of the presence of large logarithms beyond fixed order NLO. The di-photon invariant mass m γγ and the di-photon rapidity distributions η γγ are shown in Figure 4. They receive slightly larger NLO corrections with respect to the jet transverse momenta with the K-factor for the former ranging from 1.2 in the large invariant mass region to 1.7 for lower invariant masses, while for the latter NLO correction induce a roughly flat K-factor of 1.3. Figure 5 shows four distributions of angular quantities that are usually used in analyses of Higgs production in vector boson fusion (VBF), where additional cuts are imposed in order to reduce QCD backgrounds in pp → H(→ γγ) + 2j studies.
Owing to the increased phase-space available to the real radiation at NLO we see large deviations from the shapes of the leading order distributions. These features have been pointed out for the jet-pair azimuthal angle ∆φ j 1 j 2 and for the separation of the leading-photon/leading jet, in Ref. [9]. We also see increasing deviations for large rapidities of the jet pair η j 1 j 2 and even more so in the relative rapidity of the diphoton-dijet system, y * γγjj = y γγ − (y j 1 + y j 2 )/2. (3.9) Using the APPLgrid set-up described in the previous section, we have also performed a study of PDF uncertainties on pp → γγ + 2j, concentrating on the total cross section and the invariant mass distribution of the photon pair. In Table 1 we show the central value and PDF uncertainties for the total cross section evaluated at the central scale (µ R = H T /2) using four different NLO PDF sets: CT10 [45], MSTW2008 [46], ABM11 [47] and NNPDF2.3 [48]. All PDF sets are compared using the same value of α s (M Z ) = 0.118 in order to disentangle PDF and strong coupling constant uncertainties. PDF uncertainties are considerably smaller than the theoretical uncertainty estimated from scale variations and range from 1% for MSTW and NNPDF to 3.5% for ABM. We note, however that the ABM11 uncertainty includes errors associated with α s variations. In general, we find that central predictions from different PDF sets differ by amounts which are larger than the nominal PDF uncertainty of each set.
In Figure 6 we show the distributions for the invariant mass of the photon pair (m γγ ) computed at the scale µ R = H T /2 and the four PDF sets considered before at α s (M Z ) = 0.118. The upper (log scale) plot shows the absolute predictions for all sets, which are indeed extremely close to each other. In the lower plot, where we show the ratio of each  set to the central value computed using CT10, we do however notice discrepancies in the predictions of the order of a few percent. These differences are especially pronounced at high m γγ where the spread between them reaches ∼ 6%, a value larger than PDF uncertainties coming from individual sets (represented by the shaded regions). We also note that, while the differences between CT10, MSTW08 and NNPDF2.3 are mostly due to an overall normalization factor, the ABM11 predictions stand out from the others also in the shape being higher for low values of m γγ and dropping below the other predictions for higher values of the invariant mass.    Figure 6: The m γγ distributions for pp → γγ + 2j for the four PDF sets described in the text at α s M Z = 0.118. The lower plot shows the ratio of each set to CT10 with the shaded region representing the PDF uncertainty.

Results for pp → γγ + 3j
We now consider the production of a photon pair in association with three jets. As in the previous section we studied the dependence of the total cross section upon variation of renormalization and factorization scales with the choices of dynamical scales defined in Eq. (3.6). The results in Figure 7 show reasonable differences between quantities based on jets versus quantities based on partons. Overall we find a significant improvement in the uncertainty estimated from scale variations when going from LO to NLO. The envelope of predictions from all scale choices varied over the range x ∈ [0.5, 2] is around 0.67 − 0.99 pb at NLO compared to 0.46 − 1.28 pb at LO. This represents a decrease in variation from ∼ 50% at LO to ∼ 20% at NLO. As in the two jet case the scales based on Σ 2 give generally larger predictions than those based on H T . Other than the overall normalization, we find that all scales give very similar predictions for shapes of the distributions. Comparing Figure 7 with Figure 2 we see that the peak in the NLO curve for Σ 2 has moved further to the right than the H T scales which may suggest that a range of x ∈ [1,4] would be more appropriate here. Since we would like to make predictions for jet ratios we need to have as consistent description of γγ + 3j and γγ + 2j as possible and therefore we prefer the H T scales. In the following we choose to adopt the central scale of H T /2 for the total rates and distributions, though theoretical uncertainties are likely underestimated by the simple scale variations following the discussion above. For the total cross sections at LO and NLO we find, NJet + Sherpa pp → γγ + 3j @ 8 TeV The distributions for the jets transverse momenta are shown in Figure 8. The K-factor is quite flat with a value between 1.0 and 1.2 except in the low p T where it rises to around 1.4 for the leading jet p T . Again, this may suggest the presence of large logarithms in the missing higher order contributions in this region. The distribution of the 3rd jet shown in Figure 9a has a noticeably flatter K-factor than the ones for the two leading jets.
The di-photon invariant mass distribution in Figure 9b receives significant corrections to the shape at NLO with the K-factors increasing from around 1.0 at low m γγ to 1.4 at large values of the photon pair invariant mass. In Figure 10 we show the leading jet/leading photon separation R γ 1 j 1 and the azimuthal separation of the two leading jets ∆φ j 1 j 2 . Both quantities receive large NLO corrections for small values of the observable, though notably not as large as the corresponding distributions in pp → γγ + 2j where there is a more substantial increase in the available phase-space at NLO.
The PDF analysis was performed with the same set-up as for pp → γγ+2j with the four PDF sets compared at the same value of α s (M Z ) = 0.118. Table 2 shows the results for the total cross section at the central scale µ R = H T /2. Again we see that PDF uncertainties for all sets are noticeably smaller than theoretical uncertainties estimated from scale variations  (b) Figure 9: Differential cross section as a function of p T of the 3rd leading jet and di-photon invariant mass in pp → γγ + 3j. and that central predictions from different sets differ by more than the nominal uncertainty obtained with each set. Figure 11 shows the comparison of the m γγ distribution for the different PDF sets. The deviation between ABM11 and NNPDF, which again give the lower and higher limit of the predictions from different sets, is somewhat larger than that seen in pp → γγ + 2j, though the results are consistent within the theoretical uncertainties determined by scale variations.

The three-to-two jet ratio
In the context of multi-jet production studies it is interesting to look at the ratio of pp → γγ + 3j over pp → γγ + 2j which we will denote as R 3/2 . Due to the cancellation of many uncertainties both theoretical and experimental, ratios such as this one are prime observables for the determination of physical parameters (like α s ).
In the case of the production of di-photon in association with jets, we find R 3/2 to be: where the numbers in brackets refer to Monte-Carlo errors. We have checked that all scale choices are in much better agreement for R 3/2 with the range of predictions lying within ∼ 8% of this value. In Figure 12 we show the differential ratio with the p T of the leading jet. The NLO corrections become more important for p T > 100 GeV and reaching about 15% at high p T .  Figure 11: The m γγ distributions for pp → γγ + 3j for the four PDF sets described in the text at α s M Z = 0.118. The lower plot contains the ratio of each set to CT10 with the shaded region representing the PDF uncertainties. (dσ n+1 /dp T )/(d σn/dp T ) NJet + Sherpa at 8 TeV LO γγ 3/2 NLO γγ 3/2 Figure 12: The ratio of of pp → γγ + 3j over pp → γγ + 2j as a function of leading jet p T .

Conclusions
In this paper we have presented a study of di-photon production in association with jets. The first calculation of the full NLO QCD corrections to the process pp → γγ + 3j is presented and discussed, together with a variety of results for pp → γγ + 2j. We find that the inclusion of NLO QCD corrections leads to a significant reduction of theoretical uncertainties both on total cross sections and distributions. We have studied distributions for a number of observables, of particular interest are those relevant for Higgs production in VBF analyses which are used when modelling γγ +jets as a background to pp → H +jets → γγ + jets. The present study is based on the use of the Frixione smooth cone isolation criterion to define the final state photons, which provides a theoretically clean way of suppressing the fragmentation component in direct photon production. It would be also interesting to consider alternative isolation criteria, like the fixed cone isolation, which require the inclusion of fragmentation functions, and were shown to have a significant effect in lower multiplicity processes [7]. Moreover the majority of experimental analyses involving direct photon production (both in association with jets or not) rely on the cone isolation for photon identification. Nonetheless, we hope that the results presented here will be of use in future experimental analyses and look forward to direct comparisons with the LHC data.