Precise predictions for photon pair production matched to parton showers in GENEVA

We present a new calculation for the production of isolated photon pairs at the LHC with NNLL$'_{\mathcal{T}_0}$+NNLO accuracy. This is the first implementation within the GENEVA Monte Carlo framework of a process with a nontrivial Born-level definition which suffers from QED singularities. Throughout the computation we use a smooth-cone isolation algorithm to remove such divergences. The higher-order resummation of the 0-jettiness resolution variable $\mathcal{T}_0$ is based on a factorisation formula derived within Soft-Collinear Effective Theory which predicts all of the singular, virtual and real NNLO corrections. Starting from this precise parton-level prediction and by employing the GENEVA method, we provide fully showered and hadronised events using PYTHIA8, while retaining the NNLO QCD accuracy for observables which are inclusive over the additional radiation. We compare our final predictions to LHC data at 7 TeV and find good agreement.


Introduction
The production of two isolated photons is one of the most interesting processes to study at the Large Hadron Collider (LHC), both to further test the Standard Model (SM) and to search for new exotic signatures, e.g. by looking for heavy resonances in the diphoton invariant mass spectrum [1,2]. Precise experimental measurements are available due to the clean final state and the relatively high rate of production for this process. Efforts in the study of the diphoton final state were especially boosted by the discovery of a Higgs boson decaying into two photons at the LHC [3,4]. For all these reasons, recent experimental analyses of diphoton production have been carried out by the ATLAS experiment both at 7 TeV [5,6], 8 TeV [7] and 13 TeV [8], and by CMS at 7 TeV [9,10]. Given the relevance of this process for the LHC experimental analyses, precise calculations including higher-order QCD corrections are required from the theoretical side.
In this paper we study the production of "prompt" photon pairs which are directly produced in the hard scattering interaction. A second production mechanism involves the radiation of photons in a jet during the fragmentation process. This second mechanism can be suppressed by isolating the photons from the final-state hadrons with either the fixedcone or the smooth-cone (a.k.a. Frixione) isolation procedure [11]. 1 The fixed-cone isolation is widely used in the experimental analyses due to the simplicity of its implementation. It cannot, however, completely eliminate the fragmentation contribution without spoiling infrared (IR) safety. The smooth-cone isolation can achieve this in an IR-safe manner, thus simplifying the computation of the radiative corrections.
Next-to-leading-order (NLO) QCD corrections for the process pp/pp → γγ + X were first computed and implemented in the fully differential Monte Carlo program DIPHOX [13] which included the calculation of both the "direct" and the fragmentation contributions to the cross section. The leading-order (LO) calculation [14] of the gluon channel gg → γγ contribution, which is formally a next-to-next-to-leading-order (NNLO) effect, was also implemented in the DIPHOX code. The inclusion of this channel has a relatively large impact on the cross section, due to the enhanced gluon parton distributions at the LHC. The NLO corrections to the gluon channel (which amount to N 3 LO effects relative to the LO process qq → γγ) were implemented in the parton-level Monte Carlo programs 2gammaMC [15] and MCFM [16]. Recently, mass effects due to the top-quark loops were studied in the gluon fusion channel at NLO accuracy in Refs. [17,18].
In the context of smooth-cone isolation, where fragmentation functions are not required, the first NNLO calculation for diphoton production was carried out at the fully differential level using the q T subtraction method of Ref. [19]. It has been implemented in the numerical codes 2γNNLO [20] and Matrix [21]. An independent calculation at NNLO accuracy for diphoton production was also performed in Ref. [22] and implemented in the MCFM program [23]. This calculation uses the N -jettiness subtraction method [24,25] derived within Soft-Collinear Effective Theory (SCET). Phenomenological studies of photon isolation at NNLO accuracy were carried out in Ref. [26] by using the aforementioned programs 2γNNLO and Matrix. More recently, an independent NNLO analysis using the antenna subtraction method [27][28][29] was presented in Ref. [30], highlighting the importance of the choice of the isolation criterion and of the scales in estimating the theoretical uncertainties. Electroweak (EW) corrections for the diphoton production process at the LHC were computed in Refs. [31,32]. Their effect for inclusive observables is at the level of a few percent, and increases as one imposes higher transverse-momentum cuts on the photons.
Resummed calculations in the region of small transverse momentum of the photon pair have been carried out at next-to-next-to-leading logarithmic (NNLL) accuracy and are available in the programs RESBOS [33][34][35], 2γRes [36] and reSolve [37]. Very recently, the resummation was pushed to N 3 LL accuracy in CuTe-MCFM [38], extending the matching to NNLO accuracy.
Shower Monte Carlo (SMC) simulations for this process at NLO accuracy matched to the Sherpa [39] parton shower program were presented in Ref. [40]. A similar accuracy obtained via matching to the Herwig [41] parton shower using the Powheg approach [42,43] was presented in Ref. [44].
In this paper we implement a new fully differential NNLO event generator for the diphoton production process pp → γγ + X. We improve on the existing fixed-order (FO) calculations by including the NNLL T 0 resummation of the 0-jettiness resolution variable T 0 (beam thrust) within the Geneva Monte Carlo framework [45][46][47][48]. The resummation of the 0-jettiness variable is carried out by using a factorisation formula [49,50] valid at small T 0 derived in SCET. This is the first time that such a resummation has been presented in the literature for this process. We employ the smooth-cone isolation algorithm to define IR-finite events, ensuring that the application of this cut is compatible with the T 0 resummation. In addition, we provide fully showered and hadronised events by matching partonic events to the Pythia8 [51] parton shower while retaining the FO accuracy for inclusive observables. Diphoton production extends the list of NNLO+PS processes, such as Drell-Yan [47] and associated Higgs production [52] and decay [53], which were previously implemented in the Geneva Monte Carlo program. Alternative approaches to reach NNLO+PS accuracy are actively being developed [54][55][56][57][58][59].
This paper is organised as follows. We discuss the process definition and the photon isolation criteria in sec. 2. In sec. 3 we recap the theoretical background for the 0-jettiness resummation in SCET, while in sec. 4 we provide the details of the implementation in the Geneva framework. We present our results in sec. 5, including the comparison to 7 TeV LHC data. Finally we give our conclusions and outlook in sec. 6.

Process definition and photon isolation
In this paper we consider the process where the photon pair is produced through the hard scattering interaction (i.e. direct photon production). In order to avoid QED singularities when a photon is collinear to an initial-state parton, one needs to impose kinematic cuts on the transverse momenta of the photons. We prefer to use asymmetric cuts on the harder (γ h ) and the softer (γ s ) photon, i.e. p γ h T ≥ p γ h T,cut and p γs T ≥ p γs T,cut (where p γ h T ≥ p γs T ) which avoids issues in the fixed-order predictions with symmetric cuts [60][61][62]. These cuts are also employed in experimental analyses to exclude events in which the photons are emitted in the uninstrumented region parallel to the beam line.
As mentioned in the introduction, final-state photons can also be produced via the fragmentation of a quark or a gluon into a photon. Usually, photons produced in a fragmentation process are distinguishable from the direct photons since they lie inside hadronic jets. Direct photons can therefore be separated from the rest of the hadrons in the event through an isolation procedure. A possible isolation mechanism is to construct a cone with fixed radius R iso around the direction of the photon candidate. One then restricts the amount of transverse hadronic (partonic) energy E had T (R iso ) inside the cone. In particular, a photon is considered to be isolated if E had T (R iso ) is smaller than a certain value E max T , usually parameterised as a (linear) function of the transverse energy E γ T of the photon and a fixed numerical value E thres T [63]: Due to its simplicity, the fixed-cone isolation procedure is currently used in all experimental measurements of processes involving photons. This method has the theoretical drawback of being sensitive to the fragmentation contributions since configurations with a photon collinear to a final-state parton are still allowed. Indeed, to completely remove the collinear QED divergences, one should require that absolutely no hadronic energy is allowed within the cone. Unfortunately, this condition is not IR safe since it forbids emissions of soft partons inside the cone, therefore spoiling the cancellation of QCD divergences. The smooth-cone isolation procedure instead overcomes this problem by considering a continuous series of smaller sub-cones with radius r ≤ R iso , in addition to the outer cone of radius R iso . The isolation condition then requires that where the isolation function χ(r; R iso ) must be a smooth and monotonic function which vanishes when r → 0. This requirement implies that the hadronic activity is reduced in a smooth way when approaching the photon direction, and the parton radiation collinear to the photon is completely absent at r = 0. Hence, the fragmentation component is eliminated while soft radiation is still permitted in a finite region away from the collinear limit, making the cross section IR safe. The standard choice for the χ-function is where the exponent parameter n is usually set to n = 1. Other isolation functions satisfying these criteria are possible and have been employed in the literature [26]. Though the use of a smooth-cone isolation procedure has the positive effect of simplifying theoretical calculations, its implementation in the experimental analyses is complicated by the granular nature of the detector. The direct comparison of the theoretical predictions with data is therefore challenging. One possible solution is to adjust the free parameters of the smooth-cone isolation algorithm to reproduce the effects of the fixed-cone procedure, making a comparison at least feasible. A potentially better approach, which has been recently investigated in Refs. [12,64], is the introduction of a hybrid-cone isolation procedure. In this case the theoretical calculation is initially carried out using the smooth-cone isolation with a very small radius parameter R iso , such that only a tiny slice of phase space around the photon direction is removed. In a second step, the fixed-cone isolation procedure with a larger radius R R iso is applied to the events which passed the smooth-cone criterion. In other words, one initially applies very loose smooth-cone isolation cuts which are then tightened by the fixed-cone procedure. This makes the resulting events directly suitable for experimental analyses. 2 In this paper we investigate both the smooth-cone and the hybrid-cone isolation procedures. We use the former to compare to results obtained with the Matrix code [21] in sec. 4.5, while the latter is used for a direct comparison to LHC data in sec. 5.

Resummation in Soft-Collinear Effective Theory
In this work, we use the N -jettiness [65] resolution variable to discriminate between resolved emissions with different jet multiplicities. Given an M -particle phase space point Φ M with M ≥ N , the N -jettiness is defined as where the index k runs over all strongly-interacting final-state particles and whereq i = n i = (1, n i ) are light-like reference vectors parallel to the beam and jet directions, defined in the rest frame of the γγ colour-singlet state. The limit T N → 0 describes an event with N pencil-like hard jets, where the unresolved emissions can either be soft or collinear to the final-state jets or to the beams. In the case of colour-singlet production (e.g. Drell-Yan, V H, γγ, . . . ) the relevant resolution variable is the 0-jettiness (beam thrust). Starting from the general definition in eq. (3.1), the expression for the 0-jettiness can be considerably simplified to where | p kT | and η k are the transverse momentum and the rapidity of the emission p k and Y is the rapidity of the colour-singlet state. The introduction of a jet resolution variable sets a new dynamical energy scale in the problem, which can in principle differ from the other physical scales of the process. As such, large logarithms of the ratios of these scales may appear in the cross section which spoil the convergence of the perturbative expansion and must therefore be resummed. For 0-jettiness this happens in the small T 0 /Q region, with Q being a typical hard scale of the process, e.g. the invariant mass of the colour-singlet final state. In this region, the cross section differential in the Born phase space Φ 0 and T 0 obeys a factorisation formula which was originally derived in Refs. [49,50] for the Drell-Yan process, but which can be extended for diphoton production as The sum in the equation above runs over all possible qq pairs ij = {uū,ūu, dd,dd, . . .}. For this process we need to impose process-defining phase space restrictions Θ PS (Φ 0 ) in order to have a finite cross section. In particular, we apply p T cuts on each of the photons as well as isolation cuts to eliminate final-state collinear QED singularities. 3 Other cuts, such as on the photons' rapidities or on the invariant mass of the pair, can be imposed at the analysis level but are not needed to define IR-finite cross sections. The factorisation formula depends on the hard H γγ ij , soft S and beam B i,j functions which describe the square of the hard interaction Wilson coefficients, the soft real emissions between external partons and the hard emissions collinear to the beams respectively.
The hard functions H γγ ij (Q 2 , t, µ) are process-dependent objects and encode information about the Born and virtual squared matrix elements. In order to achieve NNLL accuracy, they are needed up to two-loop order. They are regular functions of the Mandelstam invariants Q 2 = s and t and can be extracted from the two-loop squared amplitude expressions [66], after subtracting the IR poles (as explained in detail in appendix A).
The B i (r, x, µ) are the inclusive (anti-)quark beam functions [49]. They depend on the virtualities r a,b of the initial-state partons i and j annihilated in the hard interaction and on the momentum fractions x a,b . These can be written in terms of the diphoton rapidity Y γγ and invariant mass M γγ as where E cm is the hadronic centre-of-mass energy. The beam functions are calculated via an operator product expansion and are defined as The perturbatively computable parts of the above equation are the matching coefficients I ik (t a , z a , µ) which describe the collinear virtual and real initial-state radiation (ISR) emissions. The function f k (ξ a , µ) represents the usual parton distribution function (PDF) for parton k with momentum fraction ξ a . The matching coefficients I ik (t a , z a , µ) were computed to NNLO accuracy in Ref. [67]. S(k, µ) is the quark hemisphere soft function for beam thrust and has been computed to the required NNLO accuracy including the scale independent terms in Refs. [68,69]. The hard, beam and soft functions which appear in eq. (3.3) are single-scale objects. This means that when they are evaluated at their own characteristic scales no large logarithmic corrections are present in their fixed-order perturbative expansions. However, eq. (3.3) must be evaluated at a single common scale µ. In order to do this, the separate components must be evolved via renormalisation group (RG) equations from their own characteristic scales to the common scale µ, which results in the large logarithms being resummed. This proceeds via convolutions of the single scale factors with the evolution functions U i (µ i , µ). The resummed formula for the T 0 spectrum is then given by where the convolution between the different functions is written in a schematic form. The scale setting procedure will be explained in the next section, where we will introduce the profile functions which are employed to switch off resummation outside its kinematical range of validity. In order to reach NNLL accuracy, we need to know the boundary conditions of the evolution, namely the hard, beam and soft functions up to NNLO accuracy, and the cusp (non-cusp) anomalous dimensions up to three-(two-)loop order. The expressions for the anomalous dimensions up to the required order can be found in Refs. [25,[70][71][72][73].

Implementation within the GENEVA framework
In this section we review the Geneva framework and present the implementation of the diphoton production process, highlighting the main differences compared to the processes that have previously been implemented. More details on the general features of the Geneva method are given in Refs. [46,47,52]. For all numerical results in this paper we evaluate all the matrix elements up to the one-loop level with the amplitudes provided by the OpenLoops 2 package [74]. We include the top-quark mass dependence in the NLO 1 calculation and set m t = 173.2 GeV. In the present implementation, however, we do not include the diagrams with a closed top-quark loop in the hard function. These effects begin to contribute at O(α 2 s ) and are currently unknown in the qq channel. 4 If not stated otherwise, in this section we use the following settings: we use the MMHT2014nnlo68cl PDF set [76] from LHAPDF [77] and the α(0) input scheme to fix EW couplings and set α −1 (0) = 137. We apply the kinematic cuts  We set µ FO = M γγ or µ FO = M T γγ ≡ M 2 γγ + p 2 T,γγ and indicate the respective choice with the results. We consider only the qq channel contribution up to and including sec. 4.7 and take into account also the loop-induced gg channel contribution thereafter.

The GENEVA method
The final aim of a fully exclusive event generator is to produce physical events, where all of the IR divergences are cancelled on an event-by-event basis. Following the Geneva method, these events are assigned a Monte Carlo (MC) cross section dσ MC N according to the value of the N -jettiness resolution parameter T N . This encodes the probability of producing a phase space point Φ N with N partonic jets: all contributions from unresolved emissions below a certain resolution cutoff T N < T cut N are included in dσ MC N . For the process at hand at NNLO, the exclusive cross sections for events with 0 and 1 jet, and the 2-jet inclusive cross section are defined by the cutoffs on the T 0 and T 1 resolution variables as Since one integrates over the unresolved emissions below the cutoffs, the definitions of the partonic jets used above depend on the phase space maps Φ N (Φ M ) (with N ≤ M ) which project the unresolved M -body phase space points onto Φ N . Using eq. (4.3) the cross section for a generic observable X is written as where M X (Φ N ) is the measurement function that computes the observable X for the N -parton final-state point Φ N . The cross section defined above is not equivalent to that obtained from a fixed-order calculation. Indeed, for any unresolved emission, the observable is computed on the projected point Φ N (Φ M ) rather than the exact Φ M point. Since the resulting difference vanishes in the limit T cut N → 0, it is advisable to choose this cutoff to be as small as possible. This, however, introduces large logarithms of T N and T cut N , which need to be resummed in order to obtain physically meaningful results.
We perform the resummation of T 0 at NNLL accuracy, matching it to a NNLO 0 5 calculation. This has the positive feature of correctly describing the spectrum both in the small-T 0 region, where the resummation dominates, and also in the large-T 0 region, where the correct dependence is given by the fixed-order result. Depending on the final-state jet multiplicity, the Geneva method requires one to evaluate the resummed and resummed-expanded terms in the cross sections on projected phase space points of lower multiplicity. Even these projected configurations are required to satisfy the cuts which avoid the QED singularities. We use the symbol Θ proj ( Φ N ) (and Θ proj ( Φ N ) for its complement) to indicate this set of phase space restrictions acting on the higher dimensional Φ N +1 phase space due to the cuts on the projected configuration Φ N . In practice, this means that when a term in the cross section, evaluated at a Φ N +1 phase space point, is multiplied by Θ proj ( Φ N ), the Φ N +1 phase space point is projected onto a Φ N point and the cuts are applied on this lower dimensional space. If the projected configuration does not pass the cuts, the initial Φ N +1 configuration is excluded. Notice that the separation realised by the introduction of Θ proj and Θ proj is not usually required in a fixed-order calculation. Nonetheless, we choose to perform it even when evaluating the FO contribution and split this into two parts. Of these, the first, singular part enters in the 0-jet cumulant, eq. (4.10), while the second, nonsingular part enters in the 1-jet cross section below the T cut 0 , eq. (4.14). Since the resummation for 0-jettiness is carried out at NNLL accuracy in Geneva, meaning that it contains all of the singular corrections in T 0 up to O(α 2 s ), we can write the 0-and 1-jet cross sections as where dσ NNLL /dΦ 0 dT 0 is the resummed T 0 spectrum and dσ NNLL /dΦ 0 (T cut 0 ) is its integral. In the above equation we introduced a splitting probability function P(Φ 1 ) which satisfies the normalisation condition to make the T 0 spectrum fully differential in Φ 1 . The nonsingular contributions are given by The terms in squared brackets are the expanded expressions to O(α 2 s ) of the resummed cumulant and spectrum. The NLO 1 term refers to the NLO corrections to the diphoton plus jet production process. The projected Φ 1 → Φ 0 point in the above equation is obtained through a FKS projection. After explicitly writing the FO contributions to the cross sections we obtain where B 1 and B 2 are the 1-parton and 2-parton tree-level contributions respectively, V 0 and V 1 correspond instead to the 0-parton and 1-parton one-loop contributions while W 0 is the two-loop contribution. Notice that the matrix elements appearing in the formula above are implicitly assumed to be UV-renormalised and IR-subtracted such that they are finite. We refrain from explicitly writing the subtraction terms in order to make the notation less cumbersome. We also introduced the notation 12) to indicate that the integration over a region of the M -body phase space is done keeping the N -body phase space and the value of some specific observable O fixed, with N ≤ M . The Θ O (Φ N ) term in the previous equation limits the integration to the phase space points included in the singular contribution for the given observable O. Since the resummed and resummed-expanded contributions are differential in T 0 , the phase space integration of the 2-parton contribution in eq. (4.11) should be parameterised in such a way that it is also differential in T 0 . Indeed the projection In this way all of the terms in the inclusive 1-jet cross section (eq. (4.11)) can be evaluated at the same value of T 0 , and the point-wise cancellation of the singular T 0 contributions is achieved. The symbol Θ T (Φ 2 ) defines the region of Φ 2 which can be projected onto the Φ 1 phase space through the map Φ T 1 (Φ 2 ). The non-projectable region of Φ 2 , which corresponds to nonsingular events, will be included in the 2-parton event cross section below. To this end we introduce the following notation. We identify the set of cuts due to the particular map used by Θ map (while its complement is defined as Θ map ). Nonsingular contributions coming from non-projectable regions (either because they fail the isolation cuts or because they result in an invalid flavour projection) are also present in the case of a Φ 1 → Φ 0 projection. We explicitly include them in our formula for the cross section below the T cut 0 as dσ mc where we indicate the failure to obtain a valid flavour projection with the symbol Θ FKS map and the failure to satisfy the isolation cuts with Θ proj ( Φ 0 ).
Starting at O(α 2 s ) a new channel opens up and contributes to the diphoton production process at the LHC. This is the loop-induced gg channel "box" contribution, which is finite and included at the order we are working, keeping the full dependence on the top-quark mass. In the present work we do not implement any T 0 resummation for this channel and only add its effects at fixed order in the δ(T 0 ) term. We leave the resummation of this channel, which formally also starts contributing at NNLL accuracy, to future work.

Profile functions
In this section we describe the scale choices for the hard, beam and soft functions in eq. (3.7). The cross section for the diphoton production process is affected by large logarithms of T 0 /Q when T 0 Q and Q = M γγ is the hard scale of the process. However, for larger values of T 0 , the logarithmic terms become numerically small, and the cross section is dominated by the FO result in that region of phase space. Therefore, we need to switch off the resummation before reaching this region. To this end, it is helpful to investigate the numerical relevance of the different terms which contribute to the cross section as a function of T 0 . In Fig. 1 we perform this comparison at different fixed orders. 6 In both panels we show the FO calculation, the leading singular terms from the expansion of the resummed formula to the relevant fixed order and the difference between the two. The latter amounts to a nonsingular contribution to the cross section in the limit T 0 → 0. As expected from the leading-power factorisation formula, each resummed-expanded contribution perfectly reproduces the singular behaviour of the corresponding FO term and acts as a subtraction to remove the IR divergences. Indeed, the nonsingular contributions vanish on a logarithmic scale as T 0 → 0. One may also notice that the nonsingular terms become of similar size to the singular contributions at small values of T 0 , in the range of a few GeV. This behaviour is something of a novelty compared to the previous Drell-Yan and V H production calculations and appears to be peculiar to the diphoton production process. Although the nonsingular terms are power suppressed for small T 0 values, the photon isolation procedure might become the primary source of power corrections, enhancing their contribution and making them numerically large [38,[78][79][80].
Resummation is carried out via RGE evolution and can be switched off by setting all of the scales equal to a common nonsingular scale µ NS = µ S = µ B = µ H . We do not evolve the hard function, which is always evaluated at µ H = µ NS [81,82]. Instead we evolve the soft and the beam functions from their characteristic scales up to the hard scale. This is achieved by employing profile scales µ B (T 0 ) and µ S (T 0 ) which ensure a smooth transition between the resummation and the FO regimes. Explicitly, where the common profile function f run (x) is given by [83] f run (x) = This functional form ensures the canonical scaling behaviour as in eq. (3.6) for values below x 1 and turns off resummation above x 3 . After considering that the invariant mass distribution peaks in the range 50-80 GeV (depending on the specific cuts that are applied) and that the nonsingular corrections in Fig. 1 become of the same size of the singular at T 0 ∼ 1 − 3 GeV, we choose the following parameters for the profile functions: In the resummation region the nonsingular scale µ NS must be of the same order as the hard scale of the process M γγ , while in the FO region it can be chosen to be any fixed or dynamical scale µ FO . One can, for example, set it either to M γγ or to the transverse mass of the photon pair M T γγ . We estimate the theoretical uncertainties for the FO predictions by varying the central choice for µ NS up and down by a factor of two and take for each observable the maximal absolute deviation from the central result as the FO uncertainty. For the resummation uncertainties, we vary the central choices for the profile scales µ B and µ S independently while keeping µ H = µ NS fixed. This gives us four independent variations. In addition, we consider two more profile functions where we shift all the x i transition points together by ±0.05 while keeping all of the scales fixed at their central values. Hence we obtain in total six profile variations. We consider the maximal absolute deviation in the results with respect to the central prediction as the resummation uncertainty. The total perturbative uncertainty is then calculated by adding the FO and the resummation uncertainties in quadrature.
As explained in detail in Refs. [47,52], the T 0 integration of the resummation formula (eq. (3.7)) and the procedure of choosing the scales are operations which do not commute with each other. The expression for the cumulant is not, therefore, exactly the same as the integral of the T 0 spectrum, since the profile scales have a functional dependence on T 0 . To obtain an expression for the resummed cumulant instead, one must first integrate the expression in eq. (3.7) for the resummed T 0 distribution and then choose the scales using the same profile scales but with the T 0 replaced by T cut 0 . For example the canonical scales assume the values The difference between the two results for inclusive quantities is formally beyond NNLL accuracy [84]. However, the numerical effect can be large and, in order to preserve the value of the NNLO 0 cross section, we follow the prescription for fully inclusive quantities discussed in Refs. [47,52] which amounts to adding the contribution where κ(T 0 ) and µ h (T 0 ) are smooth functions. The µ h (T 0 ) are new profile scales which are chosen to turn off resummation earlier than the normal profile scales in the resummed calculation, in order to maintain the accurate description of the tail of the T 0 spectrum. Indeed, in the FO region we have that µ h (T 0 ) = Q and the difference in the brackets of eq.
The value of κ(T 0 → 0) can be tuned such that after integration of the T 0 spectrum together with eq. (4.19), the correct inclusive cross section is obtained. A more detailed explanation on the precise structure of this additional term can be found in Refs. [47,52].

Comparison with standard resummation and matching
The implementation of resummation in the Geneva framework takes a very different perspective compared to the usual resummation approach. The latter is normally carried out with the purpose of improving the description of a single particular observable in the limit where it approaches very small or very large values compared to other scales present in the process. In practice, this involves a scan over the T 0 spectrum by directly evaluating the resummation formula (eq. (3.7)) on a Φ 0 phase space point. Hence the information about the physical event which generates a particular value of the resummed variable is not retained in the calculation. In an event generator such as Geneva, however, one starts by generating e.g. Φ 1 events and calculating the value of the resummed variable, in this case T 0 , resulting from that particular event configuration. At this point one has to perform a projection to the lower-multiplicity phase space to evaluate the resummation formula. This second method is more flexible because it allows one to access the event information in a fully differential way and to associate to each event a resummed weight, also allowing the matching to a FO calculation.
Infrared safety ensures that even in the presence of process-defining cuts, which are applied either to Φ 0 in the standard resummation or to Φ 1 in Geneva, these two approaches are identical in the limit T 0 → 0. Away from the limit, the two procedures will give the same result only if the quantities upon which the cuts are imposed are preserved by the Φ 1 → Φ 0 projection. For the Drell-Yan and V H production processes previously studied in Geneva, the only process-defining cut was applied on the invariant mass of the coloursinglet system. Since this quantity was preserved by every projection, the problem was avoided. In the case of diphoton production instead, the p T cuts on each photon are not preserved by our mappings. 7 Hence, these cuts applied to the projected Φ 0 configurations will effectively remove some contributions to the resummed and resummed-expanded terms of the cross section formula in eq. (4.11) which are present in the usual resummed results (i.e. without any recoil considered).
The difference between the two procedures is shown in the left plot of Fig. 2 for the resummed contribution alone, while the right plot shows the same comparison after matching to fixed order. We observe good compatibility between the two curves even at large T 0 , meaning that the difference between the two approaches is eliminated at FO by the matching procedure. We also observe larger fluctuations of the standard result at small T 0 values, likely due to combining the histogram bins of the matched calculation rather than combining the contributions on an event-by-event basis as is done in Geneva.
We are now in a position to compare the effects of the T 0 resummation matched to FO calculations by evaluating the cross section formulae in eq. (4.10) and eq. (4.11) at different accuracies. In the left panel of Fig. 3 we compare the NLL and NNLL results for the T 0 distribution in the peak region. In the same figure we also show the nonsingular contribution at NLO and NNLO in the same range of T 0 on the right.
In the peak region, the two results at different resummation accuracies do overlap, but we do not observe a substantial reduction of the resummation uncertainties. We also notice that the nonsingular contribution at very small T 0 values takes opposite signs at the different orders. Its size also looks particularly large when plotted on a linear scale, as is done in this plot (c.f. Fig. 1).
In Fig. 4 we instead plot the same resummed results matched to the appropriate FO calculation, in the peak, transition and tail regions. As a consequence of the size of the nonsingular corrections, the two curves only partially overlap for 1 < T 0 < 2 GeV, close to the peak. A similarly poor convergence was also observed for the p γγ T distribution after performing the q T resummation (see Ref. [36]). Also in the tail and transition regions, the effect of including the NNLO corrections is large and the uncertainty bands do not overlap with those at lower order. This was previously noticed in Refs. [20,22,26].

Subleading power corrections
In order to express the 0-jet cross section as in eq. (4.10), i.e. fully differential in the Φ 0 phase space and with NNLO 0 accuracy, one would need to implement a local NNLO subtraction method. This is not, however, explicitly needed in the Geneva approach, which is based on the N -jettiness subtraction [24,25].
Moreover, even if a local subtraction were provided, the predictions of an event generator would be inherently correct only for the total cross section and for observables which are left unchanged by the Φ 1 → Φ 0 and Φ 2 → Φ 1 → Φ 0 projections like, for example, the diphoton invariant mass. Hence, the presence of power corrections in T cut 0 cannot be avoided for generic observables that depend on the Φ 0 kinematics. We therefore replace the formula for the 0-jet cross section in eq. (4.10) with where the local subtraction and the expansion of the resummation formula are only needed up to O(α s ). This formula assumes that there is an exact cancellation between the FO and the resummed-expanded contribution at O(α 2 s ) below the T cut 0 . This holds for the singular contributions due to the NNLL accuracy of our resummation formula. However, this formula is only accurate at leading power in the SCET expansion parameter and fails to capture the nonsingular contributions in T cut 0 . These can be expressed as dσ nons ns , as a function of T cut 0 .
Since the functions f i (T cut 0 , Φ 0 ) contain at worst logarithmic divergences, the nonsingular cumulant vanishes in the limit T cut 0 → 0. In our calculation we include the term f 1 (T cut 0 , Φ 0 ) exactly by means of the NLO 1 FKS local subtraction. The f 2 (T cut 0 , Φ 0 ) term is instead completely neglected in eq. (4.21). This is acceptable as long as we choose T cut 0 to be very small. The effect of our approximation is shown in Fig. 5, where we plot the size of the neglected pure O(α 2 s ) terms in Σ ns (T cut 0 ) as a function of T cut 0 . The size of the missing contributions is not completely negligible and to reduce their impact we run with a default cut value of T cut 0 = 0.01 GeV. The magnitude of the missing corrections for such value of the cut is around 1.45 pb (which corresponds to ∼ 2% of the total cross section for the particular set of cuts chosen). Comparing this result to the previous Drell-Yan and V H calculations, we notice that in the diphoton case the relative size of the nonsingular corrections below the cut is larger.
One could improve on this by systematically calculating the subleading terms in the expansion parameter using a SCET formalism. Presently, only the first terms in the expansion are known, for a limited set of processes [85][86][87][88].
We eventually provide the missing nonsingular O(α 2 s ) contributions from an independent NNLO calculation obtained with Matrix [21], by simply rescaling the weights of the Φ 0 events by the cross section ratios. We remind the reader that the Matrix calculation is based on the q T subtraction method, which is very similar in spirit to 0-jettiness subtraction and thus in principle affected by the same issue. However, Matrix uses an extrapolation procedure for q cut T → 0 which provides an estimate of the cross section and its numerical error. Therefore, reweighting the Φ 0 events using the total cross section inputs from Matrix provides us with NNLO accuracy.

NNLO validation
After reweighting the Φ 0 events for the central scale choice as well as for its variations, we compare the inclusive (i.e. not probing additional radiation) distributions obtained with Geneva to the independent NNLO results obtained with Matrix [21]. This check is nontrivial since the complete dependence on the Φ 0 kinematics is not, in general, captured by the reweighting.
In Fig. 6 we show the transverse momentum of the hardest photon, the rapidity and the invariant mass of the diphoton system, and the absolute value of the cosine of the photon scattering angle in the frame of the LO partonic collision, defined as where ∆y γγ is the diphoton rapidity separation. After comparing three different choices of T cut 0 = {0.01, 0.1, 1} GeV in Geneva we conclude that the best agreement with the NNLO predictions for these inclusive distributions is obtained for T cut 0 = 0.01 GeV, as expected from the study of the missing power corrections (see sec. 4.4). We therefore set T cut 0 to this default value for all of our predictions.
As mentioned above, the Geneva predictions, despite being NNLO accurate, are not exactly equivalent to those of a NNLO calculation: indeed, they differ by power-suppressed terms as a consequence of the projective map which is used to define the Φ 0 events and by higher-order resummation effects that are not completely removed even after the inclusion of the additional terms in eq. (4.19). Nevertheless, the agreement between Geneva and Matrix for the distributions in Fig. 6 is good, within the uncertainty bands of the two calculations (representing the 3-point µ r and µ f variations).

Resumming the 1-jet/2-jet separation at LL
In order to provide an event generator which is as flexible as possible, and thus also able to provide exclusive predictions for higher-multiplicity bins, we proceed with the separation of the inclusive 1-jet cross section into an exclusive 1-jet cross section and an inclusive 2-jet cross section. We can achieve this separation by using the T 1 resolution variable. This introduces a new scale which in principle requires a simultaneous resummation of all the different ratios between the scales T 0 , T 1 and Q, in all the possible kinematic regions. A fully satisfactory treatment of this kind is still lacking, but, in the region T 1 T 0 , we can take the simpler approach explained next. We concentrate first on the T 1 resummation and start by separating where the terms dσ LL 1 and dσ LL ≥2 contain the LL resummation of the T cut 1 and T 1 dependencies respectively. The dσ match 1 and dσ match ≥2 terms contain the matching corrections to the required FO accuracy. At this point we can implement a unitary approach such that the resummed contributions take the form The evolution function 8 (or Sudakov factor) U 1 (Φ 1 , T cut 1 ) resums the T cut 1 dependence at LL accuracy in the region T cut 1 T 0 , while U 1 (Φ 1 , T 1 ) is the derivative of the evolution function with respect to T cut 1 . The latter resums the differential T 1 dependence and is evaluated at the projected configuration Φ 1 = Φ T 1 (Φ 2 ). The function P(Φ 2 ) is a normalised splitting probability which is defined similarly to P(Φ 1 ) in eq. (4.7). The quantity dσ C ≥1 /dΦ 1 , which appears both in eqs. (4.26) and (4.27), is the inclusive 1-jet cross section in the singular T 1 → 0 limit. Its NLO 1 expansion is given by where C 2 (Φ 2 ) reproduces the point-wise singular behaviour of B 2 (Φ 2 ) and acts as a local subtraction at NLO [89] with its own projection dΦ 2 ]. After requiring that dσ mc 1 and dσ mc ≥2 are accurate to NLO 1 and LO 2 respectively, the matching corrections are expressed as (4.30) 8 The explicit expressions for the evolution function U1(Φ1, T cut 1 ) up to NLL accuracy can be found in sec. 2 of Ref. [52].
After inserting eq. (4.26), eq. (4.27) and eq. (4.28) in the above equations and taking into account the appropriate phase space restrictions we find In the above expressions U 1 (Φ 1 , T cut 1 ) and U (1) 1 (Φ 1 , T 1 ) indicate the O(α s ) expansions of the evolution function U 1 (Φ 1 , T cut 1 ) and of its derivative U 1 (Φ 1 , T 1 ) respectively. Since the T 1 resummation is carried out to LL accuracy, the matching corrections still contain subleading single-logarithmic terms.
So far we have presented a NLO 1 +LL T 1 matched result, but we still need to incorporate the T 0 resummation that we discussed in the previous sections. We achieve this by requiring that the integral of the NLO 1 +LL T 1 result reproduces the T 0 -resummed result for the inclusive 1-jet MC cross section dσ mc The next step for the process at hand is nontrivial and requires some detailed explanation.
In the case of Drell-Yan or V H production, once this point was reached in the derivation one could simply proceed by summing the two contributions in the equation above. In this manner one could obtain a result which was independent of the evolution function, of its derivative and of the respective FO expansions, by exploiting the unitarity condition Unfortunately, in the case of diphoton production, this step is complicated by the presence of additional phase space restrictions on Φ 2 , due to the application of process-defining cuts, isolation cuts and projection cuts. Assuming that eq. (4.34) holds to to a good approximation despite the presence of all these cuts, we obtain (4.35) where the dots represent the remaining unitarity violating terms. We verified that, when implementing the T 1 resummation at LL, these terms are indeed numerically small both for the total cross section and at the differential level for the T 0 distribution. Therefore, we neglect them in the following. 9 By comparing the expression on the r.h.s. of eq. (4.35) with the r.h.s. of eq. (4.11) we obtain the result for dσ C ≥1 , Finally, we obtain the complete formulae for the exclusive 1-jet and the inclusive 2-jet cross sections as implemented in the Geneva code: ) could be computed numerically as where Θ PS (Φ2) indicates a set of phase space cuts. However, such a choice would have the drawback that the introduction of cuts would make the T1 resummation accuracy no longer clearly specified.
In order to simplify the notation, in the above equations we use the same symbol for two different projections. In eqs. (4.37) and (4.38) the symbol Φ 0 refers to the single projection Φ 1 → Φ 0 using the FKS mapping, while in eq. (4.40) the same symbol refers to the phase space point obtained after a double projection Φ 2 → Φ 1 → Φ 0 , where the first projection is evaluated with the T 0 preserving map and the second with the FKS map. The nonsingular contributions below T cut 0 and below T cut 1 are also explicitly written in eq. (4.37) for the exclusive 1-jet MC cross section and in eq. (4.39) for the 2-jet inclusive MC cross section respectively.

Parton shower and hadronisation
The Geneva interface to the parton shower has been extensively discussed in sec. 3 of Ref. [47]; here we summarise only the most relevant features. The shower makes the calculation fully differential at higher multiplicities by adding extra radiation to the exclusive 0and 1-jet cross sections and also further jets to the inclusive 2-jet cross section. The extra emissions are added in a recursive and unitary way. However, if additional analysis cuts are applied after the shower, large differences are usually expected also in distributions that are inclusive over the radiation.
For 0-jet events, the purpose of the shower is to restore the emissions which were integrated over when constructing the exclusive 0-jet cross section. In particular, the shower supplements events with additional emissions below the cut which are required to satisfy the constraint T 0 (Φ N ) < T cut 0 . In practice we allow for a small spillover at the level of 5%.
The showering of the 1-and 2-jet events requires a dedicated treatment to avoid significantly altering the T 0 spectrum calculated at the partonic level. The Φ 2 points generated after the first emission must satisfy the restriction T 1 (Φ 2 ) < T cut 1 as well as the projectability condition onto Φ 1 using the T 0 -preserving map T 0 (Φ 2 ) = T 0 (Φ 1 ) presented in sec. 4.1. In order to fulfil these constraints, we carry out the first emission in Geneva using the LL Sudakov factor described in sec. 4.6. The subsequent emissions are controlled by the shower and need only satisfy the constraint T 1 (Φ N ) < T cut 1 . In addition, we multiply the entire 1-jet cross section by a second Sudakov factor U 1 (T cut 1 , Λ 1 ) as described in Ref. [52]: The parameter Λ 1 determines the ultimate 1-jet resolution cutoff which we set to be much smaller than the Geneva cutoff parameter T cut 1 = 1 GeV. For the process at hand we use Λ 1 = 10 −4 GeV with the consequence that the 1-jet cross section is extremely suppressed and accounts for about 0.1% of the total cross section. The restrictions on the shower for Figure 7: Comparison of T 0 spectra between the partonic NNLO 0 +NNLL and the showered results, after interfacing to Pythia8, before the inclusion of non-perturbative effects (above). Comparison between the showered and hadronised T 0 spectra (below). The peak (left), transition (centre) and tail (right) regions are shown.
such a small contribution to the cross section can then be ignored. Regarding the showering of the partonic Φ 2 events, it was shown in Ref. [47] that the first emission of the shower acting on these events affects the T 0 distribution at order α 3 s /T 0 . From the above discussion it then follows that the showered events originate either from dσ MC 0 or dσ MC ≥2 . In the following we compare the partonic, showered and hadronised results for a selected set of distributions. We use the same inputs as above with the only difference that, for these comparisons, the fixed-order scale is set to µ FO = M T γγ . We use the parton shower program Pythia8 [51] interfaced to Geneva. At the hadronisation level, we switch off the hadron decays in order to keep the analysis simple and avoid contributions from secondary photons. For the same reason we also avoid including multi-parton interactions (MPI) and a QED shower.
We present in the upper panel of Fig. 7 the comparison between the partonic and showered results for the T 0 distribution, showing that the T 0 distribution is not modified by the shower above T cut 0 . Below T cut 0 , instead, the shape of the T 0 distribution is determined entirely by the shower; the effects are hardly visible since the cutoff is set to a very small value (T cut 0 = 10 −2 GeV).  We study the impact of hadronisation, which provides the nonperturbative effects, by comparing, in the bottom panel of Fig. 7, the showered and hadronised T 0 distributions. As expected, we notice a large difference between the two results only in the peak region, since the T 0 observable is very sensitive to additional low-energy hadronic emissions. At larger values of T 0 these corrections are instead suppressed as O(Λ QCD /Q) and their effects are lessened.
We further study effects due to the parton shower and hadronisation in Fig. 8, where we compare the partonic, showered and hadronised results for the transverse momentum of the photon pair, the transverse momentum of the hardest photon, the invariant mass of the photon pair and the pseudorapidity of the hardest photon. We first observe that for all the inclusive distributions, the NNLO accuracy is maintained at a very precise numerical level after both the showering and the hadronisation processes.
Moreover, although the transverse momentum of the photon pair, or any other exclusive observable, formally have the same logarithmic accuracy as the shower, we expect that it could also benefit from the high resummation accuracy of the T 0 distribution. We observe that the distribution is significantly modified after the shower only in the region below 10 GeV, while for larger values of p γγ T , the higher-order partonic result is practically recovered. In order to quantify the quality of our predictions for this observable we can compare with the direct resummation of p γγ T , which is performed in the Matrix+RadISH interface [90] up to N 3 LL p T +NNLO 0 accuracy.
In the left panel of Fig. 9 we show such a comparison at the partonic level, i.e. before the shower, observing a very good agreement. 10 In the the right panel of the same figure, we compare our results after showering but before hadronisation against the Ma-trix+RadISH results at both N 3 LL p T +NNLO 0 and NLL p T +NLO 0 accuracy. Even after adding the shower effects the Geneva results are in better agreement with those with higher logarithmic accuracy.

Inclusion of the gg channel contribution
The effects of including the gg channel contribution are quite large both for the total cross section (in the 6-10% range) and the differential distributions. This is a consequence of the relative size of the gluon parton distributions at the LHC.
In Fig. 10 we compare the results of Geneva with Matrix after the inclusion of the gg channel contribution for the same set of inclusive distributions presented in Fig. 6. As shown in the plots, we find very good agreement between the two calculations. We also show the effect of including the gg channel contributions by comparing to the Geneva results before its inclusion. Due to the numerical relevance of this channel, its NLO QCD corrections have been the subject of dedicated studies [15,17]. However, since these terms are formally of higher order (N 3 LO) with respect to the qq channel contribution, we neglect them in our calculation.
When showering events in the gluon fusion channel, we set the starting scale of the shower to be equal to the highest scale present in the process, which is the partonic centreof-mass energy. The reason for doing so is that we do not presently resum these contributions, whose resummation accuracy is then entirely given by the shower. A dedicated higher-accuracy resummation of this channel is of course possible but is left to future investigation.
In Fig. 11 we show the comparison between the partonic, showered and hadronised results after the inclusion of this channel for the T 0 distribution, the rapidity of the diphoton system, the transverse momentum of the photon pair and the transverse momentum of the hardest photon. We observe somewhat larger effects after the inclusion of the shower compared to the case of the qq channel alone, especially for the T 0 , p γγ T and p γ h T distributions. The y γγ distribution is instead left untouched by the shower. These are most likely due to the high scale at which we start the showering process. A similar behaviour was also observed for the V H production process in Geneva after including the gg channel contribution, as well as in the Powheg and MC@NLO implementations of similar processes [91,92].

Event generation and analysis cuts
In this subsection we study the effects of applying process-defining and isolation cuts at the generation and analysis levels, both before and after shower and hadronisation. At the generation level, we are forced to use a smooth-cone isolation procedure in order to generate 10 We compare against results at N 3 LL because the public version of Matrix+RadISH does not presently allow for NNLL accuracy. In order to have a like-for-like comparison with Geneva results we have also selected an additive scheme for the matching of the resummation to the fixed-order.   well-defined, IR-finite events, without fragmentation contributions. At the analysis level, however, when one is interested in comparing with data, a fixed-cone isolation algorithm is needed. For these reasons, in sec. 5 we will apply a hybrid isolation procedure, i.e. first imposing a very loose smooth-cone isolation cut at the generation level followed by a tighter fixed-cone isolation at the analysis level.
In order to check the consistency of this approach, we must first quantify the dependence of the results at the various levels of the analysis from the cuts imposed at generation. We separate this investigation into two parts: in the first, at parton level, we examine the power-suppressed isolation effects due to the phase-space projections below the jet resolution cutoffs; in the second, after the shower, we study the effect of the random momenta reshuffling due to recoil and hadronisation.
For the first part, we use the set of "tight" cuts introduced in eqs. (4.1) and (4.2), which we report here for convenience  We first generate the events by applying the set of loose cuts in eq. (4.44) and, as a second step, we analyse them by applying the tighter cuts of eq. (4.43) before showering. We compare these predictions to the results obtained by directly applying the set of tight cuts at generation level. This is shown in Fig. 12 for the pseudorapidity of the softer photon and the T 0 distribution, where we show the results of the calculation directly carried out with tight generation cuts together with that where we apply loose generation cuts (as in eq. (4.44)) and tighter cuts at the analysis level. The two predictions are in good agreement and this gives us confidence that our results are not strongly dependent on the generation cuts applied.
For the second part, one should expect that power-suppressed effects connected with the recoil after any emission could modify the momenta of the final-state particles and, consequently, result in a different rate of events passing the analysis cuts compared to those passing the generation cuts. This effect is particularly severe after the shower, since multiple emissions can greatly reshuffle the final-state momenta. The same applies to the reshuffle used by SMC programs to impose momentum conservation after hadronisation.
In order to quantify these effects we compare in Fig. 13 results obtained employing the loose generation cuts in eq. The figure shows reasonable agreement between the two predictions for the transverse momentum of the photon pair and the cosine of the photon angle in the Collins-Soper frame, demonstrating that the size of these effects is not large for variations of the isolation radius at generation level. However, qualitatively we did find a stronger dependence of the final results on the choice of the generation cuts on the photons' transverse momenta.

Results and comparison to LHC data
In this section we compare our predictions against 7 TeV LHC data obtained from both ATLAS [6] and CMS [10]. We employ the hybrid isolation procedure, as detailed in sec. 2 and sec. 4.9. This means that we first generate partonic events with looser smooth-cone isolation cuts, and only after the shower and hadronisation procedures do we apply the tighter analysis cuts and fixed-cone isolation algorithms which are used by the ATLAS and CMS experiments.
For these particular comparisons, we generate events using the NNPDF31 nnlo as 0118 PDF set [93]. We set the FO scale to µ FO = M T γγ and apply the following process-defining cuts at generation level: GeV, R iso = 0.1, and n = 1 . (5.1) Note that, in principle, there is no need to require a lower limit on the invariant mass of the photon pair, but, since our hard function is evaluated at µ H = M γγ in the resummation region, we set this lower cutoff so that α s (µ H ) is not evaluated at scales which are too small. We then shower and hadronise the partonic events generated with the cuts in eq. (5.1). To this end, we use Pythia8 and, in order to avoid any contamination by photons coming from hadronic jets, we prevent the decay of hadrons. For this comparison we also include the MPI and the QED shower effects. However, we do not allow QED splittings of photons into quarks or leptons. 11 In Fig. 14 we show the comparison between our predictions and the ATLAS data at 7 TeV [6] for the invariant mass of the photon pair M γγ , the transverse momentum of the diphoton system p γγ T , the azimuthal-angle separation between the two photons ∆φ γγ and the cosine of the polar angle θ * γγ in the Collins-Soper frame of the diphoton system. The comparison is carried out by applying to the showered and hadronised events the    where ∆R γγ = ∆η γγ 2 + ∆φ γγ 2 is the separation between the photons. Overall, we observe very good agreement between the theoretical predictions and the ATLAS data. For the invariant mass distribution, in the region M γγ ≥ 350 GeV, above the tt production threshold, the theoretical predictions seem to depart from data. Here we expect that the inclusion of EW corrections and of the two-loop diagrams with a closed top-quark loop in the hard function will improve the theoretical description. We also observe a deviation in the extreme region ∆φ γγ ∼ π of the ∆φ γγ distribution. Next, we compare against the CMS data at 7 TeV [10]. The CMS analysis uses a fixed-cone isolation algorithm with parameters The comparison between the predictions obtained with Geneva + Pythia8 and the CMS data is shown in Fig. 15 for the same set of distributions presented for the ATLAS case. Unfortunately, since no corresponding Rivet analysis is available, we have implemented the aforementioned cuts in the Geneva analyser. The M γγ distribution is shown with a linear scale on the abscissa up to 100 GeV and a logarithmic scale beyond. Similarly, the p γγ T distribution is shown with a linear scale up to 10 GeV and a logarithmic scale beyond. We observe a similarly good agreement for the inclusive distributions as for ATLAS. Close to ∆φ γγ ∼ π the photons' azimuthal separation shows an opposite trend compared to that of ATLAS, but our predictions in this case always agree with CMS data within experimental errors.
For both experiments, the p T,γγ prediction also shows some systematic trends, undershooting the data at the very low end of the spectrum. However, the theoretical uncertainties on our predictions for this observable are not completely exhaustive as they do not yet include e.g. the uncertainties related to the matching to the shower or to the variations of the shower parameters. Since this discrepancy seems entirely caused by shower effects (cfr. Figs. 8 and 9) we expect that including shower uncertainties and hadronisation tuning will improve the agreement with data.

Conclusions
We have presented the first calculation for the production of isolated photon pairs at the LHC resummed in the 0-jettiness resolution variable to NNLL accuracy and matched to the NNLO calculation. This has been performed within the Geneva Monte Carlo framework, allowing us to interface to the Pythia8 parton shower and hadronisation model. Our work constitutes the first NNLO event generator matched to a parton shower (NNLO+PS) for this process.
The implementation of photon pair production in an event generator is complicated by the nontrivial process definition, which suffers from QED singularities. In order to solve this problem we used a smooth-cone isolation algorithm in Geneva to remove such divergences. We studied the dependence of our results on the parameters of the isolation applied at generation level and also investigated the differences due to the recoil between the standard resummation approach and our implementation in Geneva. We have validated our calculation using the Matrix program to NNLO accuracy and found good agreement when using T cut 0 = 0.01 GeV, which sufficiently reduces the size of the neglected subleading power corrections.
We have further studied the effects of parton shower and hadronisation. We first ensured that the T 0 distribution is not affected by the shower alone if no additional phase space cuts are applied after showering. We have then quantified the size of the nonperturbative effects provided by hadronisation at the low end of the T 0 spectrum, finding the expected modifications.
We have also found that the NNLO description of the inclusive distributions is preserved by the shower and hadronisation procedures. Larger effects instead appear for more exclusive distributions. After including the gg channel at O(α 2 s ), we observed larger shower effects, both for inclusive and exclusive distributions, connected to the usage of a higher starting scale of the shower for this contribution.
Finally, after including both MPI and QED effects in Pythia8, we have used a hybrid isolation procedure to compare to LHC data at 7 TeV. We find in general a good agreement with both ATLAS and CMS, with minor tensions appearing for exclusive distributions such as p T,γγ and ∆φ γγ . These could possibly be ameliorated after the inclusion of theoretical uncertainties connected to the matching to the shower or to the variations of the shower and hadronisation parameters.
Possible directions for future work include the NLL resummation of the gluon fusion channel contribution, which is currently included only at LO+PS, the addition of electroweak corrections, and of the top-quark mass effects in the two-loop hard function. Moreover, other interesting processes with a single photon in the final state, such as Zγ and W γ, could be also implemented in the Geneva framework.
The code used for the simulations presented in this work is available upon request from the authors and will be made public in a future release of Geneva.

A Hard functions for γγ production
The hard function is one of the ingredients of the SCET factorisation formula in eq. (3.3), and, in order to achieve NNLL accuracy, it is needed up to O(α 2 s ) in perturbation theory. It can be calculated as the square of the hard interaction matching coefficients. The relevant partonic process is given by and the hard function can be extracted from the 2 → 2 matrix elements, which were computed up to two-loop level in Ref. [66]. We introduce the following Mandelstam invariants: Momentum conservation implies that they satisfy the relation s+t+u = 0, where s > 0 and t, u < 0. Hence only two of the invariants are independent and it is convenient to express the results in terms of s and the dimensionless parameter x = −t/s. The UV-renormalised amplitudes have the following perturbative expansion 12 |M qqγγ ( , s, x) = 4πα |M where = (4 − d)/2 is the dimensional regulator. Note that the perturbative coefficients on the r.h.s. depend also on the renormalisation scale µ r while the all-order amplitude on the l.h.s. is independent of it. After UV renormalisation the virtual amplitudes still contain IR poles in the dimensional regulator . 13 We subtract these poles in the MS scheme by acting on the amplitudes with a renormalisation factor Z, where the IR-finite amplitudes now depend on the renormalisation scale µ. The explicit form of the renormalisation factor Z was recently determined up to O(α 4 s ) in massless QCD in Ref. [96] by investigating the structure of the associated anomalous dimension. For our 12 In the SCET literature it is customary to express the perturbative expansions of the hard matching coefficients of certain operators in powers of αs/4π, rather than αs/2π. This introduces a 2 L conversion factor for L-loop amplitudes. 13 The precise structure of the IR poles depends on the regularisation scheme employed in the calculation.
For the relation between different regularisation schemes, see the analysis done in Ref. [95]. where for simplicity we dropped the common kinematic dependence on s, x from all terms in the above equations. The renormalised amplitudes on the l.h.s. are free of IR poles. For the diphoton production process the Z factors are defined in terms of anomalous dimension coefficients where β 0 = 11/3C A − 4/3T F n f , and for this process The perturbative coefficients of the anomalous dimensions γ cusp and γ q entering in the above equations are given by We briefly comment on these equations. It is important to notice that in the above equations the pole content of Z (2) and I (2) is different, but when all the other terms are taken into account, the subtraction removes the same poles and the result is finite. While the Z( ) factors are defined in the MS scheme and only contain poles, the I operators need to be expanded to the correct order in . For example, the term between parentheses in the second term of the second line of eq. (A.19) is finite in but must be expanded to O( 2 ) since it multiplies I (1) , which contains up to a double pole in . The squared amplitudes (summed over colours and spins) are defined as where we dropped the dependence on , s, x on the r.h.s. terms. The squared amplitude coefficients are not yet averaged over the initial spin polarisations and colours, which introduces a factor 1/4N 2 c . The squared amplitudes should also be divided by the number of identical particles in the final state, and we therefore need to multiply by an additional factor 1/2. The explicit analytic expressions in eq. (A.22), in terms of logarithms and classical polylogarithms can be found in Ref. [66]. The IR pole terms are expressed in terms of Catani's I (1) and I (2) operators. The translation into the MS Z-factors can be found in Ref. [97] for a generic QCD process with coloured particles in the final state. For our convenience one of the authors of Ref. [66] directly provided us with the explicit expressions in eq. (A.22) in the FORM format.
With this in mind we can now define the hard function for the diphoton production process as where all of the hard function coefficients on the r.h.s. are real and finite. We report here the explicit expressions for the first two perturbative coefficients of the hard function: H (1) (s, x, µ) = (N 2 c − 1) Q 4 q 1 x + 1 1 − x − 2 − ln 2 s/µ 2 + 3 ln s/µ 2 − 7 + 7π 2 6 where Q q is the electric charge of the active quark. Unfortunately, the NNLO hard function is too lengthy to be included here. The result is available upon request to the authors.