Predictions for diphoton production at the LHC through NNLO in QCD

In this paper we present a next-to-next-to-leading order (NNLO) calculation of the process pp → γγ that we have implemented into the parton level Monte Carlo code MCFM. We do not find agreement with the previous calculation of this process in the literature. In addition to the Oαs2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O}\left({\alpha}_s^2\right) $$\end{document} corrections present at NNLO, we include some effects arising at Oαs3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O}\left({\alpha}_s^3\right) $$\end{document}, namely those associated with gluon-initiated closed fermion loops. We investigate the role of this process in the context of studies of QCD at colliders and as a background for searches for new physics, paying particular attention to the diphoton invariant mass spectrum. We demonstrate that the NNLO QCD prediction for the shape of this spectrum agrees well with functional forms used in recent data-driven fits.


Introduction
The discovery of a light Higgs boson [1,2], which decays to two photons, has helped cement the diphoton process as one of the most interesting final states to study during the second run of the LHC (Run II). Experimental studies of prompt (γ) and diphoton (γγ) production at hadron colliders have been undertaken for several decades [3][4][5][6][7][8][9][10][11][12][13][14][15]. These studies are possible in part due to the high rate of production, but also because of the relative cleanliness of the experimental final state. As the energy available for collisions has increased, and as the amount of data collected has grown, so too has the region of diphoton invariant mass (m γγ ) that can be probed. At the LHC experimental data is now available up to scales of order 1 TeV, allowing for searches for new heavy resonances that may decay to photon pairs [16,17].
During Run II, the large data set will result in many measurements being performed at a level of detail that demands exquisite theoretical predictions. Therefore, in addition to the detailed experimental studies, the prompt and diphoton processes have received JHEP07(2016)148 considerable theoretical attention. The next-to-leading order (NLO) calculations embodied in the Jetphox [18] and Diphox [19] Monte Carlo codes have been extensively utilized in the experimental literature. In addition to the NLO calculation of diphoton production, gg → γγ contributions that are formally higher-order, but phenomenologically important, have also been computed [20,21]. However, existing 7 TeV analyses have already confirmed the inadequacies of NLO calculations when confronted with data [11,15]. Instead, much better agreement is found with the recently-completed next-to-next-to Leading Order (NNLO) calculation [22], that naturally subsumes the first gg → γγ contributions.
This calculation was made possible through the application of the Q T -subtraction procedure [23]. This procedure makes use of the known factorization properties at small transverse momenta of the diphoton system to efficiently handle complications arising from infrared singularities. Although a variety of other methods for regularizing and combining infrared singularities have been devised [24][25][26], and used to provide a suite of new predictions for 2 → 2 hadron collider processes [27][28][29][30][31][32], the relative simplicity of the Q Tsubtraction method is highly appealing [33][34][35][36][37][38][39]. The Q T -subtraction method generates a counter-term that regularizes the singularity as Q T → 0 but is otherwise non-local; it also naturally lends itself to implementation as a slicing method ("Q T -slicing"). A promising new development is a generalization of the Q T -based methods, which were originally only applicable to color-neutral final states, to new methods [40][41][42] based on Soft Collinear Effective Field Theory (SCET) [43][44][45][46][47]. One of these methods [41,42], based on the Njettiness global event shape [48], can in principle be applied to arbitrary processes [41,[49][50][51][52]. In its implementation as a slicing method, the N -jettiness variable (τ ) is used to split the phase space into two regions. In the region where τ > τ cut at least one of the additional partons is resolved. Therefore the calculation contains only single unresolved limits and is amenable to calculation using standard NLO techniques. For the second region, where τ < τ cut , both partons can be simultaneously unresolved. In this region a factorization theorem from SCET [48] is used to approximate the cross section to the desired perturbative accuracy. This is a natural generalization of Q T -subtraction, where a similar reasoning applies when replacing τ with Q T and SCET factorization with one based on the Collins-Soper-Sterman formalism [53].
The aim of this paper is to present a new NNLO calculation of pp → γγ using the N -jettiness slicing approach and compare it with the existing calculation of ref. [22]. Given its importance for Run II phenomenology an independent calculation is crucial. In fact we will find that we cannot reproduce the results of the literature and we believe that existing results for this process are inaccurate. This underlines the need for multiple independent calculations of processes such as this one that are of great importance for existing and future experimental analyses. We investigate the role of higher-order effects to the gg initiated closed loops of quarks, and combine this prediction with NNLO for the first time. We will also investigate the role of top quark loops at high invariant masses. Our calculation is implemented in MCFM [54][55][56] and will be released in a forthcoming version of the code.
We continue this paper by outlining the various component pieces of our calculation in section 2. In section 3 we compare our predictions to existing results from the literature and discuss the checks we performed on our calculation. In section 4 we turn our attention to JHEP07(2016)148 Figure 1. Representative Feynman diagrams for the calculation of pp → γγ at NNLO. From left to right these correspond to double virtual (calculated in ref. [57]), real-virtual and real-real corrections.
LHC phenomenology, comparing our predictions to data obtained by the CMS experiment at 7 TeV, and to the m γγ spectrum reported by ATLAS at 13 TeV. Finally, we draw our conclusions in section 5. Appendices A, B and C contain additional technical details of our calculation.

Calculation
In this section we present an overview of our calculation of diphoton production at NNLO and discuss the various contributions that are included in this paper. Before going into detail we introduce the following notation In this way ∆σ X represents the correction obtained from including the coefficient that first arises at order X in perturbation theory. We use this notation both inclusively (as written above) and for differential predictions.

Overview
We present representative Feynman diagrams for the various topologies that enter the calculation of the pp → γγ process at NNLO in figure 1. At this order in perturbation theory contributions arise from three distinct final states. The simplest is the one that also represents the Born contribution and corresponds to a 2 → 2 phase space. At NNLO this final state receives corrections from two-loop amplitudes interfered with the LO amplitude, and one-loop squared contributions. The 2 → 3 real-virtual phase space consists of treelevel and one-loop amplitudes for qqgγγ interfered with one another. Finally the largest phase space, representing a 2 → 4 process, is referred to as the double-real contribution and consists of two tree-level qqγγ + 2 parton amplitudes squared. The contributions discussed above have ultraviolet (UV) poles in the double-virtual and real-virtual phase spaces, which we renormalize in the MS scheme. Amplitudes for the double-virtual contribution can be found in ref. [57], for the real-virtual in ref. [58], and tree-level amplitudes for the real-real can be found in ref. [59]. After UV renormalization the individual component pieces of the calculation still contain singularities of infrared (IR) origin. These infrared poles must be regulated, made manifest, and combined across the different phase spaces in order to ensure that a sensible JHEP07(2016)148 prediction is obtained. As discussed in the introduction, we will use the N -jettiness slicing technique proposed in refs [41,42] for this task. This results in an above-cut contribution corresponding to the calculation of pp → γγj at NLO. The below-cut contribution requires 2-loop soft [60,61] and beam [62] functions, together with the process-dependent hard function. Various component pieces of this calculation, including explicit results for the hard function, are given in appendix A

gg initiated loops at LO and NLO
The NNLO calculation of γγ production represents the first order in perturbation theory that is sensitive to gg initial states. One class of gg configurations corresponds to real-real corrections, i.e. the gg → qqγγ matrix element that is related to the contribution shown in figure 1 (right) by crossing. These pieces are combined with contributions from the DGLAP evolution of the parton distribution functions in the real-virtual and double-virtual terms to ensure an IR-finite result. The second type of contribution is due to n F "box" loops, for which a representative Feynman diagram is shown in the top left corner of figure 2. This contribution has no tree-level analogue and is thus separately finite.
The box diagrams result in a sizeable cross section (≈ σ LO ), primarily due to the large gluon flux at LHC energies and the fact that this contribution sums over different quark flavors in the loop. In this section, we focus on n F = 5 light quark loops. Since this contribution is clearly important for phenomenology it is interesting to try to isolate and compute higher order corrections to it. We illustrate typical component pieces of these NLO corrections in the remaining diagrams in figure 2. They comprise two-loop ggγγ amplitudes, and one-loop gggγγ and gqqγγ amplitudes. A NLO calculation of gg → γγ including the two-loop and one-loop gggγγ amplitudes was presented in refs. [20,21]. An infrared-finite calculation can be obtained from the gg → γγ two-loop amplitudes and the gggγγ one-loop amplitudes, provided that a suitable modification to the quark PDFs is used (essentially using a LO evolution for the quark PDFs and a NLO evolution for the gluon PDFs). On the other hand if the qqgγγ amplitudes are included then the corresponding collinear singularity can be absorbed into the quark PDFs as normal at NLO, allowing JHEP07(2016)148 Figure 3. The ratio of the invariant mass distribution in gg → γγ computed using five light flavors and the effect of the top quark, to the calculation with n F = 5 alone. The dashed line shows the ratio of the result for n F = 6 to the one for n F = 5 that corresponds to eq. (2.2). for a fully consistent treatment. In the original calculation [20,21] (and the corresponding implementation in MCFM [55]) the first approach was taken. Here we will follow the second approach and include the qqgγγ amplitudes. Although formally an improvement, we find that the differences between the two approaches are negligible. Most of the required qqgγγ amplitudes can be found in ref. [58]. However, since that paper was concerned only with the NLO predictions for the γγj process, it did not include the one-loop amplitude that interferes with a vanishing tree-level term. In the calculation presented here this purelyrational amplitude is squared and therefore must be properly included. For completeness we present this missing amplitude in appendix C.
Since the NLO corrections to the gg initiated diagrams form a part of the N 3 LO cross section but do not represent a full calculation at that order, we define the additional cross section associated with them as ∆σ N3LO gg,n F . The subscript indicates that they are associated with gg initiated closed loops of quarks. Although by no means a complete O(α 3 s ) prediction, it is possible that the ∆σ N3LO gg,n F contribution forms a sizeable part of this correction. The impact of these terms will be discussed at length in section 4.

Impact of the top quark at high m γγ
The previous subsection outlined the calculation of gg loops for n F = 5 light quarks. While this is an excellent approximation for low invariant mass photon pairs, at higher energies this is no longer appropriate due to contributions from top quarks circulating in the loop. Current searches for physics beyond the Standard Model are sensitive to regions of large invariant mass m γγ > 2m t , so it is essential to quantify the role of the top quark in this region of phase space. This is the primary aim of this section. To that end we have computed the amplitudes for gg → γγ that proceed through a closed loop of heavy quarks, and include details of the calculation in appendix B.
We will perform a detailed phenomenological study of the high invariant mass region in section 4, but to illustrate the importance of the top quark loop we assess its impact on the relevant invariant mass spectrum in figure 3. The results have been obtained for the LHC operating at 13 TeV and under fiducial cuts inspired by the ATLAS collaboration [16] JHEP07(2016)148 that are described in section 4. We show the ratio of the result with n F = 5 light flavors (for the gg initiated pieces only) and the top quark loop included, to the result for n F = 5 light flavors alone. There is a slight decrease in the prediction below the 2m t threshold, due to the effects of a destructive interference, then a steady rise to an asymptotic value. This asymptotic value is of course the result for n F = 6 light quark flavors (without including any modification to the running of α s ) and is simply given by,

Summary
In this section we have presented an overview of the various component pieces of our calculation. For the bulk of this paper we will define our NNLO calculation to only account for five light flavors of quarks. Unless otherwise stated we do not include the NLO corrections to the gg initial state that have been discussed in section 2.2. Instead we refer to these pieces always as σ NNLO + ∆σ N3LO gg,n F . Our default scale choice for the renormalization and factorization scales will be µ = m γγ . We estimate the theoretical uncertainty by varying this central scale by a factor of two in each direction, i.e. µ = 2m γγ and µ = 0.5m γγ . This variation will be indicated by shaded bands in the figures of section 4.

Validation
In this section we compare our results for pp → γγ with those presented in ref. [22]. A summary of cross-sections that have been computed in that work is shown in table 1. To emulate their calculation we impose a series of phase space selection cuts. The cuts on the transverse momenta of the photons depend on their relative size, p hard T > 40 GeV and p soft T > 25 GeV. The photons are also required to be central, |η γ | < 2.5 and in addition we require that the invariant mass of the photon-photon system lies in the interval 20 ≤ m γγ ≤ 250 GeV. Finally at NLO and NNLO we impose the following isolation requirement [63] with n = 1, γ = 0.5 and R = 0.4. We use α = 1/137 and the remaining EW parameters are set to the default values in MCFM. The PDFs are taken from MSTW2008 [64] and are matched to the appropriate order in perturbation theory. The renormalization and factorization scales are mostly set to the invariant mass of the photon pair µ F = µ F = m γγ , although we will also present results for µ F = µ R = m γγ /2 and µ F = µ R = 2m γγ . The results that we obtain from our implementation in MCFM are presented in table 2 and should be compared with the results from ref. [22] that are shown in table 1. Whilst our LO and NLO predictions are in good accord, we find no such agreement for the NNLO cross sections, for any of the choices of scale. The discrepancy is approximately 3pb, or around 8% of the total NNLO prediction. However we do note that the size of the scale variation, i.e. the departures from the central choice, is the same for both calculations. Table 1. Cross sections reported in ref. [22]. Table 2. Cross section results obtained using MCFM. The NLO contribution is always computed using Catani-Seymour dipole subtraction; the NNLO coefficient corresponds to the τ → 0 limit of a calculation using N -jettiness regularization (cf. figure 5). In the NNLO calculation the errors are obtained by adding the fitting and NLO Monte Carlo uncertainties in quadrature.

JHEP07(2016)148
Since we therefore do not agree with the essential results of the existing literature we now describe the further checks that we have performed on our calculation. Several of the ingredients for the below-cut contribution have been reused from previous calculations where good agreement with the literature results was obtained. Specifically, the soft and beam functions have already been used to compute the Drell-Yan and associated Higgs production processes [52,65]. The MCFM predictions for these cross sections agree perfectly with the known results from the literature. The remaining below-cut contribution, the hard function, has been implemented in two independent codes that check both the SCET matching and the proper inclusion of the double-virtual results of ref. [57]. 1 Additionally we have checked that by setting µ 2 = s, and implementing the hard function for a specific scale, we can reproduce the full result by application of the renormalization group equations. This test is extremely non-trivial since the µ 2 dependence occurs both in the finite functions taken from ref. [57] (in their notation, a dependence on S) and also in the matching to the SCET formalism. This check therefore ensures that no mistakes are made in the relative normalization between the two parts of the hard function calculation. For the gg → γγ pieces we have reproduced the results of refs. [20,21], which were implemented previously in MCFM [55]. For the above-cut pieces we have compared our NLO prediction for γγj with the results presented in ref. [66], finding agreement for the isolation procedure used here ("smooth-cone"). We have also checked the analytic calculation of the helicity amplitudes for the real and virtual contributions to γγj production against an in-house implementation of the numerical D-dimensional algorithm [67].
In order to eliminate the N -jettiness slicing procedure as a cause of the difference, we have also implemented Q T -slicing in MCFM. 2 This implementation has been additionally 1 We have adjusted the results of ref. [57] to account for small typos in the manuscript, as detailed in appendix A. 2 The QT -slicing method is based on the same factorization and ingredients that were used in the previous QT -subtraction calculation [23]. Figure 4. The dependence of the NLO cross section on the slicing parameter δ cut . Results are presented using the N -jettiness (circles) (δ cut ≡ τ cut ) and Q T -slicing (triangles) (δ cut ≡ Q cut T ) methods. In both cases the results are normalized to the standard MCFM prediction obtained with Catani-Seymour dipole subtraction, which does not have a slicing parameter dependence. checked, for large values of Q cut T , with a calculation using a completely different setup. The alternate Q T -slicing calculation is implemented using the Sherpa framework [68] and uses the OpenLoops [69] and BlackHat [70][71][72] programs to evaluate the above-cut matrix elements. An obvious cause for concern in either of these slicing-based methods is the dependence on the regulating parameter. When comparing our predictions it is therefore crucial to investigate the dependence of them on this unphysical slicing parameter, either τ cut or Q cut T as appropriate. As a point of reference, we first study the dependence of the total NLO cross section on the slicing parameter in figure 4. To assess the agreement with the known result, we divide the results of these calculations with the one obtained from the existing NLO calculation of MCFM. This implementation of the pp → γγ process [55] uses Catani-Seymour dipoles [73] to regulate the infrared divergences and thus contains no dependence on a slicing parameter. The figure indicates that the slicing results approach the correct cross section, with deviations in the cross section that are O(0.1)% and smaller for τ cut 0.002 GeV or Q cut T 0.04 GeV. This agreement is an additional check of the correctness of the NNLO calculation since the one-loop hard function is also used there.

JHEP07(2016)148
Although the effect of power corrections appears to be milder for Q T -slicing than Njettiness regularization, by around a factor of 20, we note that the computational resources required to perform the calculations at these two points is similar. The resources needed for a computation of a given accuracy is dominated by the calculation of the above-cut contribution, which scales as [42,74], for the N -jettiness calculation. In this equation Q is an appropriate hard scale that is given here by the transverse momentum of the photons. A similar analysis for Q T -slicing yields the result [75], Δσ ���� (δ ��� ) [��] Figure 5. The dependence of the NNLO coefficient ∆σ NNLO on the slicing parameter δ cut . Results are presented using the N -jettiness (δ cut ≡ τ cut ) (circle) and Q T -slicing (triangles) (δ cut ≡ Q cut T ) methods. The dashed lines correspond to the errors associated with the fitting procedure.
Therefore one expects similar computational effort for the two methods when the values of τ cut and Q cut T are related by [74], For Q = 40 GeV one therefore expects the NLO calculation using Q cut T = 0.04 GeV to be as expensive as the one with τ cut = 0.0023 GeV. Figure 5 shows the δ cut dependence for the NNLO coefficient, ∆σ N N LO (cf. eq. (2.1)). It is clear that the dependence is much more pronounced than at NLO. To achieve a 1% accuracy for ∆σ N N LO requires a value of τ cut around 0.002 GeV or Q cut T smaller than about 0.02 GeV. Once again power corrections are less significant for Q T -slicing, but the computing time to achieve equivalent accuracy is comparable in both methods. This is in line with the scaling expected from eq. (3.4). The NNLO results reported in table 2 are obtained from the asymptotic τ → 0 results obtained by a fit to the τ cut dependence that is represented by the solid red line in figure 5. We observe that for values of Q cut T around 1 GeV there is a a local maximum in the NNLO coefficient, which could be mistaken for the onset of asymptotic behavior.
We have communicated our findings with the authors of ref. [22], who have acknowledged a problem with their results presented in ref [22]. The updated version of their code produces results that are consistent with ours, within Monte Carlo uncertainties.

LHC phenomenology
In this section we present results that are relevant for current LHC phenomenology. We first investigate the comparison of our calculation with existing data taken by the CMS experiment with the LHC operating at √ s = 7 TeV. Although such comparisons have already been performed, we believe that this is especially important given the disagreement with the previous NNLO calculation noted in section 3. Additionally, we are able to make the first comparison of the data to a theory prediction that includes both NNLO

JHEP07(2016)148
and ∆σ N3LO gg,n F . We then turn our attention to more recent data taken at √ s = 13 TeV and concentrate on the region of high invariant mass of the diphoton pair, which is relevant for searches for new physics. This region of phase space is particularly interesting given the recent observations of excesses in the data at around 750 GeV [16,17]. For the remainder of this paper we will use the NNLO CT14 PDF set [76] for all predictions (NNLO, NLO, and ∆σ N3LO gg,n F ). The NLO (and ∆σ N3LO gg,n F ) contributions are computed using dipole subtraction and the NNLO coefficients use jettiness regularization with a value of τ cut = 0.002 GeV. From the studies of section 3 we expect this to give us control of the power corrections at the few per-mille level in the total cross-section. We maintain the EW parameters from the previous section, namely α = 1/137.

pp → γγ as a probe of hard QCD
As a benchmark we take the recent study by CMS at 7 TeV [15]. In order to mimic the cuts applied in the experimental analysis we enforce the following phase space selection cuts, Note that the small slice of rapidity that is excluded is due to the design of the CMS detector. In addition we apply isolation cuts to the photon using the smooth cone prescription [63] that does not require an implementation of photon fragmentation. As part of their study CMS compared various smooth cone implementations to that of Diphox, which includes the fragmentation contribution, ultimately employing the following isolation prescription, with = 5 GeV, R 0 = 0.4 and n = 0.05. The rather low value of n results in a fairly weak damping of the collinear singularity present in the calculation as ∆R → 0. Therefore at the cost of deviating from the isolation requirement outlined in ref. [15], we instead use the following definition, with γ = 0.1 and n = 2. We have tuned the values of γ and n such that our NLO smooth cone cross section agrees with the theory prediction obtained at NLO with the CMS isolation experimental requirement and GdRG fragmentation functions [77]. Choosing such a value of n results in a much more efficient Monte Carlo code. A related study of photon plus jets [71] drew similar conclusions. We do not believe that the difference in isolation is a particular cause for concern [58], especially since the cross section has been tuned to a NLO calculation that includes the effects of fragmentation. In principle, the isolation procedure used by CMS in their theory predictions could be chosen in MCFM, but the  Figure 6. The pp → γγ cross section at various orders in perturbation theory, as a function of the LHC operating energy, √ s. Acceptance cuts have been applied, as described in the text. Also shown is the CMS measurement, under the same set of cuts, at 7 TeV [15]. calculation of the corresponding NNLO corrections would require significantly more Monte Carlo statistics to evaluate, with little additional benefit.
We begin by comparing the total cross section as measured by CMS to our prediction using MCFM. The value reported by CMS is, σ CM S = 17.2 ± 0.2 (stat) ± 1.9 (syst) ± 0.4 (lumi) pb , Thus, within the theoretical and experimental uncertainties, the two are in good agreement. Including the NLO corrections to the gg initiated pieces raises the theoretical prediction by around 7%, Since we do not include the full N 3 LO prediction we do not obtain any improvement in the scale variation when including the gg box contributions at NLO. As a brief aside, in figure 6 we show the cross section computed at higher center of mass energies, from the 7 TeV result discussed above to the highest design energy of the LHC, 14 TeV. In the figure we include the cross sections computed at LO, NLO, NNLO and NNLO+gg boxes at NLO. As the order in perturbation theory increases there are sizeable corrections. Going from LO to NLO the cross section increases by around a factor of 4. The corrections going from NLO to NNLO are around 1 Figure 7. The invariant mass of the photon pair m γγ at NLO and NNLO, compared with the CMS data from ref. [15]. The pure NNLO prediction is shown in the left panel, while the result that also includes gg n F contributions that enter at N 3 LO is depicted in the right panel. The lower panels present the ratio of the data and NNLO scale variations to the NNLO theory prediction obtained with the central scale.
gg contributions at NLO increases the cross section by about a further 10%. At the 13 TeV LHC the difference between σ N N LO and σ N N LO + ∆σ N3LO gg,n F is still similarly relevant and it is entirely possible that a measurement will prefer one value over the other. Note that it is not trivially true that σ N N LO +∆σ N3LO gg,n F is a better prediction than σ N N LO since the former is not a complete N 3 LO calculation. The missing pieces are not positive definite, and may reduce the cross section such that σ N 3LO lies completely within the uncertainty bands of the NNLO calculation. It will be interesting to compare the measured cross sections at 13 TeV and 14 TeV to the two predictions to see if indeed σ N N LO + ∆σ N3LO gg,n F does a better job of describing the data than σ N N LO alone.
We now turn our attention to more differential quantities, namely the invariant mass of the photon pair, m γγ (figure 7), the transverse momentum of the γγ system, p γγ T (figure 8), and the azimuthal angle between the two photons, ∆φ γγ (figure 9). We note that, of these predictions, only m γγ is non-trivial at LO since the back-to-back nature of the kinematics at LO means that p γγ T = 0 and φ γγ = π. Such distributions that are trivial at LO are particularly sensitive to higher order corrections. In the bulk of the phase space they first appear at one order higher in α s than the total inclusive cross section. Sadly, most of the distributions made publicly available by the experimental collaborations suffer from this problem. It would be interesting to additionally compare true NNLO observables, such as the transverse momenta and rapidities of the photons, in future analyses at higher energies.
We now examine the predictions for the invariant mass of the photon pair shown in figure 7 in more detail. Note that the transverse momentum cuts on the photons requires  m γγ > 80 GeV at LO, so that the region of this distribution below that value is particularly sensitive to higher order corrections. For all of the figures described here, the plots on the left hand side are obtained using a pure NNLO prediction, while those on the right represent the prediction obtained with the inclusion of the ∆σ N3LO gg,n F contributions. The NNLO prediction does a good job of describing the data obtained by CMS, although the central values are typically a little on the low side compared to data. The situation is improved in the right hand plot, after inclusion of the ∆σ N3LO gg,n F pieces. In particular in the region around 80 m γγ 150 GeV the prediction follows the shape of the data a little more closely.
In figure 8 we turn our attention to the p γγ T spectrum, using the same style as for the m γγ plots. The pure NNLO prediction again describes the data very well, even in the very soft p γγ T < 10 GeV region of phase space. Including the gg pieces at NLO improves the agreement with data in the region p γγ T > 10 GeV. In the soft region of phase space it is difficult to argue that the inclusion of the additional pieces improves the agreement with data. This is understandable since the softest bins are described only after a delicate cancellation between the various real and virtual pieces of the calculation. By only including a subset of the N 3 LO calculation we are unlikely to improve this bin. However in the bulk of the phase space we are typically interested in the types of correction that are sensitive to the staggered phase space cuts. This is exactly the places where we expect the gg → γγg contribution to be important. By including these pieces we therefore do a better job of describing the data. �����/���� Figure 9. As for figure 7, but for the azimuthal angle between the two photons, ∆φ γγ .
The situation with the ∆φ γγ distribution is similar. The NLO prediction for this observable does a very bad job of describing the CMS data. However by including the NNLO corrections we get much closer to the data, whilst still observing deviations from the experimental data of order 20%. Thus, this observable clearly requires at least a full N 3 LO prediction to match the experimental data. However, our partial prediction does not do much better. Again we are exposed to the LO phase space sensitivity in the bins around π where it is entirely possible that reasonably large corrections from the three-loop triple virtual and real-double virtual may drive the theoretical prediction down towards the data.

Studies of γγ at high invariant masses
One of the most interesting phenomenological aspects of the diphoton production channel during Run II at the LHC is its ability to search for new resonances that may manifest themselves in the m γγ spectrum. In particular a recent observation of an excess around 750 GeV in the ATLAS experiment [16], with a smaller excess in the same region reported by CMS [17], has caused considerable excitement in the theoretical community. In these analyses the Standard Model background is accounted for by using a data-driven approach that fits a smooth polynomial function to the data across the entire m γγ spectrum. A resonance might then be observed as a local excess in this spectrum, deviating from the fitted form. Although well-motivated, one might be concerned that the spectrum may not be correctly modeled at high energies, where there is little data, and that small fluctuations could unduly influence the form of the fit and result in misinterpretation of the data. Such worries could be lessened by using a first-principles theoretical prediction for the spectrum and it is this issue that we aim to address in this section.
As a concrete example, we will produce NNLO predictions for the invariant mass spectrum at high energies using cuts that are inspired by the recent ATLAS analysis [16].

JHEP07(2016)148
�/����(�� � ) Figure 10. The ratio of various different theoretical predictions to the NNLO n F = 5 differential cross section. The different predictions correspond to: the inclusion of the top quark gg → γγ box diagrams (green), the ∆σ N3LO gg,n F correction (red) and the ∆σ N3LO gg,n F and the top boxes with the ∆σ N3LO gg,n F correction re-scaled by the ratio K(m t ) described in the text (blue).
Specifically, these are: |η γ | < 2.37, excluding the region, 1.37 < |η γ | < 1.52. (4.6) We will only be interested in the region m γγ > 150 GeV, so these represent hard cuts on the photon momenta. The small region of rapidity that is removed corresponds to the transition from barrel to end-cap calorimeters. We maintain the same isolation requirements as the previous section, which again differs slightly from the treatment in the ATLAS paper. Our first concern is to address the impact of the gg pieces at NLO, represented by the contribution ∆σ N3LO gg,n F defined previously, and the contribution of the top quark loop. We summarize our results in figure 10, in which we present several different theoretical predictions, each normalized to the the default NNLO prediction with 5 light flavors. The first alternative is one in which the NNLO prediction is augmented by the inclusion of the top loops, i.e. the gg contribution corresponds to σ gg (m t + 5l f ) in the notation of section 2.3. In the second prediction we use the result for five light flavors but add the NLO corrections to the gg channel, i.e. the term ∆σ N3LO gg,n F . For the final alternative we include the top quark loop contribution and attempt to account for the NLO corrections to all gg loops by rescaling the ∆σ N3LO gg,n F result by a factor K(m t ) that is given by, This collection of predictions covers a range of theoretical options that may extend the NNLO predictions described in the previous sections. The top loops, illustrated by the

JHEP07(2016)148
green curve in the figure typically represent around a 1% effect across the invariant mass range of interest. For m γγ < 2m t there is a destructive interference, which reduces the cross section, whilst at higher energies there is a small enhancement. Therefore, although the top loops are an important contribution in terms of the n F box loops (as shown in section 2), they are not particularly important in the total rate. At this order the gg pieces reside in the Born phase space, which is particularly impacted by the staggered cuts at high m γγ .
As we found in the previous section the effects of the NLO corrections to the gg pieces are larger, however their effects are much more pronounced at lower invariant masses. By the time invariant masses of order 500 GeV are probed, the corrections are 2% or smaller. The attempt to model the combined effect of corrections to both the light-quark and top quark loops shows, as expected, the largest deviations from the NNLO(5 f ) prediction. However the deviations are still of order 3% or smaller in the high invariant mass region. Therefore, although the corrections to the gg loops and the effect of the finite top quark mass can have about a 6% effect at invariant masses around 200 GeV, the effect at higher masses is somewhat smaller. Since we aim to compare the ATLAS data, which is not corrected for fakes or identification efficiencies, to our parton-level prediction we are not concerned about effects at this level. As a result we will simply use the most consistent prediction, 3 corresponding to NNLO(5 f ), for comparison with the fitting function used by ATLAS.
We compare our NNLO prediction to the ATLAS data in figure 11. We note that to properly compare our prediction to the data requires knowledge of both the fake rate and the photon efficiencies and acceptance corrections of the ATLAS detector. To try to minimize the impact of such corrections we simply compare the shape of the ATLAS data to the shape of our NNLO prediction, i.e. we normalize our prediction to 1/σ N N LO and the ATLAS data to 1/N events . From this comparison we can draw several conclusions. First, we note that our prediction is in excellent agreement with the overall shape of the data, indicating that the theoretical prediction for the shape of the m γγ distribution could easily be used in place of the somewhat arbitrary fitting functions currently employed. Second, the excellent agreement in shape suggests either a low number of fakes, or that the fake events are distributed with a similar shape to the Standard Model prediction for the γγ spectrum. Of course a combination of these two explanations is also possible. Finally, and most excitingly, a comparison to the fitting function presented in ref. [16] illustrates that there is no significant hardening from the prediction of the SM compared to the form of the fitting function used in the ATLAS experiment. This can clearly be seen upon comparison with figure 1 in ref. [16]. For instance, both the ATLAS fit and our NNLO prediction pass directly through the data in the 1090 GeV bin, and just under the central value in the 690 GeV bin. Therefore we can conclude that the interpretation of an excess of events around 750 GeV appears to be supported by a first-principle calculation within the SM. It is not diluted by a hardening of the SM spectrum relative to the fitting function used in the analysis. If the excess is confirmed, NNLO predictions for the shape of the irreducible background will be able to significantly enhance analyses designed to discriminate between different model hypotheses, by providing predictions for the properties of background events that cannot be captured by a simple spectrum fit.

Conclusions
The process pp → γγ is a flagship process for Run II phenomenology. Besides its intrinsic interest as a tool to understand the perturbative nature of QCD at high energies, it represents an important background in studies of the Higgs boson that are a cornerstone of the Run II physics program. In addition it is a clean and well-measured final state that can be used in the search for new heavy resonances. Often these analyses require staggered photon transverse momentum cuts that induce large corrections at higher orders in perturbation theory. Essentially the NLO prediction behaves like a LO prediction since the staggered cuts are first accessible at this order. This therefore necessitates the inclusion of NNLO corrections to capture the corrections to the rate in this larger phase space and hence adequately describe data.
In this paper we have presented a NNLO calculation of the process pp → γγ and studied the phenomenology of this process at the LHC. We have used the recently-developed JHEP07(2016)148 N -jettiness slicing procedure to manage the infrared singularities present in the NNLO calculation and have implemented the calculation into the Monte Carlo MCFM. The calculation will be made available in a forthcoming release of the code. Given the signficant effect of the NNLO corrections to this process, our slicing procedure is subject to large power corrections and care must be taken to ensure that a small enough value of the slicing parameter is employed. We have compared our results to an existing calculation of the same process and found that we could not reproduce the results present in the literature, despite extensive testing and investigation. However we have communicated with the authors of ref. [22] and believe that, after correction of a bug in their numerical code, their results will be consistent with ours.
We have used our calculation to compare to data obtained at 7 TeV by CMS and to 13 TeV data collected by the ATLAS experiment. The latter is particularly exciting given the excess reported in the data at around m γγ ∼ 750 GeV. We found that the shape of our NNLO prediction does a good job of describing the experimental data, and simultaneously has a good agreement with the fitting function used by the ATLAS collaboration. We therefore do not expect further data at high energies to dramatically alter the form of the fit used by the collaboration. Furthermore, we do not believe that the excess is due to the use of a fitting function that underestimates the prediction of the SM at high invariant masses.
A Ingredients for qq → γγ at NNLO A.1 Below τ cut : hard function The virtual matrix elements needed to compute qq → γγ to O(α 2 s ) accuracy can be found in ref. [57]. In order to be utilized in N -jettiness slicing these results must be translated into the form of a SCET hard function. This can be achieved using the procedure outlined, for instance, in refs. [78,79]. We begin by defining the UV-renormalized matrix element as follows, where α s is the renormalized strong coupling, and α is the (bare) electromagnetic coupling. The matrix elements are defined in terms of Mandelstam invariants s = (p 1 + p 2 ) 2 , t = (p 1 + p 3 ) 2 and u = (p 1 + p 4 ) 2 , with s + t + u = 0. For the process under investigation p(−p 1 ) + p(−p 2 ) → γ(p 3 ) + γ(p 4 ) we have s > 0 while t, u < 0. Following the notation of JHEP07(2016)148 ref. [57] we define the matrix element squared as follows, Expanding to O(α 2 s ) we define, In terms of the matrix elements defined above we have The aim of this section is to re-write the above expressions in the SCET renormalized form, which is obtained via the following re-definitions [78,79]  I (1) ( ) and I (2) ( ) are obtained via Catani's IR-subtraction formula [80]. For the qqγγ process under investigation here I (1) ( ) and I (2) ( ) are defined as follows In the above equation the H 2 R.S. ( ) is a scheme dependent function, containing 1/ poles that for this process is defined as [80]  T f n f , (A.14) and Z is defined, for our process and order in perturbation theory as [78,79]  For brevity we present the results obtained at µ 2 = s; the full scale dependence may be obtained by inspection of the distributed MCFM routines, or analytically by appropriate usage of the renormalization group equations. The hard function for the NLO process is given bỹ where we have introduced the following notation [57] X = log − t s , Y = log − u s , (A.24)

A.2 Above τ cut
For τ > τ cut the calculation corresponds to an NLO calculation of the γγj process. An implementation of this process and the γγγ process in the MCFM framework was presented in ref. [58]. We use the results of this calculation, which corresponds to an analytic calculation using helicity amplitudes and D-dimensional unitarity methods to obtain our above-τ cut pieces. We refer the interested reader to ref. [58] for more details.
B gg → γγ: m t loops In this section we present the calculation of the gg → γγ process that proceeds through a top-quark loop. We use the spinor helicity formalism to define our amplitudes, and refer readers unfamiliar with the notation and conventions to one of the many comprehensive reviews of the topic (for instance ref. [82]). Throughout our calculation of these amplitudes we made frequent use of the Mathematica package S@M [83].
We define the partial amplitude for this process as follows,

JHEP07(2016)148
In the notation of the QCDLoop library [84], which we use to evaluate the integrals, I 4 (s, t, m 2