NNLO+PS Monte Carlo simulation of photon pair production with MiNNLOPS

We present a NNLO QCD accurate event generator for direct photon pair production at hadron colliders, based on the MiNNLOPS formalism, within the POWHEG BOX RES framework. Despite the presence of the photons requires the use of isolation criteria, our generator is built such that no technical cuts are needed at any stage of the event generation. Therefore, our predictions can be used to simulate kinematic distributions with arbitrary fiducial cuts. Furthermore, we describe a few modifications of the MiNNLOPS formalism in order to allow for a setting of the renormalization and factorization scales more similar to that of a fixed-order computation, thus reducing the numerical impact of higher-order terms beyond the nominal accuracy. Finally, we show several phenomenological distributions of physical interest obtained by showering the generated events with Pythia8, and we compare them with the 13 TeV data from the ATLAS Collaboration.


Introduction
The production of two isolated photons at hadron colliders, henceforth denoted diphoton production and abbreviated as pp → γγ, is one among the most relevant Standard Model (SM) processes, due, on one side, to the high production rate, and, on the other, to the relative cleanliness of the experimental final state. Furthermore, diphoton production is the dominant background for studies involving Higgs boson production decaying into a photon pair, and, notably, it was one of the dominant backgrounds for the Higgs boson of the event generation, we have also devised a new mapping from the pp → γγ + j to the pp → γγ kinematics, that preserves the direction of one photon with respect to the beam axis, thereby allowing for a full control of the singular regions. The combined use of the new mapping and of the method to deal with the QED divergences allowed us to simulate diphoton production at NNLO+PS accuracy without introducing any generation or technical cut, as done instead in refs. [35,68]. This gives us the advantage that, on one side, we do not have to worry about checking the independence of the fiducial physical differential cross sections from the generation and technical cuts, and, on the other side, we can use the same set of generated events for any choice of fiducial cuts.
In this paper we also propose a generalization of the MiNNLO PS method that allows for greater flexibility in the choice of the renormalization and factorization scale used in the evaluation of the non-singular term for F + 1 jet production, where F is the color-singlet system. In the original formulation of the MiNNLO PS method, such term is evaluated with scales of the order of the transverse momentum of F . The prescription we introduce here generalizes such aspect of the MiNNLO PS method, allowing one to treat this term more similarly to a fixed-order computation. Such choice turned out to be desirable in γγ production, for comparisons with fixed-order results, where this term gives an important contribution to the total cross section, unlike in previous studies, such as those in refs. [50,51], where it was small.
The paper is organized as follows: in section 2 we describe the ingredients used to build the event generator and the details of the novel aspects we introduce here for the first time, whose validation is discussed in section 3. In section 4 we show some phenomenological results, and we compare our predictions with the ATLAS diphoton data [74]. Finally we give our conclusions and outlook in section 5.

Outline of the calculation
In this section we review the theoretical formulation of our calculation. We first introduce the basic notation used throughout the paper and describe the main theoretical issues to be dealt with in the implementation of a NNLO+PS event generator for photon pair production. Then we discuss in detail the handling of the QED singularities and introduce the mapping used to project the Φ γγj kinematics onto the Φ γγ one. We conclude the section with a description of the modifications we have introduced in the MiNNLO PS formalism in order to reproduce more accurately the results of a fixed-order calculation for distributions totally inclusive with respect to the QCD radiation.

Description of the process
We consider the process of direct production of a photon pair in proton-proton scattering pp → γγ + X, (2.1) with the requirement that the two photons are isolated, i.e. each photon has a minimum transverse momentum and is isolated with respect to the other photon and to the final-state hadronic activity. These requirements are needed to make the process well-defined both from the theoretical and the experimental point of view. We call p γ 1 and p γ 2 the momenta of the hardest and next-to-hardest photons and denote the squared mass and transverse momentum of the diphoton system as The starting point for the MiNNLO PS method we use in this work is the NLO differential cross section for diphoton plus one jet production pp → γγ + j . (2.3) The matrix elements for this process at NLO in QCD have been obtained from the automated interface of the Powheg Box Res [75] to OpenLoops2 [76][77][78][79][80].
Since the goal of the MiNNLO PS formalism is to reach NNLO accuracy for inclusive observables in the diphoton system, one also needs the two-loop amplitudes for qq → γγ, that have been taken from refs. [81,82] and implemented into the code.
In this paper, we work in the approximation of five light quarks and neglect the contributions given by two-loop diagrams with the massive top quark. We consider instead the contributions given by the massive top quark in the one-loop diagrams for the process pp → γγ + j.

Handling of the QED divergences
In this section, we describe a general way to deal with processes that present QED divergences at the Born level within the Powheg formalism. 7 Before illustrating the method, we briefly review the relevant parts of the Powheg Box framework.
We start by introducing the general form of the Powheg [43] differential cross section for a 2 → n Born process, using the notation of refs. [84,85]. Indicating with dΦ n the phase space for the n-body final state, we write the (n + 1)-body phase space dΦ n+1 in terms of dΦ n and three additional radiation variables, that we collectively label Φ rad dΦ n+1 = dΦ n dΦ rad . (2.4) We indicate with B(Φ n ), V (Φ n ) and R(Φ n+1 ) the Born, virtual and real amplitudes, convoluted with the corresponding parton distribution functions (PDFs) and multiplied by the flux factor, and we split R into two positive contributions: R s (Φ n+1 ), that contains all the QCD singularities, and a QCD-finite term, In general we achieve this separation with the help of a suitable function F (Φ n+1 ), such that The way the Powheg Box deals with QCD divergences at Born level has already been illustrated in several papers, e.g. ref. [83]. 8 Please notice that R f (Φn+1) could well be set to zero, in which case Rs(Φn+1) = R(Φn+1), the whole real contribution.
-5 - The detailed form for F (Φ n+1 ) will be discussed in section 2.2.1. When R(Φ n+1 ) is split as in eq. (2.5), the Powheg differential cross section can be written as where k T is the transverse momentum of the Powheg jet, 8) and the Powheg Sudakov form factor contains only the QCD singular part of the real contributions R s (Φ n+1 ). In the expressions above k min T acts as an infrared cutoff for unresolved radiation. If the Born amplitude is divergent, the Powheg Box applies a suppression factor S B (Φ n ) toB(Φ n ) such that the product S B (Φ n )B(Φ n ) is integrable over the entire Born phase space, and the Φ n kinematics can be generated from the distribution S B (Φ n )B(Φ n ) without the need of any hard cut. At the end, each event is given a weight 1/S B (Φ n ) in order to compensate for the suppression factor.
In general, this way of regularizing divergences by means of a suppression factor that depends on Φ n can be used every time QCD and/or QED singularities are present at the Born level. While applying a suppression factor onB(Φ n ) is enough when only QCD singularities are present, for diphoton production we also have to deal with the QED singularities appearing inside R f (Φ n+1 ). We cannot simply suppress them through a second suppression factor (that would be a function of Φ n+1 ) since this term is integrated in the radiation variables insideB(Φ n ), preventing one to compensate for such suppression factor, after the events have been generated.
We then proceed as follows: we take advantage of the possibility to separate the real contributions into two terms (see eq. (2.5)) in such a way that R s (Φ n+1 ) contains all the QCD singularities but no QED ones. Since the generation of Φ n+1 according to R f (Φ n+1 ) in eq. (2.7) is performed with a hit-and-miss technique, we apply a QED suppression factor that is now integrable. Finally, we multiply the weights of the generated events by the factor 1/S R (Φ n+1 ). 9

The damping function F
According to the FKS method [86,87], the real contribution R is partitioned into a sum of terms R αr , each of them having at most one collinear and one soft singularity associated with one parton R = αr R αr . (2.10) In each α r region, the radiated parton (r) can then become soft or be collinear to an emitting one (e), and we can define the damping function Tj if i is an initial-state particle, if i and j are final-state particles, (2.12) and the sum in the denominator runs over the n c massless charged particles and the n γ photons. In eq. (2.12), p Tj and E j are the transverse momentum and energy of the particle j, and θ ij the angle between the particles i and j, computed in the partonic center-of-mass frame. 10 Every contribution R αr to the real amplitude is then split into two terms (2.14) In the limit where the radiated parton is soft or collinear to the emitter d e,r is small, and F αr tends to 1, so that all the QCD singularities are in R αr s . When instead the photon γ j becomes soft or collinear to a massless charged particle c i , the term d c i ,γ j is small, and F αr tends to 0, so that R αr s is free from QED singularities.

The suppression factors S B and S R
We have chosen as suppression factor S B (Φ n ) introduced in section 2.2 the following expression where p Ti is the transverse momentum of the particle i with respect to the beam axis, R ij the angular distance between the particles i and j in the azimuth-pseudorapidity plane 16) and the barred quantities are arbitrary constant parameters needed to define the transition point (from zero to one) of each factor in the product. The first three term of eq. (2.15) suppress the limit where the final-state parton or the two photons are soft or collinear to the beam axis, while the last two terms suppress the regions where the photons are collinear to the final-state parton. The power a does not need to be the same for all the terms, but we have chosen the common value a = 1 for simplicity. The expression for S R (Φ n+1 ) can be obtained by generalizing eq. (2.15), for the case with an extra jet. When the MiNNLO PS formalism is used, the QCD initial-state singularity is already regularized by the Sudakov form factor from the resummation formalism, and the first term in eq. (2.15) is no longer needed. We give more details in section 2.4.1.

The Φ γγj → Φ γγ mapping
When applying the MiNNLO PS formalism, we need a mapping from the γγj phase space to the γγ one. In fact, as it will be recalled in the next section, the Sudakov form factor (eq. (2.19)) and the D term (eq. (2.25)), at the core of the MiNNLO PS method, are functions of Φ γγ , as they contain the qq → γγ amplitudes. One needs then to evaluate these quantities while integrating over the Φ γγj phase space upon which the PowhegB function depends, thereby requiring a smooth Φ γγj → Φ γγ mapping.
The qq → γγ amplitudes are singular when the photons are collinear to the beam axis. From the physics point of view, one does not expect to evaluate such amplitudes arbitrarily close to their singularities, as such phase space points are discarded by the request of the presence of two isolated photons in the final state. Since in the MiNNLO PS formalism the Φ γγ kinematics is obtained by a mapping from the Φ γγj one, we need to make sure that we never get too close to the singular regions, and without having to introduce an explicit cut on the transverse momentum of the single photons in the Φ γγ kinematics.
We have built such a mapping without having to introduce any cutoff or isolation criteria, at any stage of the event generation, as done in other Monte Carlo simulations dealing with photons [35,68]. The mapping that we propose is such that the direction of one of the photons with respect to the beam axis in the laboratory frame, for a given phase space point Φ γγj , is preserved in the projected point Φ γγ . As a consequence, a point in Φ γγ with small p T always comes from a projection of a point in Φ γγj where at least one photon is close to the beam axis, and this configuration is suppressed by the factor S B of eq. (2.15). A detailed derivation of such a mapping is given in appendix A.

MiNNLO PS differential cross section
In this section, after briefly recalling the general method, we discuss the modifications that we have introduced in the definition of the MiNNLO PS differential cross section given in refs. [50,51]. Following the conventions introduced in those papers, we write the p T spectrum of the cross section for diphoton production as dσ dΦ γγ dp T = dσ s dΦ γγ dp T + dσ f dΦ γγ dp T , (2.17) where σ s (also called singular contribution in the following) is obtained from the resummation of logarithmically enhanced contributions at small-p T , while σ f (also called non-singular term) is the difference between the fixed-order differential cross section and the truncated perturbative α S expansion of σ s , till the second order. The singular contribution can be written as dσ s dΦ γγ dp T = d dp where the Sudakov form factor is given bỹ and L is a function of the PDFs, of the qq → γγ Born, one-and two-loop matrix elements, and of the collinear coefficient functions (see ref. [88]). The functions A andB in eq. (2.19) can be written as where H (1) is the ratio between the one-loop and the Born qq → γγ amplitudes in the MS subtraction scheme, and introduces a dependence on the phase-space kinematics inB (2) . The coefficients A (1) , A (2) , A (3) , B (1) , and B (2) for quark-initiated processes are collected, for example, in ref. [50]. More details can be found in appendix B. The non-singular contribution dσ f is instead given by the difference where dσ NLO γγj is the NLO differential cross section for γγj production, and [dσ s ] (i) is the i-th order of the expansion of dσ s in the strong coupling constant. In this paper, we use the notation [X] (i) for the coefficient of the i-th term in the perturbative expansion of the quantity X. The above difference is integrable, since the expansion of σ s cancels the non-integrable terms of dσ NLO γγj in the p T → 0 limit. At variance with refs. [50,51,68,69,72], we have chosen to set the renormalization and factorization scales µ R and µ F in dσ f to Q, instead of p T . This is why we have added the subscript Q to the quantities appearing in eq. (2.23). While formally any scale choice would be legitimate for evaluating this term, since the difference would be of order O(α 3 S ), beyond the accuracy of our calculations, setting the scale to Q allows to follow more closely what is typically adopted in fixed-order calculations, thereby allowing for a more accurate comparison with the latter. Moreover, the choice of the central scale also plays a role in the estimation of the theoretical uncertainties. In this case, since the size of the non-singular contribution is numerically relevant, setting the scale to p T would result in an unnatural enhancement of the scale-variation bands. We will give more quantitative details in section 3.2, where we compare the new µ = Q and the µ = p T scale choices for dσ f for the previously discussed cases of Drell-Yan and Higgs boson production. Following refs. [50,51], the singular contribution of eq. (2.18) can be rewritten as where all the terms are evaluated at µ R = µ F = Q.
The D(Φ γγ , p T ) term in eq. (2.24) has a formal expansion and, at difference with eqs. (2.26) and (2.27), all the terms are evaluated at µ R = µ F = p T . 11 The explicit expression of these terms are collected in appendix B, together with the expressions of the other ingredients needed in the MiNNLO PS formulae. We would like to stress that the cancellation of the divergences associated with the small p T limit in eq. (2.23) is numerically challenging. For this reason, and following the MiNLO original approach, we have chosen to multiply dσ s by a Sudakov form factor, after adding an additional term to preserve the nominal α 2 S accuracy . (2.29) 11 We stress the fact that we do not compute separately the [D(Φγγ, pT)] (i) terms, but we compute numerically the whole function D, as discussed, for the first time, in ref. [51].

S(p
We stress that the two approaches are equivalent up to O α 2 S . Summarizing, the MiNNLO PS differential cross section for pp → γγ + X production can be written as where we have introduced the symbol dσ LO γγj /dΦ γγj to denote the LO differential cross section for pp → γγj production, and the function F corr (Φ γγj ), needed to spread the contributions proportional to the D(Φ γγ , p T ) terms over the entire Φ γγj phase space. It has the property that, given an arbitrary function G(Φ γγ , p T ), For the explicit expression we have used for F corr (Φ γγj ) we refer to ref. [50]. Equation (2.31) is the main result of this section and it summarizes the novel aspects we introduced to the MiNNLO PS method, as can be evinced by comparing it against, e.g., eq. (3.7) of ref. [50]. In the rest of the manuscript, we will denote the MiNNLO PS results obtained according to eq. (2.31) with the acronym FOatQ. We also recall that the mapping introduced in section 2.3 guarantees that all the Φ γγ -dependent terms in eq. (2.31) are evaluated in kinematic points away from the singular regions.
By rewriting the NLO differential cross section for γγj production in the following compact form and introducing the MiNNLO PS improvedB function our final expression for the differential cross section dσ γγ reads

The suppression factors with MiNNLO PS
As already discussed in section 2.2.2, we do not need to suppress the small p T region of the first jet while the Powheg Box Res integrates theB(Φ γγj ) MiNNLO PS function over the whole phase space, since the presence of the MiNNLO PS Sudakov form factor suppresses such region. In this case, the Born suppression factor S B in eq. (2.15) can be replaced by (2.36)

Scale settings in the small-p T limit and modified logarithms
We freeze the renormalization and factorization scales at values below Q 0 = 2 GeV to avoid the Landau pole and further non-perturbative effects, connected with the PDF evolution to lower scales. We stress that this scale does not act as a cutoff in the integration over the physical space but only enters in the evaluation of the singular contribution in eq. (2.24), since all the other terms in our formulation of the MiNNLO PS differential cross section are evaluated at the hard scale Q. We also highlight that, due to the presence of an overall Sudakov form factor, the dependence of the differential cross section on Q 0 is strongly suppressed. 12 In addition, following what was done in ref. [50,51], we smoothly turn off the contribution of the logarithms in the D functions at large transverse momentum with the replacement so that the p T → 0 limit remains unaffected, while at p T > Q, the modified logarithm tends to zero. In our simulation we have set p = 6. 13

Validation
In this section we discuss the validation of our MiNNLO PS predictions. We first briefly present the setting of the physical and technical parameters of the calculation and the isolation criterion used to define the diphoton process. We then study the effects that the modifications to the MiNNLO PS differential cross section described in section 2.4 have on two previous implementations of the method (i.e. the Drell-Yan and Higgs production processes). Finally, we present a validation of the MiNNLO PS results for diphoton production against the fixed-order distributions produced with the public version of the Matrix code.

Physical and technical parameters and photon isolation criterion
The phenomenological results presented in this paper have been obtained for a protonproton collider with a hadronic center-of-mass energy of 13 TeV. We have used the Lhapdf [89] parton distribution function set NNPDF31_nnlo_as_0118 and the evolution of α S provided by the same package. The electromagnetic coupling for the final-state photons has been set to α = 1/137, and the mass of the top quark to m top = 173.2 GeV. The fixed-order results have been obtained with the public version of the Matrix code [11,77,81,88,[90][91][92], setting the central renormalization and factorization scales equal to the invariant mass of the diphoton system µ R = µ F = Q. The uncertainty band has been estimated via a seven-point scale variation obtained by multiplying and dividing the central renormalization and factorization scales by a factor of 2.
We apply the photon isolation prescription of ref. [9] to the two final-state photons. For each photon, we compute the angular distance R iγ with respect to the i-th final-state parton. We discard the event unless, for every photon and every R < R c , where n i is the number of final-state partons, p i T is the transverse momentum of i with respect to the beam, and In our analysis, we have set In addition, the two photons have to fulfill where p Tγ 1 and p Tγ 2 are the transverse momenta of the hardest and next-to-hardest photons, m γγ is the diphoton mass, and The values of the barred quantities in eq. (2.36) and of the power a have been set tō In addition, the Powheg Box Res code has been run with the flag doublefsr set to 1. This parameter was first introduced in ref. [93], and we refer the interested reader to this paper for further details.

Validation for Drell-Yan and Higgs boson production
In this section, we compare the results obtained with the original MiNNLO PS method used in ref. [51] for Drell-Yan and Higgs boson production against those obtained with the new formulation spelled in section 2.4. As discussed in that section, one expects the new formulation to yield results compatible with the original method for processes where the size of non-singular corrections is small compared to the total cross section. An estimate of the size of the non-singular correction can be obtained by comparing the total NNLO cross section against the integral of the first term on the right-hand side of eq. (2.31) over the full phase space. The results obtained for diphoton, Drell-Yan and Higgs boson production are the following For diphoton production, σ γγ s contributes to only about one third of the total cross section, at difference with Drell-Yan and Higgs boson production, thereby justifying the choices made in this paper. In the left pane we show the rapidity distribution of the dilepton system in Drell-Yan production, whereas in the right pane we show the Higgs boson rapidity. The ratios with respect to the Matrix results are also shown in the lower panel.
As a further validation, in figure 1 we compare the rapidity distribution of the color singlet for Drell-Yan and Higgs boson production. We show the distributions obtained with the original formulation, where the finite terms are evaluated at µ = p T (labelled as no FOatQ), and the new one, where such terms are evaluated at µ = Q, with Q the invariant mass of the color singlet (labelled as FOatQ). The MiNNLO PS results shown in the figure are obtained after generating the Powheg hardest emission, i.e. at the "Les Houches" level.
In the plots we also show the NNLO results from Matrix, obtained setting µ = Q. For this comparison only, we use the same PDF sets as those used in ref. [51]. In the lower inserts of the figure we plot the ratio of the displayed distributions with respect to the Matrix result.
The curves show a very good agreement between the NNLO and the MiNNLO PS results obtained with either formulations, both for the central scale and the uncertainty band. In particular, since the Drell-Yan process features a very small perturbative uncertainty band, the remarkable agreement between the NNLO and FOatQ curves displayed in the left pane of figure 1 is a robust validation of the new formulation of the MiNNLO PS method.

Diphoton production: comparison with Matrix
In this section we validate the differential cross section of eq. (2.31) against the fixed-order NNLO one implemented within the public version of the Matrix code. We expect the two results to agree up to terms beyond the NNLO accuracy. The Matrix results have been obtained setting the slicing parameter r cut = 0.0005 and for a central scale choice µ R = µ F = Q. Using the isolation criteria and the fiducial cuts reported in section 3.1, the total cross section computed by Matrix and MiNNLO PS are in agreement within the statistical errors, and they read respectively In figure 2 we plot the distributions for the invariant mass and rapidity of the diphoton system computed with Matrix and MiNNLO PS , in the top panes, and the ratio of the two curves in the lower ones. The figure shows a reasonably good agreement between the two curves within the statistical errors. We ascribe the mild difference in the high invariantmass region to effects beyond the nominal accuracy of our result, which can differ from the purely fixed-order Matrix one by higher-order effects, present in eq. (2.31). We have explicitly verified that, by applying a transverse momentum cut on the diphoton system, the same trend is also present when we compare the exact fixed-order NLO result for γγj production against the MiNNLO PS one: when we move away from the region of small transverse momentum of the diphoton system that is affected by large logarithms, we still observe a mild difference between the central values of the two distributions, but still within the scale-variation bands.    14 Matrix provides also extrapolated values for the total cross section down to rcut = 0, that, for the central scale, is equal to σ Matrix tot = 153.9 ± 1.9 pb. We notice that, with the extrapolation process, the associated statistical error is larger than in the non-extrapolated one. For comparison with eq. (3.9), the extrapolated scale variation band is ±4%. and are in good agreement. In figures 3 and 4 we compare the scale variation bands obtained from the two codes for the invariant mass and rapidity distributions of the diphoton system and the transverse momentum of the hardest and next-to-hardest photon. We find an overall good agreement among the different curves, and compatible size for the scalevariation bands.

Phenomenological results
After validating the implementation of the MiNNLO PS differential cross section for diphoton production given in eq. (2.31), in this section we present some distributions of physical interest obtained from the generated events, before and after passing them through a parton shower. We have generated about 16 million events without any generation cut, except for a minimum invariant mass of the diphoton system of 10 GeV. 15 As stressed previously, except for this last constraint on the invariant mass, no other cuts have been imposed, so that these events can be used to make predictions with arbitrary fiducial cuts.

Results at partonic level
In this section we compare the Matrix results against those obtained after generating the Powheg hardest emission of eq. (2.35), often denoted as "Les Houches events" (LHE), and labelled in the figures as "MiNNLO PS (LHE)".  In figure 5 we show results for the invariant mass and the rapidity of the diphoton system, whereas in figure 6 we display the transverse momentum of the hardest and next-tohardest photon. We find good agreement between the Matrix and the MiNNLO PS (LHE) predictions, both for the central value and for the uncertainty bands due to scale variations. Notice that, at variance with similar comparisons for other processes not involving photons 15 We point out that such generation cut has no effect on the final kinematic distributions if the fiducial cut on the diphoton invariant mass at the analysis stage is greater than it.  Figure 6. Comparison between the distributions obtained from the fixed-order calculation from Matrix and the results after the generation of the hardest emission, for the transverse momentum of the hardest and next-to-hardest photon.
(like Drell-Yan or massive diboson production), this is not trivial due to the presence of an isolation criterion in the definition of the process.

Results after parton showering
In this section we show results obtained after the completion of the parton shower performed by Pythia8 [94,95]. We rely on the Pythia8 interface to the Powheg Box Res, provided by the main31 configuration file, distributed with Pythia8. The results presented in this section are obtained after switching off multi-parton interactions, QED radiation and hadronization effects, using the Monash tune [96]. We set the Pythia8 parameter POWHEG:pThard to 2 (i.e. we use the prescription introduced in section 4 of ref. [93]), and the SpaceShower:dipoleRecoil to 1.  transverse momenta of the two photons. We ascribe this behaviour to the fact that, with the increased multiplicity of the partonic/hadronic activity due to the Pythia8 shower, the photons in the showered events are less likely to satisfy the isolation criteria, giving a smaller cross section. We observed such pattern also for similar distributions computed with a completely independent code for diphoton production, built with the standard Powheg NLO+PS formalism, reassuring us that the above interpretation is correct, and that this behaviour is not due to some MiNNLO PS features. We leave further investigation on this aspect for future studies.

Comparison with ATLAS results
In this section we compare our results with those presented by the ATLAS Collaboration in ref. [74] obtained at 13 TeV. In order to make a fair comparison with the experimental data, we have added to the MiNNLO PS diphoton contribution, presented in this paper, the leading-order contribution coming from the gluon-induced production of a photon pair through a closed quark loop (denoted by gg → γγ in the following), which is of the same order of the NNLO corrections and is further enhanced by the particularly sizable gluon luminosity. As we consider this contribution only at leading order, we simulate it at LO+PS accuracy, setting the renormalization and factorization scales, and the upper limit for the transverse momentum of the subsequent shower evolution, equal to the invariant mass of the diphoton system. The analytic amplitudes for this process have been taken from ref. [18] and implemented in the POWHEG BOX RES framework, neglecting the top-quark loop contribution [17]. This contribution amounts for at most a few percent in certain regions of the kinematic distributions we are plotting. Its inclusion would improve the agreement with data. The analysis of our events is performed using Rivet [97], with the same Pythia8 settings as those in the previous section. The fiducial volume is defined by the following requirements p Tγ 1 > 40 GeV, p Tγ 2 > 30 GeV, ∆R γγ > 0.4 , |y γ | < 1.37 , 1.52 < |y γ | < 2.37. (4.1) We do not report here the exact photon-isolation criteria which are described in detail in section 4.1 of ref. [74].
In figures 9 and 10 we show the invariant mass and the transverse momentum of the  Figure 9.
Comparison between the ATLAS data and the distributions obtained with MiNNLO PS +Pythia8, for the diphoton invariant mass and its transverse momentum. Scalevariation bands are also shown, together with the statistical errors of the central value.  Comparison between the ATLAS data and the distributions obtained with MiNNLO PS +Pythia8, for the transverse momentum of the leading and subleading photon. Scalevariation bands are also shown, together with the statistical errors of the central value.
diphoton system, and the transverse momentum of the leading and subleading photon, respectively. Overall we find a reasonably good agreement between data and theoretical predictions.
For the invariant mass distribution m γγ we observe a good agreement in the bulk of the cross section. Given the cuts of eq. (4.1), the region where m γγ < 80 GeV is populated only by γγ + jet(s) events, therefore our results are only NLO+PS accurate, as confirmed also by the wider uncertainty bands. For m γγ 40 GeV, MiNNLO PS results overshoot ATLAS data by an amount compatible with what has been observed, for other predictions of similar accuracy, in ref. [74]. This is particularly true for the first bin. The fact that this region is characterized by a large NLO K-factor [13] hints at the possibility that the inclusion of higher-order corrections will improve the agreement with data. At large m γγ values we observe differences up to about 15%. Such differences might be due to the absence of higher-order effects. Top-quark mass effects above the threshold m γγ 2m t , that we are neglecting in the quark-induced 2-loop amplitudes as well as in the gg → γγ channel, can also induce differences at a few percent level.
As far as the transverse-momentum spectrum of the diphoton pair is concerned, our predictions are compatible with experimental data across most of the p Tγγ spectrum, within the estimated theoretical uncertainties. Due to the all-order treatment of soft and collinear radiation encoded in the MiNNLO PS Sudakov form factor and in the subsequent matching to Pythia, the region of small p Tγγ is in a fair agreement with the data, within the statistical error and the scale-variation band. For large values of p Tγγ , data and theoretical predictions agree well in shape but display a flat offset of about 10%.
The comparison of our results against ATLAS data for the p T spectrum of the leading and subleading photon is shown in figure 10. Data and theoretical predictions agree quite well over the full range of the transverse momentum.
We would also like to point out that the choices for the central value of the renormalization and factorization scales, in a fixed-order computation, can have a sizable impact on the theory predictions at large invariant mass and transverse momentum of the diphoton system, as discussed in detail e.g. in ref. [14]. The size of these effects are a consequence of a slowly convergent perturbative series. Assessing this type of effects by exploring different scale choices in our MiNNLO PS simulation goes beyond the scope of this paper, although this is a-priori possible, both for the handling of the large-p T limit of the singular component (first term of eq. (2.31)), as discussed for instance in ref. [73], as well as for the scales used in the non-singular contribution (last term in the same equation).
In figure 11 we show the differential cross section as a function of two angular observables: the azimuthal separation ∆φ γγ of the two photons and φ * η , defined as where ∆y γγ is the azimuthal angular separation of the two photons. Such variable, first introduced for Drell-Yan processes in ref. [98], is known to be sensitive to the same dynamics governing the p Tγγ spectrum, but it allows for a better resolution at small values of p Tγγ . The agreement between data and the theoretical predictions is rather good on the whole Comparison between the ATLAS data and the distributions obtained with MiNNLO PS +Pythia8, for the azimuthal separation of the two photons ∆φ γγ and for φ * η , as defined in eq. (4.2). Scale-variation bands are also shown, together with the statistical errors of the central value.
range within the statistical and scale-variation bands. We stress that, for ∆φ γγ → π, the NNLO fixed-order results for the above distributions are divergent, while we get a finite distributions due to the MiNNLO Sudakov form factor.

Conclusions
In this paper we have presented the implementation of a new MiNNLO PS event generator accurate at NNLO in QCD for the production of a photon pair in hadronic collisions, within the Powheg Box Res framework. This implementation is based on the MiNNLO PS formalism presented in refs. [50,51], that we have modified in order to deal with the new features needed to describe a process with photons at the Born level. We have devised and implemented a general method to deal with the presence of QED singularities at the Born level within the Powheg Box Res, and we have presented a new mapping from γγj to the γγ. The combined use of the new mapping and of the method to deal with the QED divergences allowed us to simulate diphoton production at NNLO+PS accuracy without introducing any generation or technical cut. Our generator produces weighted events that cover the entire phase space, thereby allowing the generated sample to be used with arbitrary fiducial cuts. Furthermore, we have introduced a few modifications to the MiNNLO PS formalism, that do not change the formal accuracy of the method and reduce the numerical impact of higher-order terms, thus leading to better agreement with fixed-order computations. We have also verified that the new formulation of MiNNLO PS , when applied to the original MiNNLO PS implementations with massive bosons, gives results fully compatible with the original one.
Finally, we have presented some distributions of physical interest obtained from the generated events both before and after passing them through the parton shower, and compared them with the ATLAS results at 13 TeV. Using the standard settings of Pythia8 when interfacing it with the Powheg Box Res, we observe a reduction of about 5-10% of the inclusive differential cross sections compared to the NNLO fixed-order ones. We ascribe this effect to the increased multiplicity of the partonic/hadronic activity coming from the shower causing the photons to less likely satisfy the isolation criteria. Further investigation of the matching procedure between Powheg Box Res+MiNNLO PS and Pythia8 is left for future studies.

Acknowledgments
We would like to thank S. Alioli, A. Broggio and J. Lindert for useful discussions. We also thank J.P. Guillet, P.F. Monni, and M. Wiesemann for carefully reading the manuscript and for their useful comments.
E.R. thanks P.F. Monni and M. Wiesemann for several discussions about theoretical aspects related to the MiNNLO PS method. In particular, the relevant aspects of the formulation of the MiNNLO PS method with the scale choice discussed in section 2.4, and their initial implementation in the code, were developed in collaboration with them.
The work of A.G. was supported by the ERC Starting Grant REINVENT-714788.
A The Φ γγj → Φ γγ phase space projection In section 2.3 it was shown that one needs a mapping from the Φ γγj phase space to the Φ γγ one, and this mapping should be smooth when p T → 0, i.e. when the final-state jet is soft or collinear with the incoming beams. This projection plays a more important role for the process under study than in previous processes dealt within the MiNNLO PS procedure. In fact, for diphoton production, we would like to avoid that the kinematics of a γγj event that is away from any collinear divergence is projected on a kinematics for γγ production, where the photons are quasi-collinear with the incoming beams, which would give divergent contributions for the pp → γγ matrix elements needed in the MiNNLO PS D terms (see section 2.4). This projection is determined by the requirements that it preserves the mass and rapidity of the color singlet, and the direction of one of the photons in the laboratory frame.
As far as the MiNNLO PS formulae are concerned, we just need to define a Φ γγj → Φ γγ mapping that projects a Φ γγj phase-space point to a Φ γγ one. In addition to this, we also provide the inverse mapping, that will be used to compute the Jacobian of the transformation in section A.3.

A.1 The Φ γγj kinematics
In the partonic center-of-mass frame, momentum conservation reads and s is the squared center-of-mass energy. In this section, at difference with the rest of the paper, we denote with p γ 1 and p γ 2 the momenta of the two photons, irrespectively of which one is the hardest or the softest. Using the standard FKS [86] radiation variables ξ, y and φ j , we write the jet momentum as where cos θ j = y. Using eq. (A.1), the momentum of the diphoton system is then given by that has invariant mass and rapidity given by In order to preserve the direction of one photon, we need to express the momenta of the involved particles in two frames where the diphoton system has the same rapidity. We choose to work in the center-of-mass frame of the γγj system because, in this frame, the momentum of the final-state parton has a simple representation in terms of the FKS variables, as given in eq. (A.3). We parametrize the momenta of the two photons in their center-of-mass frame as and we then perform a longitudinal boost with rapidity y CM γγ of eq. (A.5) to obtain We now impose that the momentum p γ 1 in the partonic center-of-mass frame of γγj has the same direction asp γ 1 . In this way, if the photon γ 1 of Φ γγ becomes too close to the incoming beams, then also the photon γ 1 of Φ γγj will be close to the beams, and will be suppressed by the Born suppression factors of eq. (2.36). To achieve this, we change the energy ofp γ 1 , introducing the dimensionless parameter E, such that the momentum of γ 1 in the center-of-mass frame kinematics of Φ γγj is given by We fix the value of E by imposing that the second photon stays massless. In fact, from the momentum conservation of eq. (A.1) and by imposing p 2 γ 2 = 0 we have Once the momenta p j , p γ 1 and p γ 2 , given by eqs. (A.3), (A.9) and (A.10) are known in their partonic center-of-mass frame, we can obtain their expression in the laboratory system by performing a longitudinal boost with the rapidity of the center-of-mass system.

A.3 The Jacobian of the mapping
In order to conclude the procedure to build Φ γγj from Φ γγ , we need to compute the Jacobian of the transformation we have outlined in the previous section. Following the Powheg notation [84,85], we write the infinitesimal phase space volume in terms of the Bjorken x as dΦ 3 ≡ dx ⊕ dx dΦ γγj = dx ⊕ dx dΦ rad dΦ γγ , is written using the FKS radiation variables. We would like to express dΦ γγ using the polar and azimuthal angle of one photon in the center-of-mass frame of the diphoton system, i.e. θ and ϕ , where it can be written as in terms of the angles of our mapping, i.e.θ andφ, that we have used to write eqs. (A.9) and (A.10). Using the fact that dΦ γγ is a Lorentz scalar, we need to compute the Jacobian J that relates the change of variables dΦ γγ = J 32π 2 d cosθ dφ .
(A. 15) In order to do this, we find the sequence of boosts from the partonic center-of-mass frame of γγj to the center-of-mass frame of the diphoton system. We first perform a longitudinal boost of the momenta p γ 1 and p γ 2 of eqs. (A.9) and (A.10) to the frame where the diphoton system has zero rapidity. The longitudinal rapidity of the boost is given by eq. (A.5) and we obtain p L γ 1 = E √ s 1, sinθ sinφ, sinθ cosφ, cosθ , Then we perform a transverse boost to the diphoton rest system with transverse velocity and we get In the diphoton rest system, p T γ 1 , written as a function of θ and ϕ , has components p T γ 1 ≡ sin θ sin ϕ , sin θ cos ϕ , cos θ .
and from these expressions, we can compute the Jacobian J This completes the calculation of the direct and inverse mapping Φ γγj ↔ Φ γγ and the construction of the three-body phase space, from one side, and of its inverse mapping, on the other.

B MiNNLO PS formulae
In this section, we provide the explicit expressions for the ingredients needed in the MiNNLO PS formulae, focussing in particular on their renormalization and factorization scale dependence. We use the same notation as ref. [50], and we refer the interested reader to this paper for further details. For ease of notation, we drop all the dependences on the kinematic variables, and we retain only that on the scales, which we assume to be equal, µ R = µ F = µ. In section 2.4, a physical quantity G computed at scale µ was indicated with [G(Φ γγ , p T )] µ . In the compact notation of this section, is is now simply written as G(µ). The explicit expression for the luminosity L, evaluated at scales p T , is where |M cd→γγ | 2 is the Born squared amplitude for the process cd → γγ (in this case c and d are a quark-antiquark pair),H cd→γγ encodes its virtual corrections up to two loops, 16 theC ij are the collinear coefficient functions and f [a] i is the PDF for the parton i in the hadron a. BothH cd→γγ andC ij can be written as an expansion in α S as In the above formula H (1) andH (2) are the one-and two-loop virtual corrections for the process qq → γγ written in the MS renormalization and subtraction schemes.
We stress again that the renormalization and factorization scales in eq. (B.3) are both set to p T , as indicated by the argument of the luminosity L. Finally, the convolution operator is defined, as usual, as (B.9) The expansion of D at a scale µ can be written as .

(B.12)
All the needed ingredients to compute eqs. (B.11) and (B.12) can be obtained by applying the DGLAP evolution equations to compute the scale dependence of the parton distribution function 14) and the renormalization group equation to compute the running of α S We then obtain where the parton distribution functions on the right-hand side are evaluated at scale µ F . In addition, from the corresponding expressions evaluated at scale p T given in refs. [50,51], we can compute dL(µ) dp T (1) and dS(µ) dp T (1) dS(µ) dp T = A (2) log Q 2 p 2 T +B (2) + β 0 dS(µ) dp T (1) (B.23)