Diphoton production at the LHC: a QCD study up to NNLO

We consider the production of prompt-photon pairs at the LHC and we report on a study of QCD radiative corrections up to the next-to-next-to-leading order (NNLO). We present a detailed comparison of next-to-leading order (NLO) results obtained within the standard and smooth cone isolation criteria, by studying the dependence on the isolation parameters. We highlight the role of different partonic subprocesses within the two isolation criteria, and we show that they produce large radiative corrections for both criteria. Smooth cone isolation is a consistent procedure to compute QCD radiative corrections at NLO and beyond. If photon isolation is sufficiently tight, we show that the NLO results for the two isolation procedures are consistent with each other within their perturbative uncertainties. We then extend our study to NNLO by using smooth cone isolation. We discuss the impact of the NNLO corrections and the corresponding perturbative uncertainties for both fiducial cross sections and distributions, and we comment on the comparison with some LHC data. Throughout our study we remark the main features that are produced by the kinematical selection cuts that are applied to the photons. In particular, we examine soft-gluon singularities that appear in the perturbative computations of the invariant mass distribution of the photon pair, the transverse-momentum spectra of the photons, and the fiducial cross section with asymmetric and symmetric photon transverse-momentum cuts, and we present their behaviour in analytic form.


Introduction
The production of photon pairs (diphotons) with high invariant mass at high-energy hadron colliders is a very relevant process in the context of both Standard Model (SM) studies and searches for new-physics signals.
Experimentally a pair of photons is a very clean final state, and photon energies and momenta can be measured with high precision in modern electromagnetic calorimeters. Since photons do not interact strongly with other final-state particles, prompt photons represent ideal probes to test the properties of the SM and corresponding theoretical predictions (see Refs. [1]- [8] for recent experimental analyses of Tevatron and LHC diphoton data). Measurements involving a pair of isolated photons have played a crucial role in the discovery at the LHC [9,10] of a Higgs boson, whose properties are compatible with those of the SM one. Studies of the Higgs boson properties in the diphoton decay mode have been performed [11,12]. Diphoton measurements (see, e.g., Refs. [13]- [18]) are also important in many new-physics scenarios, including searches for extra dimensions or supersymmetry. The relevance of LHC measurements of the diphoton invariant-mass spectrum is highlighted by the recent observation [19]- [23] of an excess of events with invariant mass of about 750 GeV that might have indicated the presence of resonances over the diphoton SM background. That observation raised a great deal of attention from December 2015 till the time of the 2016 Summer Conferences [24,25].
Owing to its physics relevance, the study of diphoton production requires accurate theoretical calculations which, in particular, include QCD radiative corrections at high perturbative orders. In high-energy collisions, final-state prompt photons with high transverse momentum can be originated through direct production from hard-scattering subprocesses and through fragmentation subprocesses of QCD partons. The theoretical computation of fragmentation subprocesses requires non-perturbative information, in the form of parton fragmentation functions of the photon, which typically has large associated uncertainties. However, the effect of fragmentation contributions is significantly reduced by the photon isolation criteria that are necessarily applied in hadron collider experiments to suppress the very large reducible background of 'non-prompt' photons (e.g., photons that are faked by jets or produced by hadron decays). Two such criteria are the so-called 'standard' cone isolation and the 'smooth' cone isolation proposed by Frixione [26]. The standard cone isolation is the criterion that is typically used by experimental analyses. This criterion can be experimentally implemented in a relatively straightforward manner, but it only suppresses part of the fragmentation contribution. By contrast, the smooth cone isolation (formally) eliminates the entire fragmentation contribution, although, due to the finite granularity of the detectors, it cannot be directly applied at the experimental level in its original form. Owing to the absence of fragmentation contributions, theoretical calculations are much simplified by using the smooth cone isolation, and it is relatively simpler to compute radiative corrections at high perturbative orders. Considering calculations at the next-to-leading order (NLO) in the QCD coupling α S , in Ref. [27] it was shown that, if the isolation is 'tight enough', the two isolation criteria lead to theoretical NLO results that are quantitatively very similar for various observables in diphoton production processes at high-energy hadron colliders.
In the present paper we deal with the diphoton production process pp → γγX, where p is a colliding proton and X denotes the inclusive final state that accompanies the γγ pair. QCD radiative corrections to this process were first computed up to the NLO in Ref. [28]. This is a complete NLO calculation of both the direct and fragmentation (specifically, single-and double-fragmentation) components of the cross section. The calculation is implemented at the fully-differential level in the numerical Monte Carlo code DIPHOX, which can be used to perform computations with any infrared and collinear safe isolation criteria (including the standard and smooth cone criteria). The DIPHOX calculation also includes the so-called box contribution [29] to the partonic channel gg → γγ (which is formally a contribution of next-to-next-to-leading order (NNLO) type), since this contribution can be quantitatively enhanced by the large gluon-gluon parton luminosity at high-energy hadron-hadron colliders. The next-order gluonic corrections to the box contribution (which are part of the N 3 LO QCD corrections to diphoton production) were computed in Ref. [30] (and implemented in the numerical program GAMMA2MC) and found to have a moderate quantitative effect. An independent NLO calculation [31] is implemented in the code MCFM, which, however, includes the fragmentation component only at the leading order (LO), while the box contribution is treated according to the GAMMA2MC code. A complete calculation at the NNLO of both the direct and fragmentation components still nowadays remains computationally very challenging.
Fragmentation contributions are absent by considering smooth cone isolation. In the context of smooth cone isolation, diphoton production at the LO (i.e. at O(α 0 S )) emerges via the quark-antiquark annihilation subprocess qq → γγ. The NLO QCD corrections are due to quarkantiquark annihilation and to new partonic channels (via the subprocess qg → γγq andqg → γγq) with an initial-state colliding gluon. At the NNLO the gg channel starts to contribute and it (the entire channel and not only its box contribution) can be fully consistently included in the perturbative QCD calculation. The large value of the gluon parton distribution function (PDF) of the colliding proton makes the gluon initiated channels important, especially at small and intermediate values of the diphoton invariant mass. The scattering amplitudes that are needed to evaluate the NNLO QCD corrections were computed in Ref. [29,32,33,34], and they were used in Ref. [35] to perform the first NNLO calculation of diphoton production in hadronic collisions. The NNLO calculation at the fully-differential level is based on the q T subtraction method [36] and it was implemented [35] in the numerical code 2γNNLO. More recently, an independent calculation of diphoton production at the NNLO has been performed in Ref. [37] by using the method of N-jettiness subtraction [38,39]. The NNLO calculation of Ref. [37] is implemented [40] in the MCFM program. An independent NNLO computation of diphoton production at the NNLO, based on the q T subtraction method, has been implemented in the general purpose NNLO generator MATRIX [41].
The transverse-momentum spectrum of the photon pair is sensitive to logarithmically-enhanced contributions at high perturbative orders. Transverse-momentum resummation at next-to-nextto-leading logarithmic (NNLL) accuracy for inclusive diphoton production was implemented [42] in the ResBos code, and more recently the complete NNLL calculation has been combined with the NNLO contributions and implemented in the numerical program 2γRes [43].
Lowest-order electroweak radiative corrections to diphoton production at LHC energies have been computed in Refs. [55,56]. They produce small effects on NNLO QCD results for inclusive diphoton production. The effects can increase by selecting photons with transverse momenta in the TeV region.
In this paper we present studies of NLO and NNLO QCD radiative corrections to inclusive diphoton production at LHC energies. The results at the NNLO are based on the theoretical calculation of Ref. [35], as implemented in the fully-differential Monte Carlo program 2γNNLO, and on the independent NNLO calculation implemented in MATRIX [41]. The first version of 2γNNLO, which was originally used in Ref. [35], had a numerical implementation bug that was subsequently corrected. A main result of Ref. [35] is that NNLO radiative corrections are substantial for diphoton kinematical configurations of interest at high-energy hadron colliders. The selected quantitative results that were presented in Ref. [35] are mostly related to diphoton kinematical configurations that are typically examined in Higgs boson searches and studies at the LHC. In this paper we consider more general kinematical configurations. In particular, we present NLO and NNLO studies related to photon isolation, and we discuss various detailed features of the theoretical results up to NNLO.
The results of the code 2γNNLO have been used by experimental collaborations in some of their data/theory comparisons. The analyses performed at the Tevatron ( √ s = 1.96 TeV) by the CDF [5] and D0 [6] collaborations and at the LHC ( √ s = 7 TeV and 8 TeV) by the ATLAS [4,8] and CMS [7] collaborations show that the inclusion of the NNLO corrections greatly improves the description of diphoton production data. The NNLO predictions have also been used in data analyses related to new-physics searches in high-mass diphoton events at √ s = 13 GeV [20].
In summary, these measurements prove that the NNLO results are important to understand phenomenological aspects of diphoton production.
The paper is organized as follows. Section 2 is devoted to a comprehensive study of photon isolation. In particular, we perform a detailed comparison of the standard and smooth isolation criteria in the context of perturbative QCD results at LO and NLO. We discuss the role of different partonic subprocesses within the two different isolation criteria. We also remark on the effects that are produced by the isolation parameters and by the kinematical selection cuts that are applied to the photons. Specifically, in Sect. 2.1 we introduce the photon isolation criteria and comment on some of their features. The QCD calculations that are used in our study are briefly described in Sect. 2.2. The quantitative results for total and differential cross sections are presented in Sect. 2.3. Section 2.3.1 is devoted to the LO and NLO results for total cross sections. Results for differential cross sections at LO and NLO are presented in Sects. 2.3.2 and 2.3.3, respectively. In Sect. 3 we present detailed NNLO results for diphoton production within the smooth isolation prescription. We study the perturbative stability of the results and we discuss related theoretical uncertainties by considering several observables that are relevant in diphoton production at hadronic colliders. We also discuss the comparison of the NNLO predictions to recent LHC data [4]. Specifically, results for total cross sections and differential cross sections are presented in Sects. 3.1 and 3.2, respectively. In Sect. 3.3 we present some results on the dependence on the isolation parameters. The comparison with LHC data is discussed in Sect. 3.4. In Sect. 3.5 we discuss the effects of (asymmetric and symmetric) photon transverse-momentum cuts on total and differential cross sections and, in particular, we comment on related perturbative (soft-gluon) instabilities. Finally, in Sect. 4 we summarize our results.

Photon isolation 2.1 Isolation criteria
Hadron collider experiments at the Tevatron and the LHC do not perform inclusive photon measurements. The background of secondary photons coming from the decays of π 0 , η, etc. overwhelms the signal by several orders of magnitude and the experimental selection of prompt diphotons requires isolation cuts (or criteria) to reject this background. The standard cone isolation and the smooth cone isolation are two of these criteria. Both criteria consider the amount of hadronic (partonic) transverse energy ‡ E had T (r) = i E had T i Θ(r − R iγ ) inside a cone of radius r around the direction of the photon momentum p γ . Then the isolated photons are selected by limiting the value of E had T (r).
The standard cone isolation criterion fixes the size R of the radius of the isolation cone and it requires where the isolation parameter E T max can be either a fixed value of transverse energy or a function of the photon transverse momentum p T γ (i.e., E T max = ǫ p T γ with a fixed parameter ǫ). A combination of these two options is also possible: for instance, E T max = 0.05 p T γ + 6 GeV is used in the study of Refs. [19,22].
Provided E T max is finite (not vanishing) standard cone isolation leads to infrared-safe cross sections [57] in QCD perturbation theory. Parton radiation exactly collinear with the direction of the photon momentum is allowed by the constraint in Eq. (1) and, as a consequence, the treatment of standard cone isolation within perturbative QCD requires the introduction of parton to photon fragmentation functions. Decreasing the value of E T max reduces and suppresses the effect of the fragmentation function (and of the corresponding partonic subprocesses).
The smooth cone isolation criterion [26] (see also Refs. [58,59]) also fixes the size R of the isolation cone and it requires E had T (r) ≤ E T max χ(r; R) , in all cones with r ≤ R , with a suitable choice of the r dependence of the isolation function χ(r; R). The two key properties [26] of the isolation function are: χ(r; R) has to smoothly vanish as the cone radius r vanishes (χ(r) → 0 , if r → 0 ), and it has to fulfil the condition 0 < χ(r; R) ≤ 1 (in particular, χ must not vanish) for any finite (non-vanishing) value of r. Since E had T (r) does not increase by decreasing r, in practice the requirement in Eq. (2) is effective only if χ(r; R) monotonically decreases as r decreases.
The smooth cone isolation criterion implies that, closer to the photon, less hadronic activity is allowed. The amount of energy deposited by parton radiation at angular distance r = 0 from the photon is required to be exactly equal to zero, and the fragmentation component (which has a purely collinear origin in perturbative QCD) of the cross section vanishes completely. The ‡ For each four-momentum p µ i , the corresponding transverse momentum (p T i ), transverse energy (E T i ), rapidity (y i ) and azimuthal angle (Φ i ) are defined in the centre-of-mass frame of the colliding hadrons. Angular distances R iγ are defined in rapidity-azimuthal angle space ( cancellation of perturbative QCD soft divergences still takes place as in ordinary infrared-safe cross sections, since parton radiation is not forbidden in any finite region of the phase space [26]. It is also preferable to choose isolation functions χ(r; R) with a sufficiently smooth dependence on r over the entire range 0 < r < R. In particular, large discontinuities of χ(r; R) at finite values of r are potential sources of instabilities [60] in fixed-order perturbative calculations. Small discontinuities of the function χ(r; R) (such as those in the discretized version [61] of smooth cone isolation) are instead acceptable.
A customary choice of the isolation function χ(r; R) is where the value of the power n is typically set to n = 1. We also consider the following isolation function: whose value depends on the ratio r/R (rather than r and R, independently). The two functions in Eqs. (3) and (4) are equal at the isolation cone boundary r → R (χ(r; R) → 1) and they behave similarly as r → 0 (χ(r; R) ∝ r 2n ).
Comparing the isolation requirements in Eqs. (1) and (2) by using the same values of R and E T max in both equations, we see that smooth cone isolation is more restrictive than standard cone isolation. Therefore, the following physical constraint applies: where dσ generically denotes total cross sections and differential cross sections with respect to photon kinematical variables, and the subscripts 'smooth' and 'standard' refer to smooth and standard isolation, respectively. Note that the isolation parameters R and E T max are set at the same values in the two isolated cross sections, dσ smooth and dσ standard , that are compared in the inequality (5) (e.g., the inequality is not necessarily valid if smooth isolation at a given value of E T max is compared with standard isolation at a different and smaller value of E T max ). An analogous reasoning applies to the cross section dependence on the isolation parameters E T max and R, since the isolation requirement can become more or less restrictive by varying these parameters. Therefore, we have the following physical behaviour: dσ is (R; E T max ) monotonically increases as R decreases (E T max fixed) , dσ smooth (R; E T max ; n) monotonically decreases as n increases (R and E T max fixed) , (8) and the subscript 'is' equally applies to both isolation criteria (e.g., 'is'='smooth' or 'is'='standard'). The relation (8) refers to the dependence on the power n in the case of the isolation function in Eqs.
The standard cone isolation criterion is simpler and, as stated in the Introduction, it is the criterion that is used in experimental analyses at hadron colliders (the actual experimental selection of isolated photons, including isolation requirements, is definitely much more involved than the relatively straightforward implementation of the criterion). The experimental implementation of smooth cone isolation (in its strict original form) is complicated § especially by the finite granularity of the Tevatron and LHC detectors.
Independently of the intrinsic differences between the standard and smooth cone isolation criteria, Eq. (5) implies that a reliable (theoretical or experimental) cross section result with smooth cone isolation represents a lower bound on the corresponding (i.e., with the same values of E T max and R) cross section result with standard cone isolation (this statement is valid within reliable estimates of related theoretical or experimental uncertainties). In particular, this also implies that a comparison between theoretical smooth isolation results and experimental standard isolation results leads to meaningful and valuable information.
The comparison between smooth and standard isolation can be sharpened by considering tight isolation requirements. As expected on general grounds, the differences between the two isolation criteria should be reduced in the case of tight isolation cuts. This expectation is confirmed by the diphoton studies in Refs. [27,67], which show theoretical NLO results that are similar for the two isolation criteria if the isolation parameters are tight enough (e.g., R = 0.4 and E T max < 5 GeV or E T max = ǫ p T γ with ǫ < 0.1 as in Sect. 11 of Ref. [27] or Sect. 4.3.1 of Ref. [67]).
In Sect. 2.3 we present detailed results for smooth and standard isolation. We study the dependence on the isolation parameter E T max and its effects on both total cross sections and differential cross sections in various kinematical regions.

Perturbative QCD calculations
In this paper we present quantitative results on QCD radiative corrections to diphoton production by using both smooth cone and standard cone isolations. We consider total cross sections (more precisely, fiducial cross sections) and differential cross sections as functions of various kinematical variables of the prompt photons. The kinematical variables that we use are the diphoton invariant mass M γγ , the photon polar angle θ * in the Collins-Soper rest frame [68] of the diphoton system, the azimuthal angle separation ∆Φ γγ between the two photons, the diphoton transverse momentum p T γγ and the transverse momenta p hard T γ and p sof t T γ (p hard T γ > p sof t T γ ) of the harder and softer photon of the diphoton pair. As previously specified, azimuthal angles and transverse momenta are defined in the centre-of-mass frame of the colliding hadrons.
The QCD results on smooth cone isolation are obtained by using the numerical programs 2γNNLO [35] and MATRIX [41], which include perturbative QCD contributions up to NNLO. The NNLO calculations are based on the q T subtraction method [36], which is applicable to hadroproduction processes of generic high-mass systems of colourless particles and it has been already applied to explicit NNLO calculations of several specific processes (Higgs boson [36,69] and vector boson [70,71] production, associated production of a Higgs boson and a vector boson [72,73], diboson production processes such as Zγ [51,52], W ± γ [52], ZZ [74,75], W + W − [76,77] and W ± Z [78] production, Higgs boson pair production [79], associated production of a Higgs boson § There is activity related to the experimental implementation [61,62,63,64] of the discretized version of smooth cone isolation. An experimental implementation of the smooth isolation criterion was done by the OPAL collaboration [65] in a study of isolated-photon production in photon-photon collisions at LEP. A discretized version of smooth isolation is used in the QCD calculation of Ref. [66]. pair and a W or a Z boson [80]). The q T subtraction method has recently been extended to heavy-quark production and applied to the NNLO calculation of top quark pair production [81].
In the case of diphoton production, the method combines the NLO cross section dσ γγ+jets for the production of the photon pair plus one or two jets (partons) with an appropriate subtraction counterterm (its explicit expression is given in Ref. [82]) and a hard-virtual contribution (see Ref. [83]) for diphoton production with no additional final-state jets (partons). In the code 2γNNLO the cross section dσ γγ+jets is computed up to NLO by using the dipole-subtraction method [84] as implemented in the diphoton calculation of Ref. [44]. In MATRIX the cross section dσ γγ+jets is computed at NLO by using the fully automated implementation of the dipole-subtraction method within the Monte Carlo program MUNICH ‡ , and all (spin-and colour-correlated) tree-level and one-loop amplitudes are obtained from OpenLoops [85]. The combination of dσ γγ+jets and its counterterm is numerically finite, although the two contributions are separately divergent in the limit of vanishing transverse momentum p T γγ of the photon pair. In 2γNNLO the cancellation of the separate divergences is numerically treated by introducing a lower limit q T cut on p T γγ (p T γγ > q T cut ) and considering small values of q T cut (formally performing the limit q T cut → 0). Decreasing the value of q T cut reduces the size of systematical errors (due to the finite value of q T cut ) at the expense of increasing computational time and resources to handle numerical instabilities. As a practical compromise, in our study we use a finite value of q T cut in the range q T cut = 0.1 GeV-0.2 GeV. The NNLO generator MATRIX implements instead a lower limit r cut on the ratio p T γγ /M γγ (p T γγ > r cut M γγ ) [41], and we use values in the range r cut = 0.08%-0.15%.
Owing to the finite values of q T cut or r cut , a systematical uncertainty of about ±O(1%) affects all the NNLO results § presented in this paper. The quoted systematical uncertainty is sufficient for the purpose of the studies that we present throughout the paper. It is substantially smaller than the size of additional (NLO) NNLO theoretical uncertainties that we find, for instance, by performing customary variations of the factorization and renormalization scales at (NLO) NNLO. We also note that, at fixed value of q T cut (r cut ), the NNLO systematical error on differential cross sections tends to be larger than the corresponding error on total cross sections. In particular, some specific observables (and, more precisely, their value in restricted kinematical regions) that are very sensitive to the detailed shape of the p T γγ distribution in the limit p T γγ → 0 can be affected by larger NNLO systematical errors due to the finite value of q T cut (r cut ): these observables may require refined numerical studies with smaller values of q T cut (r cut ). Nonetheless, these same specific observables are (expected to be) affected by increased theoretical uncertainties due to large higher-order perturbative corrections (instabilities). Examples of these observables are the differential cross section dσ/d∆Φ γγ at ∆Φ γγ ≃ π, and the differential cross sections dσ/dM γγ , dσ/dp hard T γ and dσ/dp sof t T γ (or related integrated quantities) in soft-gluon sensitive kinematical regions (see Sect. 3.5).
We use 2γNNLO and MATRIX to obtain results at LO, NLO and NNLO in QCD perturbation theory. At each order our calculation includes all the terms (and only those terms) that contribute to the total cross section at the corresponding perturbative order according to the formal expansion in powers of α S . Therefore (in the context of smooth cone isolation), only the initial-state qq ‡ MUNICH is the abbreviation of "MUlti-chaNnel Integrator at Swiss (CH) precision" -an automated parton-level NLO generator by S. Kallweit. § At NLO the generator MATRIX can also use the dipole-subtraction method [84], which does not require the regularization parameter r cut . Our diphoton NLO results have been quantitatively cross-checked by numerical comparisons of the calculations that use q T subtraction (2γNNLO, MATRIX) and dipole subtraction (MATRIX).
partonic channel contributes at LO, the initial-state qg andqg channels start to contribute at NLO, and the initial-state gg channel starts to contribute at NNLO. In particular the box contribution gg → γγ [29] only enters at NNLO, together with all the partonic subprocesses (e.g., including the gluon initiated subprocess gg → γγqq) that contribute at the same order. We do not include the NLO correction [30] to the box contribution, since it is part of the complete (and still unknown) N 3 LO corrections to inclusive diphoton production. All the partonic subprocesses are treated in the massless-quark framework with N f = 5 quark flavours. In particular, we do not include NNLO contributions with virtual top quarks since they are not yet fully known (e.g., the loop top quark contribution to the two-loop scattering amplitude qq → γγ has not yet been computed in complete form [86]). The effect of including the NLO correction [30] and the top quark correction to the (massless-quark) box contribution gg → γγ has been considered in the diphoton calculation of Ref. [37]. The top quark correction to the box contribution is also studied in Refs. [87] and [88].
The results on standard cone isolation are obtained by using the DIPHOX calculation [28], which includes QCD radiative corrections up to the NLO. The NLO calculation that is implemented in DIPHOX is based on a combined use [28,89] of the subtraction and phase-space slicing methods. The slicing approximation is applied [28] to the phase-space region where the minimum transversemomentum of the three final-state partons at NLO is smaller than p T m (formally considering the limit p T m → 0). Our numerical results are obtained by using p T m = 0.1 GeV, which is the default value of the slicing parameter (regulator) p T m that is suggested in the DIPHOX program. Since we are interested in an order-by-order comparison with smooth cone isolation results, the box contribution gg → γγ is not included in the DIPHOX NLO results. Note, however, that at LO and NLO all the initial-state partonic channels (qq, qg,qg, gg) contribute in DIPHOX because of the presence of a non-vanishing fragmentation component (though the fragmentation component is quantitatively suppressed by the isolation procedure). Owing to the double-fragmentation component, even the initial-state gg channel is not vanishing at LO (at NLO the initial-state gg channel contributes also through the single-fragmentation component).
We note that (due to the limited number of final-state partons in fixed-order computations) LO calculations do not cover the entire kinematical region for inclusive diphoton production. LO calculations give non-vanishing cross sections only in limited LO kinematical subregions. Outside these subregions, only the NLO results start to give non-vanishing cross section contributions. Therefore, outside the LO kinematical subregions the NLO (NNLO) results effectively represent an LO (NLO) perturbative QCD prediction. In spite of the effective meaning of the results, we always use the default labels LO, NLO and NNLO according to the perturbative order in which the results contribute to the inclusive (total) cross section. For example, in the case of the ∆Φ γγ distribution, the LO kinematical subregion has ∆Φ γγ = π: the region of small values of ∆Φ γγ receives contributions only from the NLO and NNLO results that represent effective LO and NLO predictions, respectively. We also note that the LO kinematical subregions can be different in the context of smooth cone and standard cone isolation. For instance, at LO p T γγ = 0 and p hard T γ = p sof t T γ in the case of smooth cone isolation, whereas p T γγ can be different from zero and p hard T γ ≥ p sof t T γ ≥ p hard T γ − E T max in the case of standard cone isolation (non-vanishing values of p T γγ are produced through the LO fragmentation component of the cross section). We also comment about the dependence on the isolation parameters R and E T max (the comment has similarities with our previous comment on LO kinematical subregions). The LO results are independent of the size R of the isolation cone. The LO cross section depends on E T max in the case of standard cone isolation (the E T max dependence is due to the fragmentation component), while it is independent of the value of E T max in the case of smooth cone isolation.
In Eqs. (5)- (8) we have listed some constraints on isolated photon cross sections. These are physical constraints (properties) in the sense that they apply to physical ('positive definite') events: if the isolation requirements are more (less) restrictive, the selected number of isolated photons decreases (increases). Such properties do not 'a priori' apply to theoretical calculations based on perturbative QCD (of course, if the constraints are not fulfilled, the perturbative QCD calculation is not a reliable approximation of the physical result) since, beyond the effective LO approximation, they involve negative-weight contributions (which are due to virtual radiative corrections and to subtraction contributions related to the factorization procedure of parton distribution functions and fragmentation functions).
As is well known, fixed-order perturbative results can show unphysical behaviour in specific kinematical regions. In the case of diphoton production, for instance, it is known [28] that the differential cross sections at small values of p T γγ (p T γγ → 0) or large values of ∆Φ γγ (∆Φ γγ → π) cannot be reliably computed in fixed-order perturbation theory: the disease is due to large logarithmic corrections that have to be resummed to all perturbative orders [42,43] to recover the physical behaviour and predictivity.
In the presence of photon isolation cuts, perturbative computations of total cross sections can also show some misbehaviour. Isolated (with both smooth and standard isolation) photon cross sections fulfil the physical constraint σ is (R; E T max ) < σ inc , where σ inc is the inclusive (nonisolated) cross section. Nonetheless, in the case of standard cone isolation at NLO this constraint is violated [57] at very small values (R ∼ < 0.1) of the radius R of the isolation cone. The violation is due to large logarithmic (ln R) corrections, and the physical behaviour is recovered through a corresponding all-order resummation procedure [90]. We expect a (qualitatively) similar ln R dependence in the case of smooth cone isolation, though in the present paper we do not consider very small values of R.
Violation of expected properties is not necessarily related to logarithmically-enhanced QCD corrections: it can be simply a consequence of a slowly convergent (toward a reliable estimate of the physical result) perturbative expansion. The LO, NLO and NNLO results obtained in Ref. [35] with smooth cone isolation certainly indicate that we are not dealing with a fastly convergent perturbative expansion in the case of diphoton production at high-energy hadron colliders. Additional complications can occur in direct comparisons (as in Eq. (5)) between calculations with smooth and standard isolation. The complications follow from the fact that such a comparison does not involve ingredients that are in one-to-one correspondence. Owing to the presence of the fragmentation component, at each perturbative order the standard isolation result depends on the photon fragmentation function, on the corresponding factorization scale µ f rag and on related partonic subprocesses. The poorly known fragmentation function certainly affects the standard isolation results and, especially, the comparison with smooth isolation results (which have no analogue of the fragmentation function).
Throughout the paper we comment on these and additional points related to physical behaviour and perturbative QCD calculations

Quantitative results
In our theoretical study of standard and smooth isolation we consider isolated diphoton production in pp collisions at the centre-of-mass energy √ s = 7 TeV. We apply the following kinematical cuts on photon transverse momenta and rapidities: p hard T γ ≥ 25 GeV, p soft T γ ≥ 22 GeV and the rapidity of both photons is limited in the range |y γ | < 2.37. The minimum angular distance between the two photons is R min γγ = 0.4.
These are basically the kinematical cuts used in the ATLAS Collaboration study of Ref. [4]. The analysis of Ref. [4] is restricted to a smaller rapidity region since it excludes the rapidity interval 1.37 < |y γ | < 1.52, which is outside the acceptance of the electromagnetic calorimeter. For the sake of simplicity, in this section we do not consider such additional rapidity restriction; the rapidity restriction is instead applied in the results of the following Sect. 3.
In the perturbative calculation, the QED coupling constant α is fixed at the value α = 1/137. We use the MMHT 2014 sets [91] of parton distribution functions (PDFs), with parton densities and α S evaluated at each corresponding perturbative order (i.e., we use the (k + 1)-loop running α S at N k LO, with k = 0, 1, 2). We use the NLO photon fragmentation functions of Ref. [92], and specifically the BFG set II (we have checked that BFG set I leads to very small quantitative differences). The central value µ 0 of the renormalization scale (µ R ), PDF factorization scale (µ F ) and fragmentation function scale (µ f rag ) is set to be equal to the invariant mass of the diphoton system, µ 0 = M γγ . We compute scale variation uncertainties by considering the two asymmetric scale configurations with More precisely, we have considered independent scale variations by a factor of two up and down with respect to the central value µ 0 . We find a common overall behaviour of the cross sections: they increase by decreasing µ R and decrease by decreasing µ F or µ f rag . Therefore the two asymmetric scale configurations are those that maximize the scale dependence within the considered scale variation range. The sole exception regards the invariant mass cross section at large values of M γγ (M γγ ∼ > 200 GeV): this kinematical region is sensitive to larger values of the parton momentum fraction x in the PDFs and, as a consequence, it turns out that the cross section decreases by increasing µ F .
The radius R of the isolation cone is fixed at R = 0.4. We study the isolation parameter dependence by considering values of E T max in the range between 2 GeV and 10 GeV (in this section we mainly show results at the extreme values in this range). In the case of smooth cone isolation, we use the isolation function χ(r; R) in Eq. (4) and we examine the cross section dependence on the power n that controls the shape of χ(r; R). In Fig. 1 we show the r dependence of the isolation function for some selected values of the power n and the fixed value R = 0.4 of the radius of the isolation cone. We note that the value R = 0.4 is sufficiently small so that the two isolation functions in Eqs. (3) and (4) quantitatively coincide at the percent level in the case with n = 1.
We present perturbative QCD results for both standard and smooth isolation and, in particular, we perform comparisons by using the same value of the isolation parameter E T max for both criteria. The comparison between the two criteria can also be performed differently. For instance, having fixed the values of R and E T max for standard isolation, one can use different values of isolation parameters (R, E T max and n) for smooth isolation to the purpose of trying to obtain similar quantitative results for the two criteria, as a pragmatic approach to mimic the standard cone isolation that is used in experimental conditions. We think that our comparison with the same value of E T max (and R) is more informative to investigate and understand differences and similarities between perturbative QCD results for the two criteria.
The QCD results on standard cone isolation depend on the parton-to-photon fragmentation function D a/γ (z; µ f rag ) (a = q,q, g), z being the photon momentum fraction with respect to the momentum of the fragmenting parton a. Owing to the isolation procedure, the value of z is bounded by a minimum value z min (1 ≥ z ≥ z min ), and this leads to a quantitative suppression of the fragmentation component of the diphoton cross section. The typical value of z min is z min ∼ p T γ /(p T γ + E T max ), p T γ being the transverse momentum of the photon that is involved in the fragmentation process. In our quantitative study we use relatively-large values of p T γ (i.e., typically, p T γ > 22 GeV) and relatively-small values of E T max . Therefore, z min is always large (z min ∼ > 0.9 at E T max = 2 GeV, and still z min ∼ > 0.7 at E T max = 10 GeV), and the suppression factor ‡ due to α S α D a/γ is sizeable (roughly one order of magnitude or more, depending on E T max ) [92]. We note that at such high values of z the quark (or antiquark) fragmentation function D q/γ (or Dq /γ ) is much larger (roughly by more than a factor of ten) than the gluon fragmentation function D g/γ [92]. In our calculation we consistently (according to the formal perturbative expansion) include all the fragmentation functions. However, due to the dominance of D q/γ and Dq /γ , in all our qualitative (or semi-quantitative) comments we neglect the effect of D g/γ (i.e., we can assume that only D q/γ and Dq /γ contribute). We also note that, because of QCD scaling violation, at high values of z, D a/γ (z; µ f rag ) increases (although weakly) by increasing µ f rag . ‡ At the formal level αS α D a/γ is the order of magnitude of the ratio between the fragmentation component and the direct component.

Total cross sections at LO and NLO
We begin the presentation of our quantitative results by considering the total cross section (namely, the fiducial cross section for the applied kinematical cuts on the photons). Values of total cross sections at LO and NLO for both smooth and standard isolation are reported in Table 1.
Using smooth cone isolation, the total cross section at LO is σ LO smooth = 10.56 pb +10.7 % −12.0 % (scale), where the percentage variation refers to the scale dependence of the result § . We note that σ LO smooth is independent of E T max and of the isolation function χ(r; R). The cross section is produced only by the initial-state qq channel through the partonic subprocess qq → γγ, and the scale dependence of σ LO smooth (which is entirely due to variations of µ F in the quark and antiquark PDFs) is quite small. All these features are a very crude approximation of the diphoton production dynamics.
Using standard cone isolation, the LO cross section σ LO standard depends on E T max and it also depends on all the three scales µ F , µ R and µ f rag (µ R and µ f rag enter through the fragmentation component). The scale dependence of σ LO standard is relatively similar to that of σ LO smooth : we find +14.5 % −14.3 % (scale) and +25.0 % −20.8 % (scale) with E T max = 2 GeV and E T max = 10 GeV, respectively. The E T max dependence of σ LO standard is instead more 'surprising'. Considering E T max = 2 GeV (very tight isolation), σ LO standard is slightly larger, by about 15%, than σ LO smooth (actually the two cross sections are very similar within the corresponding scale dependence). However, σ LO standard increases by a factor of about 1.6 in going from E T max = 2 GeV to E T max = 10 GeV and, at E T max = 10 GeV, σ LO standard is roughly a factor of 2 larger than σ LO smooth .
Since σ LO smooth does not depend on the isolation parameters (R, n, E T max ), there is obviously no way to approximate (or, mimic) the quantitative value of σ LO standard at E T max = 10 GeV by using smooth cone isolation with different isolation parameters (e.g., a value of E T max larger than 10 GeV).
The value E T max = 10 GeV is not particularly large and it cannot be regarded as a very loose isolation parameter. Therefore, on physical grounds, we do not expect two actual features of the LO result: the large difference between σ LO standard and σ LO smooth at E T max = 10 GeV, and the large E T max dependence of σ LO standard in going from E T max = 2 GeV to E T max = 10 GeV. We mean that physical events are not expected to have these features, since such features cannot be regarded as a physical consequence of the hadronic activity inside the photon isolation cones and of its detailed structure (which leads to an ensuing sensitivity to the isolation criteria). These observed LO features require some explanation, which we are going to discuss.
At the LO, in the context of standard isolation, the distinction between direct and fragmentation components is unambiguous, and the direct component exactly coincides with the entire contribution to the smooth isolation result. The double-fragmentation component always gives a small contribution (few percent at E T max = 10 GeV and few permill at E T max = 2 GeV) to σ LO standard . Therefore, the E T max dependence of σ LO standard is due to the single-fragmentation component, whose contribution to σ LO standard is small (of the order of 10%) at E T max = 2 GeV (because of the suppression due to the small value of E T max ) and sizeable at E T max = 10 GeV. At the larger value of E T max the direct and single-fragmentation components have the same size, but this is not due to the fact that the fragmentation function is particularly large. At the LO the direct component only involves the partonic process whereas the (single) fragmentation component also involves the partonic process where the notation q → γ + X denotes the fragmentation of the final-state quark q into a photon through D q/γ (a similar notation is used in subsequent equations for partonic subprocesses). The PDF luminosity (the detailed definition of PDF luminosities is presented in Eq. (17)) L qq of the initial-state direct subprocess is sizeably smaller than the luminosity L qg of the initial-state fragmentation subprocess (this follows from the smaller size of the antiquark PDF with respect to the gluon PDF at small values, such as x ∼ 10 −2 , of parton momentum fraction x in the PDFs), and the suppression due to the isolated fragmentation function D q/γ is compensated by the increased size of L qg with respect to L qq (the suppression from D q/γ is much stronger by decreasing E T max , because z > z min and z min increases by decreasing E T max ).
According to our discussion, the presence of the fragmentation process of Eq. (10) in the LO standard cone isolation explains the quantitative dependence of σ LO standard on E T max and the quantitative differences between σ LO standard and σ LO smooth . At the same time, our discussion is useful to anticipate expected features of the NLO results. The NLO calculation within smooth cone isolation includes the partonic process qg → γγq .
In the kinematical configuration where the final-state quark is inside the isolation cone of one photon, the process in Eq. (11) roughly corresponds to the LO fragmentation process of Eq. (10): therefore we expect that the NLO smooth isolation result receives a large NLO correction from this kinematical configuration of this process, in such a way to reduce the observed LO 'deficit' with respect to standard cone isolation (in other words, qg → γγq is suppressed with respect to qq → γγ by an extra power of α S , but this suppression is compensated by the increased luminosity of the partonic initial state). Moreover, the kinematical configuration where the final-state quark is outside the isolation cones of both photons contributes through the process of Eq. (11) to the NLO calculation of both smooth and standard isolations. This is not a very limited kinematical region (the size R of the isolation cones is not large) and, due to the large value of the L qg luminosity, it gives a sizeable and E T max independent NLO contribution to both isolation prescriptions. It follows that, for both isolation criteria, we expect large NLO corrections and a much reduced E T max dependence with respect to the LO result. As we are going to show shortly, the expectations of our discussion are confirmed by the actual NLO results.
We note that the LO results that we have presented are obtained by using LO PDFs and the NLO BFG fragmentation functions, since LO fragmentation functions are not readily available in the default setup of the DIPHOX code. This mismatch (at the strictly formal level) of perturbative order in the fragmentation functions should not strongly affect the main features of the LO standard isolation results. More generally, standard cone isolation results are certainly affected by an additional uncertainty, which is due to the poorly known fragmentation functions, that is difficult to be estimated at the quantitative level. The recent Ref. [93] presents a very brief overview on prospects for improving the determination of the photon fragmentation functions.
The value of the NLO total cross section σ N LO , including its corresponding scale variation dependence, is reported in Table 1 Table 1 and Fig. 2): the differences are at most at the level of 2-3% and they are sizeably smaller than the scale dependence of σ N LO . The E T max dependence is small: the NLO cross section for smooth (standard) isolation increases by a factor of 1.06 (1.07) in going from E T max = 2 GeV to E T max = 10 GeV. These features are in qualitative agreement with physical expectations.
The NLO radiative corrections are large (as expected from our previous discussion). Considering the ratio K N LO = σ N LO /σ LO , the value of K N LO smooth is roughly 3 and the values of K N LO standard are approximately 2.6 (E T max = 2 GeV) and 1.7 (E T max = 10 GeV). In the case of the smooth isolation criterion, a sizeable part (roughly 50%) of the NLO total cross section is due to the qg initial-state partonic channel (which is absent at the LO), and also the LO qq channel receives sizeable (roughly 50%) NLO corrections that increase the total cross section. In the case of standard cone isolation, at the NLO the distinction between direct and fragmentation components is no longer unambiguous ('physical'): it depends on the factorization scheme and on the fragmentation scale µ f rag . Within the MS factorization scheme and the scale variation range that we use, the direct component contributes about 45-65% (85-90%) of σ N LO standard with E T max = 10 GeV (E T max = 2 GeV). Since σ N LO standard and σ N LO smooth are very similar, the similarity is the consequence of a non-trivial interplay between the direct and fragmentation contributions to σ N LO (especially in the case with E T max = 10 GeV) and, in particular, the LO equivalence between the smooth cone result and the direct component of the standard cone result is definitely lost at the NLO.
We note that the scale dependence of the total cross section has a similar size at LO and at NLO, and it is much smaller than the size of the NLO corrections. This implies that the scale dependence of σ N LO cannot be consistently regarded as a reliable estimate of uncalculated higherorder radiative corrections to σ: the 'true' theoretical uncertainty of σ N LO is certainly larger than the NLO scale dependence that we have computed. A similar comment applies to the scale dependence of the results for the differential cross sections that we present in the following.
We now discuss the dependence of the smooth isolation cross section σ N LO smooth on the power n of the isolation function χ(r; R) = (r/R) 2n (see Eq. (4)). The NLO results of the total cross section for some selected values of n in the range 0.1 ≤ n ≤ 4 are reported in Fig. 2. We note that the n dependence of σ N LO smooth is small (in particular, it is smaller than the scale dependence) within this range of values of n. Specifically, considering the interval 1/2 < n < 2, the central value of σ N LO smooth varies by about 3% (10%) if E T max = 2 GeV (E T max = 10 GeV). Most of the qualitative features of the results in Fig.2 are consistent with physical expectations. The photons are more isolated by increasing n at fixed E T max and, consequently, σ N LO smooth monotonically decreases (in agreement with the physical behaviour in Eq. (8)). Moreover, by decreasing E T max the photons are more isolated and consequently the total cross section is less sensitive to variations of the power n. Nonetheless, we note that, by decreasing the value of n, σ N LO smooth tends to become larger than σ N LO standard , thus violating the physical constraint in Eq. (5). This feature deserves some comments, which are presented below.
The perturbative dependence of σ smooth at very small (n ≪ 1) or very large (n ≫ 1) values of n can be understood in relatively-simple terms. If n is very small, the isolation function χ(r; R) = (r/R) 2n is approximately equal to unity with the exception of the angular region of very small values of r: therefore, the n dependence of σ smooth is very sensitive to radiation of partons that are collinear to the photon direction. If n is very large, the isolation function very strongly suppresses parton radiation inside the photon isolation cone: therefore, the n dependence of σ smooth is very sensitive to soft parton radiation. The dominant effects of soft and collinear radiation can be easily computed at NLO (see Refs. [26,57,94]). We consider the NLO correction δ N LO = (σ N LO − σ LO )/σ LO , and we limit ourselves to sketch the dominant n dependence of the soft (δ N LO,soft smooth ) and collinear (δ N LO,coll smooth ) contribution to δ N LO within smooth isolation. We have where Q is the typical hard scale of the cross section (the scale is of the order of the minimum value of p hard T γ ) and we have considered small values of E T max (by neglecting relative corrections of O(E T max /Q)).
The proportionality factor that is not explicitly denoted on the right-hand side of Eq. (12) depends on the LO cross section for the partonic process qq → γγ. The soft contribution in Eq. (12) is negative. It is due to a strong kinematical mismatch between (negative) soft virtual radiation (one-loop corrections in the subprocess qq → γγ), which is not affected by isolation, and (positive) soft real radiation (the subprocess qq → gγγ at the tree level), which is strongly suppressed by isolation. We note that δ N LO,soft smooth is proportional to n, so that eventually σ N LO smooth diverges to '− ∞' in the limit n → +∞. This NLO divergence is the perturbative signal of the infrared unsafety of the isolated cross section in the limit of completely isolated photons (no accompanying transverse energy inside the isolation cone). We observe that δ N LO,soft smooth is proportional to R 2 , so that the soft contribution is strongly suppressed if the photon isolation cone has a small radius R. We also note that δ N LO,soft smooth is due to subprocesses with a qq initial state. The subprocess qg → qγγ is formally subdominant in the soft limit (n ≫ 1), but it represents a sizeable quantitative correction to δ N LO,soft smooth because of the increased PDF luminosity of the qg initial state. The R 2 -suppressed dependence of δ N LO,soft smooth and the large size of the (positive) correction to it from the qg initiated subprocess explain why the results for σ N LO smooth in Fig. 2 have a small dependence on n at relativelylarge values of n (e.g., n ≃ 4).
The collinear contribution in Eq. (13) is relevant to discuss the n dependence of the results in Fig. 2 at small values of n. We note that the NLO contributions from the initial-state qq channel (e.g., qq → gγγ) are subdominant in the limit n ≪ 1. The contribution in Eq. (13) is due to real radiation of a collinear quark or antiquark inside the photon isolation cone through the partonic processes qg → qγγ andqg →qγγ (the proportionality factor that is not explicitly denoted in the right-hand side of Eq. (13) depends on the LO cross section for the partonic processes qg → qγ and qg →qγ). Therefore, δ N LO,coll smooth is positive and independent of R. Moreover, δ N LO,coll smooth is proportional to E T max , so that its induced n dependence is reduced by decreasing E T max (in agreement with the results in Fig. 2 at small n). Owing to its dependence on n, δ N LO,coll smooth sizeably increases by decreasing n at fixed E T max and, eventually, δ N LO,coll smooth (and, consequently, σ N LO smooth ) diverges to '+∞' in the limit n → 0. Since σ N LO smooth becomes arbitrarily large by decreasing n, it is obvious that at sufficiently small values of n the physical requirement σ smooth < σ standard (see Eq. (5)) is unavoidably violated. This misbehaviour of σ N LO smooth at small values of n implies that beyond-NLO contributions are relevant. Indeed, each higher-order contribution is equally misbehaved at n ≪ 1: the N k LO correction is proportional to (α S /n) k (because of multiple collinear radiation inside the photon isolation cone) an it cannot be regarded as a small correction if n ∼ < α S (i.e., n ∼ < 0.1). In principle, the perturbative treatment at small values of n can be improved by a proper all-order resummation of these collinear contributions of O((α S /n) k ). However, such resummation treatment cannot be pursued for arbitrarily small values of n since it unavoidably fails in the limit n → 0, because of non-perturbative photon fragmentation effects (smooth isolation with n = 0 requires photon fragmentation functions since it is equivalent to standard isolation).
We can draw some conclusions from our discussion on small values of n. Owing to the physical requirements in Eqs. (5) and (8), in principle cross sections for standard and smooth isolation tend to agree at very small values of n. However, fixed-order QCD computations for smooth isolation are not reliable if n ≪ 1 (they are affected by large higher-order corrections) and, in particular, they can violate the physical constraint in Eq. (5). In practice, to the purpose of approximating the standard isolation criterion within fixed-order QCD calculations it is more appropriate to consider smooth isolation with values of n that are not too small. From the results in Fig. 2, we can conclude that the total cross sections σ N LO smooth and σ N LO standard quantitatively agree if n ≃ 1.
In the following we consider differential cross sections with respect to various kinematical variables and we limit ourselves to present smooth isolation results with n = 1. We have checked that the shape of the various NLO kinematical distributions is very little affected by variations of n within the range 0.1 ∼ < n ∼ < 4. At the NLO, variations of n basically produce overall normalization effects, whose size corresponds to the n dependence of σ N LO smooth that is observed in Fig. 2. Figure 3: The M γγ differential cross section for E T max = 2 GeV (left panel) and E T max = 10 GeV (right panel) and with the same photon kinematical cuts as in Fig. 2. The scale variation bands of the LO and NLO results for smooth and standard isolation are as follows: LO smooth isolation (black dashed), LO standard isolation (light-blue dashed), NLO smooth isolation (red solid) and NLO standard isolation (blue solid).

Differential cross sections at the LO
In Fig. 3 we present the LO and NLO results (including their scale variation dependence) for the differential cross section with respect to the diphoton invariant mass M γγ . We consider two different values, 2 GeV ( Fig. 3-left) and 10 GeV (Fig. 3-right), of E T max and we use both standard and smooth isolation. In Fig. 4 we present the analogous results for the differential cross section with respect to the angular variable cos θ * in the Collins-Soper rest frame. The results in Figs. 3 and 4 are obtained by numerical integration over small bins in M γγ and cos θ * , respectively: we use bins of constant size equal to 2 GeV for M γγ and 0.08 for cos θ * . In the following we discuss LO and NLO differential cross sections in turn. At the LO, we preliminarily note that standard and smooth isolation results for the differential cross sections have qualitatively similar shapes, with differences of the overall normalization that are quite similar to the quantitative differences (which we have previously discussed) of the corresponding LO total cross sections. Figure 4: The cos θ * differential cross section for E T max = 2 GeV (left panel) and E T max = 10 GeV (right panel) and with the same photon kinematical cuts as in Fig. 2. The scale variation bands of the LO and NLO results for smooth and standard isolation are as follows: LO smooth isolation (black dashed), LO standard isolation (light-blue dashed), NLO smooth isolation (red solid) and NLO standard isolation (blue solid).
We first consider smooth cone isolation. Owing to transverse-momentum conservation in the corresponding LO partonic process qq → γγ, the diphoton azimuthal separation is ∆Φ γγ = π and the transverse momentum of the diphoton system is p T γγ = 0. The corresponding differential cross sections dσ LO smooth /dx, with x = ∆Φ γγ or x = p T γγ , are simply proportional to a δ-function (δ(∆Φ γγ − π) or δ(p T γγ )) and the proportionality factor is the LO total cross section σ LO smooth . The double differential cross section with respect to M γγ and the scattering angle θ S is where √ s is the centre-of-mass energy of the hadronic collision (e.g., √ s = 7 TeV as in Fig. 3), τ = M 2 γγ /s and the photon scattering angle θ S is defined in the centre-of-mass frame of the LO partonic collision, and it is related to the diphoton rapidity separation ∆y γγ : We remark that, in the case of the LO smooth isolation cross section in Eq. (14), θ S and the Collins-Soper polar angle θ * actually coincides (|cos θ S | = |cos θ * |), since p T γγ = 0. In the righthand side of Eq. (14) we have not denoted an overall proportionality factor of order unity, which is independent of the kinematical variables s, M γγ and θ S . The angular dependent function and it is specific for the angular distribution that is dynamically produced by the partonic process qq → γγ. The function L qq (τ ; µ F ) is the qq PDF luminosity and e q is the quark electric charge in units of the positron charge (e q = 2/3 for up-type quarks). The PDF luminosity for the collision of two partons a and b is defined as where f a/h (x, µ F ) is the PDF of the parton a in the hadron h.
A distinctive feature of the right-hand side of Eq. (14) and, hence, of the double differential cross section is its exactly factorized dependence on M γγ and cos θ S . The M γγ dependence is controlled by the luminosity L qq : by decreasing M γγ , L qq rapidly increases. The cos θ S (or cos θ * ) dependence is controlled by g qq (θ S ): by increasing |cos θ S |, g qq sharply increases and it becomes singular at |cos θ S | = 1 (because of the 'unphysical' behaviour of 2 → 2 'massless' parton scattering). The computation of the single differential cross sections dσ/dM γγ and dσ/d cos θ * requires the application of kinematical cuts to select the hard-scattering regime (i.e., values of M γγ in the perturbative region and values of cos θ * sufficiently far from the forward/backward scattering singularity). The simplest type of kinematical cuts is a minimum value of M γγ and a maximum value of |cos θ S |. These kinematical cuts, which preserve the factorized structure of Eq. (14) with respect to the M γγ and cos θ S dependence, lead to differential cross sections dσ/dM γγ and dσ/d cos θ * that are simply proportional to L qq (τ ) and g qq (θ * ), respectively. The shapes of these differential cross sections are different (especially in the case of the cos θ * distribution) from those observed in Figs. 3 and 4. The differences originate from the kinematical cuts on the photon transverse momenta and rapidities that are described at the beginning of Sect. 2.3, and that are actually used in the computation of the differential cross sections of Figs. 3 and 4.
We first discuss the effect of the transverse-momentum (p T ) cuts p hard T γ ≥ p H (specifically p H = 25 GeV) and p sof t T γ ≥ p S (specifically p S = 22 GeV). Since p T γγ = 0 (i.e., p hard T γ = p sof t T γ ), only the value of p H is effective, and we have the LO constraint where M LO dir = 2p H (specifically M LO dir = 50 GeV) and we have simply used 2p hard T γ = M γγ sin θ S . Since sin θ S < 1 and M 2 γγ < s, the constraint in Eq. (18) implies a lower boundary M γγ ≥ M LO dir on M γγ (the LO smooth isolation cross section in Fig. 3 vanishes for M γγ < 50 GeV) and an upper boundary |cos θ S | < 1 − (M LO dir ) 2 /s on |cos θ S | (this corresponds to 1 − |cos θ S | ∼ < 3 · 10 −5 in Fig. 4). More importantly, the constraint in Eq. (18) correlates the M γγ and cos θ S dependencies, thus destroying the factorized structure in the right-hand side of Eq. (14) and leading to relevant effects on the shape of the single differential cross sections.
This LO behaviour of dσ/dM γγ is visible in the smooth isolation results of Fig. 3. Actually, the vanishing behaviour of dσ LO /dM γγ in the limit M γγ → M LO dir is so steep that it is not clearly recognizable in the invariant-mass bin 50 GeV< M γγ < 52 GeV closest to M LO dir =50 GeV. This vanishing behaviour is more evident by using M γγ -bins with a much smaller bin size (see Fig. 10 at the end of Sect. 2.3.3).
The effect of the p T cut constraint in Eq. (18) on dσ/d cos θ * is even much more relevant than the effect on dσ/dM γγ . In the computation of dσ LO smooth /d cos θ S from Eq. (14), the PDF luminosity L qq (τ = M 2 γγ /s) is integrated over M γγ up to a lower limit, M γγ > M LO dir / sin θ S , that depends on θ S : therefore large values of |cos θ S | are suppressed, while small values of |cos θ S | are relatively enhanced by the increasing size of the PDF luminosity for decreasing values of M γγ . This PDF modulation of the cos θ S dependence has exactly the opposite qualitative effect with respect to the cos θ S dependence of the partonic angular distribution g qq (θ S ). It turns out that the PDF modulation effect quantitatively dominates, and dσ/d cos θ * has a bell (concave) shape (see Fig. 4) rather than the inverse-bell (convex) shape of the angular distribution g qq (θ S ) of the underlying partonic process. We note that, varying cos θ S over a wide range around the central region (say, |cos θ S | ∼ < 0.5), the lower limit on the PDF integration varies in a restricted range this implies that the results of Fig. 4 for dσ/d cos θ * in the region |cos θ S | ∼ < 0.5 are quite sensitive to the behaviour of the corresponding double differential cross section in a very limited range (50 GeV< M γγ ∼ < 60 GeV) of M γγ . Figure 5: The cos θ * differential cross section for E T max = 10 GeV and with the constraint 200 GeV < M γγ < 800 GeV in addition to the same photon kinematical cuts as in Fig. 4. The results, which are obtained with the scales µ R = M γγ /2 and µ F = µ f rag = 2M γγ , are as follows: LO smooth isolation (black dotted), LO standard isolation (light-blue dash-dotted), NLO smooth isolation (red solid) and NLO standard isolation (blue dashed).
From our discussion it follows that the PDF modulation effect on the shape of dσ/d cos θ * can be reduced by applying an additional kinematical cut on M γγ , namely, M γγ > M min γγ with a fixed (θ S independent) value M min γγ , and by selecting sufficiently large values of M min γγ . In particular, increasing M min γγ one can eventually recover the qualitative θ S dependence due to g qq (θ S ). For illustrative purpose, in Fig. 5 we present the results for dσ/d cos θ * with the same kinematical cuts as in Fig. 4-right (E T max = 10 GeV) and the additional constraint 200 GeV< M γγ < 800 GeV (i.e., M min γγ = 200 GeV). We see that the shape of dσ/d cos θ * in Fig. 5 is much different from that in Fig. 4-right and it is qualitatively more similar to the shape of g qq (θ S ). The constraint M γγ < 800 GeV has a negligible quantitative effect on the shape of dσ/d cos θ * . At the LO, the additional constraint M γγ > 200 GeV implies that the p T cuts have no effect on the shape of dσ LO smooth /d cos θ * in the region where |cos θ * | ∼ < 0.97: within this cos θ * region, dσ LO smooth /d cos θ * follows the shape of g qq (θ S ), modulo a PDF effect that is due to the kinematical cuts on the photon rapidities (the effect of the rapidity cuts is discussed below). We can also comment on the diphoton production study that is presented in Ref. [27] (see Sect. III.11 therein). That study uses the kinematical cuts p hard T γ ≥ 40 GeV, p sof t T γ ≥ 30 GeV and 100 GeV< M γγ < 160 GeV, which correspond to M LO dir = 80 GeV and M min γγ = 100 GeV (note that M min γγ is much closer to M LO dir with respect to the cuts considered in Fig. 5). The corresponding differential cross section dσ/d cos θ * (see Fig. III.50 in Ref. [27]) has a maximum value at |cos θ * | ∼ 0.6, and a shape that is somehow intermediate between those in Figs. 4 and 5: this behaviour is in agreement with the expectation from our discussion.
Our discussion and the results in Figs. 4 and 5 evidently show that the shape of dσ/d cos θ * can be strongly affected by the applied kinematical cuts as the consequence of a non-trivial interplay between underlying hard-scattering dynamics and PDF behaviour.
The results in Figs. 3 and 4 are obtained by also including the photon rapidity cut |y γ | < y M (specifically y M = 2.37) in addition to the p T cuts (see Eq. (18)) that we have just discussed. The photon rapidity cut reduces the size of the cross sections but, since the value of y M is sufficiently large, the overall qualitative shape of the differential cross sections is basically unchanged. More precisely, the rapidity cut leads to the LO upper boundaries |cos θ S | < tanh(y M ) (|cos θ S | ∼ < 0.98 in Fig. 4) and M γγ < e −y M √ s (M γγ ∼ < 650 GeV with y M = 2.37 and √ s = 7 TeV) on cos θ S and M γγ , respectively, and it modifies the form of the PDF luminosity contribution in Eq. (14). The modification amounts to the replacement L ab (τ ) → L ab (τ ; y max ), where the customary PDF luminosity L ab (τ ) is replaced by a PDF luminosity with 'rapidity restriction'. The rapidity restricted luminosity L ab (τ ; y max ) is simply obtained by inserting the constraint | ln(x 1 /x 2 )| < 2y max in the {x 1 , x 2 } integration region of Eq. (17) (at the LO, the rapidity y γγ of the diphoton system is 2y γγ = ln(x 1 /x 2 )). The value of y max is related to the photon rapidity cut, y max = y M − |∆y γγ |/2, and it depends on the diphoton rapidity separation ∆y γγ and, hence, on cos θ S (through Eq. (15)). Therefore, the rapidity restriction produces a suppression of the PDF luminosity contribution, and the suppression is larger at larger values of |cos θ S |.
We now consider LO kinematical distributions within standard cone isolation. The LO differential cross sections dσ LO standard are obtained by combining the direct and fragmentation components, Many features of the fragmentation component contribution dσ LO frag are similar to those of dσ LO smooth , and we only note the main differences. In the fragmentation component the photon is accompanied by collinear hadronic fragments and, therefore, we have ∆Φ γγ = π (as in the case of smooth cone isolation at LO), while p T γγ = 0 (at variance with respect to smooth cone isolation). Moreover, due to the isolation procedure, we have p T γγ < E T max . Since p T γγ = 0, cos θ S (see Eq. (15)) is not exactly equal to the Collins-Soper variable cos θ * : the relation between the two angular variables is cos θ * = M 2 γγ /(M 2 γγ + p 2 T γγ ) cos θ S (which is valid for ∆Φ γγ = π). Since we are considering relatively-small values of E T max , we still have p T γγ ≪ M γγ and cos θ * ≃ cos θ S . The LO double differential cross section for the fragmentation component is where As in the case of Eq. (14), the right-hand side of Eq. (19) does not include an overall proportionality factor of order unity and we have explicitly written only the single-fragmentation contribution due to the initial-state qg partonic channel (the dots in the right-hand side of Eq. (19) stand for all the other single-fragmentation and double-fragmentation contributions). As previously remarked in the context of our discussion of the total cross section σ LO standard , when the fragmentation component is large (i.e., with a similar size as the direct component) the single-fragmentation contribution from the qg initial state gives the bulk of the entire fragmentation component. The angular distribution g qg (θ S ) in Eq. (20) is due to qg → qγ scattering and its cos θ S dependence is similar to that of g qq (θ S ) in Eq. (16).
The function D qg in Eq. (21) is an 'effective' (with isolation) partonic luminosity, which is obtained by convoluting the qg PDF luminosity L qg with the quark-to-photon fragmentation function D q/γ . The boundary value z min in the convolution is due to photon isolation (in the case of the single fragmentation component, z min is due to the isolation requirement z > p sof t T γ /(p sof t T γ + E T max ) ) and it increases as E T max decreases, thus leading to an increasing suppression effect of the fragmentation component as E T max decreases (see Figs. 3 and 4).
Despite the isolation suppression, we have already remarked that the effective partonic luminosity α S α D qg (τ ; z min ) can still be quantitatively similar to L qq (and, hence, dσ LO frag and dσ LO dir = dσ LO smooth can have similar size) because L qg is larger than L qq . Increasing the value of M γγ , the photon transverse momenta and, consequently, z min tend to increase (unless |∆y γγ | and, correspondingly, |cos θ S | have large values), therefore reducing the size of the fragmentation component. The effect is visible in the LO results of Fig. 3, which show that the relative difference between dσ LO smooth and dσ LO standard is reduced at high values of M γγ . The effect is also visible in the comparison between the LO results of Figs. 4-right and 5: the invariant-mass cut M γγ > 200 GeV strongly reduces the relative contribution of the fragmentation component, unless |cos θ * | is large. We note that, increasing the value of M γγ , the relative effect of the fragmentation component also decreases because L qg and L qq become quantitatively more similar for increasing values, of the parton momentum fraction x.
, both values, p H and p S , of the p T cuts are effective in the case of the fragmentation component. In particular, they still lead to an LO kinematical boundary, The shape of dσ LO frag /dM γγ near the LO threshold is qualitatively similar to that of dσ LO smooth /dM γγ : the maximum value and the sharp decrease of dσ LO frag /dM γγ for M γγ ∼ M LO frag are produced by the p T cuts through the same kinematical mechanism that we have described in the case of the smooth isolation result. In the case of the single-fragmentation component, using 2p sof t T γ = 2zp hard T γ = √ zM γγ sin θ S , we can express z min as a function of E T max /(M γγ sin θ S ) and the p T cuts lead to the constraints M γγ sin θ S > 2p H √ z and √ zM γγ sin θ S > 2p S . Note that the integration region over the photon momentum fraction z is limited also by the effect of the p T cuts. In the vicinity of the LO invariant-mass threshold, M γγ ∼ M LO frag , the phase space integration region over z is strongly suppressed by these cuts, and the M γγ distribution vanishes proportionally to (M γγ − M LO frag ) 3/2 . In particular, this vanishing behaviour is stronger and smoother (by a factor of M γγ − M LO frag ) than the LO vanishing behaviour of dσ/dM γγ for smooth cone isolation.
The effect of the rapidity cut |y γ | < y M is analogous to the case of smooth isolation: it leads to the replacement L qg (τ ) → L qg (τ ; y max ) (y max = y M − |∆y γγ |/2) in the right-hand side of Eq. (21). Figure 6: The M γγ differential cross section for E T max = 10 GeV with the same photon kinematical cuts as in Fig. 3-right and an increased cut on the minimum value of p hard T γ (p hard T γ > 32 GeV). The scale variation bands of the LO and NLO results for smooth and standard isolation are as follows: LO smooth isolation (black dashed), LO standard isolation (light-blue dashed), NLO smooth isolation (red solid) and NLO standard isolation (blue solid).

Differential cross sections at the NLO
We now move to discuss the NLO results for the differential cross sections. In addition to the M γγ and cos θ * distributions, we present the NLO results for the differential cross sections with respect to the diphoton azimuthal separation ∆Φ γγ (Fig. 7) and to the transverse momentum p T γγ of the photon pair (Fig. 8). The NLO results in Figs. 7 and 8 are obtained by using the reference kinematical cuts described at the beginning of Sect. 2.3 (as in the case of Figs. 3 and 4). As we have previously noticed, the LO calculation leads to non-vanishing differential cross section only in specific LO kinematical subregions. Therefore, outside these LO kinematical subregions (i.e., if M γγ < M LO , ∆Φ γγ = π or p T γγ > E T max ), the NLO results presented in Figs. 3, 6, 7 and 8 actually represent 'effective' LO predictions for the corresponding differential cross sections. We also note that, dealing with 'effective' LO predictions, the distinction between direct and fragmentation components is unambiguous in the context of standard cone isolation.
We first discuss the results for the invariant-mass distribution (Fig. 3). It is convenient to consider three different regions: the region of intermediate values of M γγ (say, 45 GeV< M γγ < 65 GeV) around the LO kinematical threshold at M γγ ∼ M LO , and the regions of higher and lower values of M γγ . For the purposes of the subsequent discussions, we also define Note that M LO is equal to (or smaller than) the minimum between the thresholds M LO frag and M LO dir .
In the high-mass region where M γγ > 65 GeV (Fig. 3), the NLO results for smooth and standard isolation are quantitatively very similar, with a scale dependence that is comparable to that of the corresponding NLO total cross sections. The NLO corrections are large for both isolation criteria and for both values of E T max = 2 GeV and 10 GeV considered in Fig. 3. All these features are similar (both qualitatively and quantitatively) to those of the NLO and LO total cross sections and they have exactly the same origin, which we have already remarked in our discussion on the total cross sections. We do not repeat such a discussion on the role of the qg initial-state channel and of the corresponding PDF luminosity at different perturbative orders and within the two different isolation criteria. The NLO total cross section receives a negligible contribution from dσ N LO /dM γγ in the lowmass region (M γγ < 45 GeV): the contribution is smaller than the scale dependence of the total cross section. The LO kinematical boundary on M γγ is unphysical: it is due to the p T cuts on the photons but it is an artifact of the LO kinematics, which implies ∆Φ γγ = π. Physical diphoton events (and also corresponding partonic contributions beyond the LO) can have ∆Φ γγ < π: they produce non-vanishing values of dσ/dM γγ at M γγ < M LO , although this kinematical region is strongly suppressed by the photon p T cuts. Owing to energy conservation and the presence of the p T cuts, the low-mass region selects diphoton events with small values of ∆Φ γγ . Owing to transverse-momentum conservation, in the low-mass region these p T cuts effectively act also as a lower limit on p T γγ or, equivalently, on the total transverse momentum of the hadronic (partonic) final-state system. Roughly speaking, low values of M γγ imply small values of ∆Φ γγ and, in turn, relatively-large values of p T γγ .
In general, kinematics leads to the minimal constraint which implies that decreasing values of M γγ necessarily require decreasing values of ∆Φ γγ . The constraint in Eq. (23) is obtained by setting ∆y γγ = 0; larger values of |∆y γγ | further reduce the value of ∆Φ γγ at fixed value of M γγ . Eventually the kinematical lower limit ‡ on M γγ is obtained by the combined effect of the p T cuts and the cut on the minimum angular separation, ‡ If ∆Φ γγ < R min γγ , kinematics leads to the replacement sin(∆Φ γγ /2) → sinh( (R min γγ /2) 2 − (∆Φ γγ /2) 2 ) in the right-hand side of Eq. (23).
R γγ > R min γγ = 0.4, between the photons. Since the value of R min γγ is small, the lower limit on M γγ is M γγ ∼ > M LO R min γγ /2 ≃ 10 GeV.
Using kinematical considerations, the region of small values of ∆Φ γγ can also be more directly related to p T γγ . As a consequence of transverse-momentum conservation and of the photon p T cuts, provided ∆Φ γγ < π/2 we have Therefore, if p H = 25 GeV and p S = 22 GeV, the region where ∆Φ γγ ≃ 0 does not contribute to the p T γγ spectrum unless p T γγ ∼ > 47 GeV.
In the low-mass region (Fig. 3), the scale dependence of the NLO differential cross section is larger than the corresponding dependence of the NLO total cross section, as expected from an effective LO prediction. In the case of smooth isolation, the scale dependence slightly increases by decreasing M γγ ; at M γγ ∼ 20 GeV the scale dependence is roughly a factor of 2 larger than the scale dependence of the NLO total cross section. The scale dependence of the standard isolation result is larger, and it increases by either decreasing M γγ or increasing E T max ; at M γγ ∼ 20 GeV and E T max = 10 GeV the scale dependence of the NLO differential cross section is roughly a factor of 3.6 larger than the scale dependence of the NLO total cross section. At the lower value of E T max (2 GeV), smooth and standard isolations give similar NLO differential cross sections, within the corresponding scale variation uncertainties. At the higher value of E T max (10 GeV), the smooth isolation result is systematically smaller than the standard isolation result and the relative difference increases by decreasing M γγ : the NLO results for standard isolation is roughly a factor of 3.8 (2.5) larger than the corresponding result for smooth isolation at M γγ ∼ 20 GeV (M γγ ∼ 30 GeV).
The observed NLO differences between smooth and standard isolation in the low-mass region (analogously to the corresponding LO differences in the high-mass region and to the LO differences of the total cross section) deserve specific comments. The NLO calculation for smooth isolation has two photons and one parton in the final state. Owing to transverse-momentum conservation, the parton can be inside the photon isolation cones only if ∆Φ γγ > π −R ≃ 2.7, which corresponds to M γγ ∼ > 46 GeV in view of the constraint in Eq. (23). Therefore, in the entire low-mass region the NLO result for smooth isolation is exactly independent of E T max , and that represents a much simplified approximation of the expected physical behaviour. The independence of E T max also implies that the smooth isolation result is exactly equal to the result of the direct component for standard isolation. Therefore, the observed differences between the two isolation criteria are entirely due to the fragmentation component of the standard isolation calculation. The NLO result for smooth isolation (or, equivalently, for the direct component) is due to the partonic processes qq → gγγ and qg → qγγ (orqg →qγγ), where the collimated diphoton system recoils against the final-state parton, and the initial-state qg process gives the dominant contribution because of the larger L qg PDF luminosity. The large NLO effect of the fragmentation component (especially at the higher value of E T max , which leads to a smaller suppression effect from isolation) is due to its numerous partonic processes, which, moreover, also include the gg and qq initial states: the corresponding partonic cross sections (although they are suppressed by isolation) can be enhanced by the size of the PDF luminosity (L gg and L qg have a comparable size). In particular, two of these partonic processes are and where the final-state quark (or the antiquark in the case of the gg channel of Eq. (27)) is collimated with the γ and fragments into a second photon. At low values of M γγ these two partonic processes are enhanced by the relative factor (M LO /M γγ ) 2 (its value is about 2.4 and 5.5 at M γγ ∼ 30 GeV and M γγ ∼ 20 GeV, respectively), which originates from the final-state perturbative singularity in the qγ collinear limit (orqγ collinear limit in the case of the gg channel of Eq. (27)). Analogous processes, namely, qg → gqγγ and gg →qqγγ , where a final-state quark (or antiquark) is inside the isolation cone of one photon, contribute (and their contribution depends on E T max ) to the NNLO calculation for smooth isolation. We thus expect that these processes lead to large radiative corrections § (this expectation is confirmed by the NNLO results presented and discussed in Sect. 3.2; see, in particular, Fig. 11-left) and that the corresponding NNLO result removes the large differences with respect to the standard isolation result.
In summary, according to our discussion, the sizeable NLO differences between standard and smooth isolation results that are observed in the low-mass region (Fig. 3) are more an artifact of the NLO calculation than a physical effect due to the two different isolation criteria. Certainly, the achievement of high-precision QCD predictions in the low-mass region is challenging, basically because of the relatively-low value of the characteristic hard-scattering scale M γγ .
The region of intermediate values of M γγ (45 GeV < M γγ < 65 GeV in Fig. 3) is the region where the shape of dσ/dM γγ is directly most sensitive to the p T cuts. The NLO results for smooth and standard isolation are quite similar in this region. Independently of the isolation criterion, the differential cross section starts to rapidly increase at M γγ ≃ M LO ≃ 47 GeV and the position of the peak (the maximum of dσ/dM γγ ) is basically unchanged with respect to that of the LO result for standard isolation. It is noticeable that dσ/dM γγ starts to increase at M γγ ≃ M LO (M LO ≃ M LO frag ) despite the presence of two different LO threshold, M LO frag and M LO dir , that are displaced (the LO § At the strictly formal level, we note that the enhancing factor (M LO /M γγ ) 2 in the NLO calculation for standard isolation is 'unphysical' in the limit M γγ → 0, and the (M γγ ) −2 behaviour is softened by higher-order radiative corrections for both standard and smooth isolation. threshold at M LO dir is somehow more 'unphysical': it is insensitive to the value of p S since absolutely no partonic transverse momentum is allowed in the final state of the corresponding LO calculation). The NLO corrections are obviously large for M γγ < M LO dir = 50 GeV, and they are sizeable also at M γγ > 50 GeV (for the same reason as for the high-mass region). In the case with E T max = 10 GeV we note that the NLO result for smooth isolation tends to be larger than the corresponding result for standard isolation, in disagreement with the physical constraint in Eq. (5). We do not regard this behaviour as particularly worrisome since the values of dσ/dM γγ for the two isolation criteria are basically consistent with each other within the computed scale variation uncertainties ¶ . Moreover, these scale variation effects are expected to underestimate higher-order perturbative uncertainties. Indeed, perturbative calculations in regions around unphysical fixedorder thresholds are known [60] to be generally affected by perturbative instabilities at higher orders (further comments on this point are postponed to the end of this section). In our specific case, dσ/dM γγ is very steep in the region around the LO kinematical threshold and even the effect of little instabilities is amplified by the large slope of dσ/dM γγ . The ensuing effect on the comparison between smooth and standard isolation results can be relevant since the slopes of the two results are both large and they are different. We also add that the standard isolation results do not include the uncertainty (which is difficult to be quantified) that is due to the limited knowledge of the photon fragmentation functions (an increased value of the fragmentation functions can reduce the differences between smooth and standard isolation results in this M γγ region).
We briefly comment on the results in Fig. 6, which are analogous to those in Fig. 3-right apart from having more asymmetric p T cuts. The value of p H increases in going from Fig. 3-right to  Fig. 3). This similarity between the NLO results also confirms our observation in the previous paragraph that the 'threshold' at M γγ ∼ M LO is somehow more 'physical'. Even in the case of more asymmetric cuts, the NLO corrections are able to remove the relative 'deficit' of the LO results for smooth isolation in the vicinity of the LO threshold. Comparing Figs. 3-right and 6 we see that the slope of dσ/dM γγ is smaller in the case of more asymmetric cuts, and the NLO results for smooth and standard isolation are more similar (in particular, the smooth isolation result does not tend to be larger than the standard isolation result, consistently with our previous comment on the behaviour observed in Fig. 3). Additional comments on symmetric and asymmetric p T cuts are presented in Sect. 3.
In Fig. 7 we present the NLO results for the differential cross section with respect to the azimuthal angle separation ∆Φ γγ of the two photons. At the LO, ∆Φ γγ = π for both isolation criteria. As is well known [28], at higher perturbative orders the computation of the differential cross section is affected by large logarithmic corrections in the region near ∆Φ γγ = π. In this ¶ We also note that in this M γγ region the scale dependence of the NLO standard isolation result at E T max = 10 GeV is quite small, in particular, if it is compared with the corresponding scale dependence at E T max = 2 GeV. Such a small NLO scale dependence seems accidental. Figure 7: The NLO results (scale variation bands) for the ∆Φ γγ differential cross section that are obtained by using the smooth (red solid band) and standard (blue dashed band) cone isolation criteria with E T max = 2 GeV (left panel) and E T max = 10 GeV (right panel). The photon kinematical cuts are the same as in Fig. 3.
region, any fixed-order QCD result is physically not reliable, and reliable quantitative predictions for the detailed shape of the ∆Φ γγ distribution can be recovered only through all-order perturbative resummation [42,43] of these large logarithmic contributions. Owing to this reason, in the results of Fig. 7 we have excluded the region around ∆Φ γγ = π. This also implies that the NLO results in Fig. 7 (with ∆Φ γγ < π) actually represent 'effective' LO predictions for the ∆Φ γγ distribution.
Independently of the value of E T max , we note that the ∆Φ γγ distribution is sharply peaked at large values of ∆Φ γγ : the cross section decreases by more than one order of magnitude by decreasing ∆Φ γγ from ∆Φ γγ ∼ 3 to ∆Φ γγ ∼ 2.5, and it still decreases by about one order of magnitude in going from ∆Φ γγ ∼ 2.5 to ∆Φ γγ ∼ 0.5. Standard and smooth isolation results have a similar scale dependence that increases from roughly ±20% at ∆Φ γγ ∼ 2.5 to roughly ±30% at ∆Φ γγ ∼ 0.5 (in the case of standard isolation the scale dependence also slightly increases by increasing the value of E T max ). As expected from an effective LO prediction, this scale dependence is larger (by about a factor of 2-3) than the scale dependence of the NLO total cross section (the bulk of the NLO total cross section is due to the region near ∆Φ γγ ∼ π, where the scale dependence of the ∆Φ γγ differential cross section is reduced).
The main quantitative differences between standard and smooth isolation results appear by examining the dependence on E T max of the ∆Φ γγ distribution. At very small values of E T max (E T max = 2 GeV) the two isolation criteria lead to very similar results, within the corresponding scale variation uncertainties. At larger values of E T max (E T max = 10 GeV), smooth cone results tend to be smaller than standard cone results, and the differences increase by decreasing the value of ∆Φ γγ . At very small values of ∆Φ γγ , the differences are sizeable; for instance, standard NLO results are roughly a factor of 2.1 larger than smooth NLO results at ∆Φ γγ ∼ 0.5. In the case of smooth isolation, the qg initial-state partonic channel contributes more than the qq channel (because of the larger L qg PDF luminosity). At high values of ∆Φ γγ the smooth isolation result depends very weakly on E T max and it is very similar (though it is not exactly equal) to the direct component of the standard isolation result. As previously mentioned, if ∆Φ γγ < π − R ≃ 2.7 the smooth isolation result is exactly independent of E T max and it coincides with the direct component of the standard isolation cross section. Therefore, at low values of ∆Φ γγ the NLO differences between the two isolation criteria are entirely due to the fragmentation component of the standard isolation calculation. Figure 8: The NLO results (scale variation bands) for the p T γγ differential cross section that are obtained by using the smooth (red solid band) and standard (blue dashed band) cone isolation criteria with E T max = 2 GeV (left panel) and E T max = 10 GeV (right panel). The photon kinematical cuts are the same as in Fig. 3.
In Fig. 8 we present the NLO results for the differential cross section with respect to the transverse momentum p T γγ of the photon pair. At LO, as already mentioned, p T γγ = 0 for smooth isolation while p T γγ < E T max for standard isolation. Analogously to the case of dσ/d∆Φ γγ at ∆Φ γγ ≃ π, the perturbative computation of dσ/dp T γγ at small values of p T γγ is affected by large logarithmic corrections [28] that have to be treated by an all-order resummation procedure [42,43] to obtain reliable QCD predictions. Therefore, in Fig. 8 we do not show the NLO results at small values of p T γγ . Actually, we consider only the region with p T γγ > E T max , where the results have to be regarded as 'effective' LO predictions.
The scale dependence of the NLO result for dσ/dp T γγ in the case of smooth isolation is similar to that of the NLO total cross section, while the scale dependence in the case of standard isolation is larger. At small values of E T max (e.g., E T max = 2 GeV) the NLO results for smooth and standard isolation are very similar, within the corresponding scale dependence. At larger values of E T max (E T max = 10 GeV) the standard isolation result tends to be larger than the smooth isolation result: however, the ratio of the two results is always smaller than approximately 1.8. In the case of smooth isolation, the qg initial-state channel gives the dominant contribution. The partonic final state of the NLO calculation for smooth isolation includes the two photons and a parton. Owing to transverse-momentum conservation, the transverse momentum of the parton is equal to p T γγ and, therefore, if p T γγ > E T max (as in Fig. 8) the parton is not allowed (by isolation) to be inside the photon isolation cones. Consequently, if p T γγ > E T max the NLO smooth isolation result is independent ‡ of E T max and it exactly coincides with the result of the direct component for ‡ Parton radiation inside the isolation cone is kinematically allowed. If p T γγ > E T max the smooth isolation result is independent of E T max but it still depends on the photon isolation criterion: events with p T γγ > E T max are kinematically allowed but they are rejected by the isolation requirement. standard isolation. Therefore, the differences between the smooth and standard isolation results in Fig. 8 are entirely due to the fragmentation component of the standard isolation result.
We comment on the comparison between the smooth and standard isolation results in Figs. 7 and 8 at large values of E T max (at small values of E T max , the results are very similar since the fragmentation component is much suppressed). As previously discussed, the fragmentation component of the NLO calculation is dynamically enhanced at small values of M γγ (mainly because of the presence of the fragmentation processes in Eqs. (26) and (27)). The low-∆Φ γγ region in Fig. 7 receives contributions from both small and large values of M γγ . The same happens in the large-p T γγ region considered in Fig. 8 and, moreover, at large values of p T γγ even diphoton events with small ∆Φ γγ tend to have larger values of M γγ (at fixed ∆Φ γγ , M γγ kinematically increases by increasing p T γγ ). Therefore, the relative enhancement effect of the NLO fragmentation component at low values of M γγ continuously decreases in going to the low-∆Φ γγ region and to the large-p T γγ region. As a consequence of this reasoning, the relative difference between standard and smooth isolation results correspondingly and continuously decreases in going from dσ/dM γγ at small M γγ (Fig. 3-right), to dσ/d∆Φ γγ (Fig. 7-right) and to dσ/dp T γγ (Fig. 8-right): this behaviour is indeed observed in the NLO results of the corresponding figures. Since the partonic processes that enhance the NLO fragmentation component contribute to the smooth isolation result at NNLO (as discussed in the accompanying comments to Eqs. (26)-(29)), we expect NNLO corrections for smooth isolation that are large in the low-M γγ region and that continuously decrease in going to small values of ∆Φ γγ , to large values of p T γγ and to high values of M γγ . The actual NNLO results that are presented in Sect. 3.2 basically confirm the expectation from our discussion.
In the context of standard cone isolation, the relevance of the fragmentation processes in Eqs. (26) and (27) for the differential cross sections dσ/d∆Φ γγ and dσ/dp T γγ was already remarked in Ref. [95]. In the case of dσ/d∆Φ γγ at small values of ∆Φ γγ , the authors of Ref. [95] noticed that these fragmentation processes give a NLO contribution that is sizeably larger than that of the direct component of the NLO cross section, as we have also remarked. The authors of Ref. [95] also nicely discussed how these fragmentation processes have impact on dσ/dp T γγ at relativelylarge values of p T γγ , thus producing a shoulder-type shape of the p T γγ distribution. The p T γγ shoulder is clearly visible in the standard isolation results with E T max = 10 GeV in Fig. 8-right (the p T γγ shoulder is much less visible in Fig. 8-left, since at very small values of E T max , such as E T max = 2 GeV, the fragmentation component is highly suppressed). The NLO shape of dσ/dp T γγ for standard isolation changes (it flattens out) in the region where 40 GeV ∼ < p T γγ ∼ < 50 GeV, and dσ/dp T γγ is larger than the corresponding NLO smooth isolation result (which is exactly equal to the direct component of the standard isolation result at NLO) for increasing values of p T γγ . According to our previous discussion, the smooth isolation NNLO processes in Eqs. (28) and (29) have a dynamical role which is analogous to that of the NLO fragmentation processes in Eqs. (26) and (27) (and to part of their NNLO radiative corrections). Therefore, we expect a pronounced p T γγ shoulder in the smooth isolation results at NNLO. This expectation is confirmed by the results that we present and further discuss in Sect. 3.2.
We now turn to consider the differential cross section with respect to the polar angle θ * . Figure 4 presents the NLO results by using the customary kinematical cuts of this subsection. Since the LO result covers the entire range of cos θ * , the corresponding NLO results are effectively NLO QCD predictions, and the main kinematical effects that we have discussed at the LO level are unchanged at the NLO level.
The main features of the comparison between the LO and NLO results in Fig. 4 are very similar to those of the comparison between the corresponding total cross sections. At variance with other kinematical distributions, the shape of the NLO results for dσ/d cos θ * has very little dependence on E T max : the increase of E T max from 2 GeV to 10 GeV has mainly an effect on the overall normalization, whose size increases analogously to the value of the corresponding NLO total cross section. The scale dependence of the NLO results is very similar for the two isolation criteria: its value is approximately ±15% at cos θ * = 0, and it is slightly smaller at larger values of |cos θ * |. Smooth and standard isolation results at NLO are very similar (with overlapping scale variation bands) for |cos θ * | ∼ < 0.7. At larger values of |cos θ * | the NLO result for standard isolation is systematically larger than the corresponding result for smooth isolation: the ratio between the standard and smooth results is approximately 1.3 at |cos θ * | ∼ 0.9. In the region of central values of cos θ * (|cos θ * | ∼ < 0.5), we note that the NLO result for smooth isolation tends to be larger than the corresponding result for standard isolation, although the two NLO results are consistent with each other within the computed scale variation dependence. A similar tendency has been already noticed in the case of the differential cross section dσ/dM γγ in the region where M γγ ∼ M LO (Fig. 3). The two effects, for dσ/d cos θ * and dσ/dM γγ , are certainly related since, as a consequence of the photon p T cuts (as previously discussed in Sect. 2.3.2), the differential cross section dσ/d cos θ * in the region where |cos θ * | ∼ < 0.5 is very sensitive to the M γγ dependence of dσ/dM γγ d cos θ * in the region 50 GeV ∼ < M γγ ∼ < 60 GeV.
We also comment on the cos θ * distribution in the region of large values of |cos θ * |. Comparing the NLO calculations for the standard and smooth isolation criteria, standard isolation involves many more partonic processes. In particular, the single-fragmentation component receives contribution from the partonic processes and qq → qqγ (q → γ + X) , in which the photon and a final-state fermion can have a small relative angle and the other final-state fermion fragments into a second photon. At large values of |cos θ * | (which roughly correspond to large rapidity separations |∆y γγ |), these two processes are dominated by the effect of the exchange of one gluon in the t-channel of the 2 → 2 fermion scattering subprocess: the gluon exchange leads to an angular distribution that is proportional to (1 − cos 2 θ * ) −2 . Therefore, the gluon exchange subprocesses are dynamically enhanced by the relative factor (1 − cos 2 θ * ) −1 with respect to the fermion exchange subprocesses (see, e.g., Eqs. (16) and (20)) that contribute to the NLO calculation for smooth isolation. Although the gluon exchange processes are suppressed by the isolation requirements, their effect on the NLO results for standard isolation is not negligible at large values of |cos θ * |. We also note that corresponding gluon exchange processes enter the calculation for smooth isolation at the NNLO. Such processes are and qq → qqγγ , in which the two photons are produced at large rapidity separation and each photon is at small relative angle with respect to one final-state fermion. These processes can enhance the size of the NNLO correction within smooth isolation at large values of |cos θ * |.
We briefly comment on the NLO results for dσ/d cos θ * in Fig. 5, which are obtained by applying the additional kinematical cut 200 GeV< M γγ < 800 GeV. For the sake of simplicity, in Fig. 5 we present LO and NLO results without considering scale variation dependence and we simply use µ R = M γγ /2 and µ F = µ f rag = 2M γγ . As already remarked in our discussion of the LO results, the constraint M γγ > 200 GeV has a major effect on the shape of dσ/d cos θ * , which qualitatively follows the shape of the angular distribution of the underlying partonic processes. This behaviour is not affected by including NLO corrections. Comparing the results in Fig. 4-right and 5, we see that the high-mass constraint M γγ > 200 GeV sizeably (and obviously) reduces the values of the cross sections and it also modifies the size of the NLO corrections. The NLO results for smooth and standard isolation are very similar also at large values of |cos θ * | (e.g., |cos θ * | ∼ 0.9). At central values of |cos θ * |, the NLO result for standard isolation is (slightly) higher than the NLO result for smooth isolation: this behaviour is different from that in the NLO results of Fig. 4 and it is consistent with our discussion on the relevance of the region of intermediate values of M γγ for the behaviour of the NLO results in Fig. 4. Table 1 and Figs. 2-4, standard and smooth isolations lead to QCD results in good quantitative agreement for physical observables that are effectively computed up to the NLO. In view of this agreement we present some additional investigations on the role of the fragmentation component in the QCD computation for the standard isolation criterion. A perturbative scheme (approximation) that is sometimes used in photon isolation computations at the NLO (see, e.g., Refs. [31,45,46]) consists in combining the evaluation of the direct component at the NLO with that of the fragmentation component at the LO. This has the practical advantage of avoiding the computation of the more cumbersome photon fragmentation subprocesses at the NLO. Within this context, one can 'equivalently' [46] use LO and NLO parton-to-photon fragmentation functions. Such a scheme is applied in the diphoton production calculation of Ref. [31], which we use (as implemented in the MCFM program) in our subsequent numerical investigation (we do not use the DIPHOX program, by removing the NLO corrections to the fragmentation component, because LO fragmentation functions are not readily available in the default setup of the DIPHOX code).

As shown in
We consider standard cone isolation and the photon kinematical cuts used in Table 1 and Figs. 2-4. The QCD calculation that includes NLO direct+LO fragmentation components is carried out by using the MCFM program [31], and we use either the NLO fragmentation functions of the BFG set II [92] (as in our customary NLO calculations) or the LO fragmentation functions of the GdRG LO set [96]. At the central value of the scales and using BFG fragmentation functions, we obtain the following values of total cross sections: σ M CF M standard = 29.66 pb with E T max = 2 GeV, and σ M CF M standard = 28.46 pb with E T max = 10 GeV (the GdRG LO set of fragmentation functions leads to results that are larger by approximately 1 pb). The corresponding results for dσ/d cos θ * are presented in Fig. 9-left. We note that these results are similar and consistent with the complete NLO results, within the corresponding scale variation dependence. This implies that the NLO corrections to the fragmentation component are not particularly sizeable. However, we also note that σ M CF M standard and σ N LO standard differ by about 15% at E T max = 10 GeV and, moreover, σ M CF M standard decreases by increasing the value of E T max from 2 GeV to 10 GeV (we have checked that the value of σ M CF M standard at E T max = 4 GeV is intermediate between the values at E T max = 2 GeV and E T max = 10 GeV). This E T max dependence of the MCFM results, which occurs for both the total cross section and dσ/d cos θ * at fixed values of cos θ * (see Fig. 9-left), violates the expected physical behaviour (see Eq. (6)). The use of BFG or GdRG LO fragmentation functions does not change the qualitative dependence on E T max . In contrast, the complete NLO result for standard isolation (see Table 1, Fig. 4 and Fig. 9-right) and also the NLO result for smooth isolation (see Table 1 and Fig. 4) have the expected physical dependence on E T max (total and differential cross sections increase by increasing E T max ). Figure 9: The cos θ * differential cross section for standard cone isolation with two different values of E T max (2 GeV and 10 GeV) and the same photon kinematical cuts as in Fig. 4. The QCD results are obtained at the central value of the scales (µ F = µ R = µ f rag = µ 0 ≡ M γγ ). The results with NLO direct + LO fragmentation components (left panel) use BFG and GdRG LO fragmentation functions. The NLO results (right panel) use BFG fragmentation functions.
We interpret the unphysical dependence on E T max of the MCFM results as the effect of a mismatch between the perturbative orders in the direct and fragmentation components. Indeed the E T max dependence of the MCFM results is mostly produced by the NLO corrections to the direct component through the partonic subprocess qg → qγγ at the tree level. The tree-level contribution of this partonic subprocess is formally positive definite and, as such, it would lead to an increasing (physical) dependence on E T max by increasing E T max . However, this contribution is divergent in the phase space region where the final-state q is collinear to one of the produced photons. The collinear divergence (which is in turn absorbed and factorized in the non-perturbative definition of the quark-to-photon fragmentation function) is removed through its subtraction (which is actually performed in the MS factorization scheme [31]) from the direct component contribution at NLO. Since the collinear-divergent term that is subtracted is formally positive definite, the final NLO correction to the direct component is finite but it is not necessarily positive definite and, consequently, it can have an unphysical dependence on E T max (this unphysical dependence is actually also visible in the direct component of the NLO results in Table 1). It turns out that the result of the calculation with NLO direct + LO fragmentation components has an unphysical dependence on E T max .
In the case of the complete NLO result for standard isolation, the NLO corrections to the fragmentation component (corrections to both the partonic cross sections and the scale evolution of the fragmentation functions) are consistently (at the formal level) included within the MS factorization scheme. Although these corrections are not particularly sizeable in absolute terms, their inclusion leads to NLO results with a qualitative dependence on E T max that agrees with physical expectation.
In the context of smooth isolation, the entire E T max dependence at NLO is due to partonic subprocesses with two photons and a parton in the final state, which are evaluated at the tree level. The contribution of these tree-level subprocesses is positive definite (with no collinear divergences to be subtracted) and it cannot produce any unphysical dependence on E T max : NLO smooth isolation cross sections cannot decrease by increasing E T max . Incidentally, we note that this elementary reasoning, based on positivity, cannot be applied at the NNLO, since the E T max dependence at NNLO receives contributions from both tree-level and loop-level partonic processes. Nonetheless, as shown in Sect. 3.3 (see Fig. 16 and the accompanying comments) the NNLO results for smooth isolation do not show unphysical dependence on E T max .
As anticipated in our comments on the results in Fig. 3, we present a more detailed discussion on the differential cross section dσ/dM γγ in the region where M γγ is close to the LO threshold (M γγ ∼ M LO ). The numerical results in Fig. 3 are obtained by using a (constant) bin size of 2 GeV in M γγ . To quantitatively examine the detailed shape of dσ/dM γγ , we perform the numerical calculation with a finer resolution in M γγ and we use M γγ bins with a constant size of 0.1 GeV (which is 20 times smaller than that used in Fig. 3). The results of the calculation for smooth isolation at both LO and NLO are presented in Fig. 10. Specifically, we use E T max = 10 GeV and, at each perturbative order, we report the two results with the scale choices {µ R = M γγ /2, µ F = 2M γγ } and {µ R = 2M γγ , µ F = M γγ /2} (the region enclosed by these two scale-dependent results corresponds to the scale variation band in Fig. 3). We remark that the smooth isolation results in Fig. 10 exactly refer to the same quantity (and to the same theoretical setup) as the corresponding results in Fig. 3-right, the only difference being the much smaller M γγ bin size used in Fig. 10. Various shape details that are visible in Fig. 10   We note that the qualitative shape of dσ/dM γγ in Fig. 10 is independent of the scale choice, which only affects the size of dσ/dM γγ . The behaviour of the LO result has been discussed in the paragraph below Eq. (18). The LO result has a threshold at M γγ = M LO dir = 50 GeV and, close to the threshold, it vanishes as where ǫ M = |1 − (M γγ /M LO dir )|. The threshold and the square-root behaviours near threshold are visible in the LO results of Fig. 10. At the NLO, dσ/dM γγ is not vanishing at M γγ = M LO dir (the threshold disappears) and it has a cusp behaviour in the vicinity of the LO threshold. Following the general analysis in Ref. [60], we have examined the NLO shape of dσ/dM γγ at M γγ ∼ M LO dir in analytic form and we find the dominant behaviour Considering the NLO result with the scale configuration {µ R = M γγ /2, µ F = 2M γγ } in Fig. 10, we also notice that the value of dσ/dM γγ at M γγ = M LO dir (i.e., the height of the cusp) is quite large and, in particular, it is larger than the value in the peak region at M γγ ∼ 57 GeV. Independently of the scale configuration, we observe that the cusp behaviour is located in a tiny region of M γγ around M γγ = M LO dir and, consequently, such behaviour quantitatively disappears by increasing the bin size in M γγ (see Fig. 3).
The shape of dσ N LO /dM γγ at M γγ ∼ M LO dir (Fig. 10 and Eq. (35)) is definitely unphysical and it deserves additional comments. This shape (and the expression in Eq. (35)) follows from the general discussion and results of Ref. [60]. According to Ref. [60], an observable that has a discontinuity at some point x = x 0 (inside the physical region) at the LO necessarily has a logarithmic divergent (though integrable) discontinuity at the same point at the NLO. In our specific case, the slope of dσ LO /dM γγ has a discontinuity at M γγ = M LO dir (though dσ LO /dM γγ is continuous at the same point): therefore, the slope of dσ N LO /dM γγ has a logarithmically-enhanced discontinuity (see Eqs. (34) and (35)) at M γγ = M LO dir (though dσ N LO /dM γγ remains continuous at the same point).
The logarithmic enhancement of the discontinuity is due to soft-gluon radiation at the NLO [60]. In our specific case, the relevant NLO partonic processes are the real emission contribution qq → γγg (36) at the tree level (here, the final-state gluon is soft and collinear to one of the initial-state partons), and the virtual contribution to qq → γγ at the one-loop level. A non-smooth kinematical mismatch (such as that produced by a discontinuous observable) between the real and virtual contributions produces the logarithmic enhancement [60]. In our specific case, the NLO virtual process qq → γγ still fulfils the same kinematical constraints as at the LO (in particular, transverse-momentum conservation implies p sof t T γ = p hard T γ ≥ p H = 25 GeV, independently of the value of p S ) and, consequently, it only contributes at M γγ ≥ M LO dir . In contrast, the soft gluon radiated in the process of Eq. (36) produces diphoton events with M γγ < M LO dir , since the softer photon absorbs the softgluon momentum recoil ‡ thus decreasing its transverse momentum p sof t T γ below its LO kinematical limit (i.e., we have p sof t T γ < p H ≤ p hard T γ , although p sof t T γ ≥ p S ). It follows that real soft-gluon radiation is completely unbalanced by virtual radiation in the region just below the LO threshold, and this produces the corresponding logarithmically-enhanced cusp at M γγ < M LO dir in Fig. 10 and Eq. (35). Therefore, the shape of dσ N LO /dM γγ just below the LO threshold is exactly a consequence of the reasoning and results of Ref. [60]. Instead, the logarithmic enhancement of the slope of dσ N LO /dM γγ just above the LO threshold (M γγ > M LO dir ) can be somehow regarded as a corollary to the reasoning in Ref. [60]. Both real and virtual terms contribute above the threshold, but here the slope of dσ/dM γγ is already divergent (see Eq. (34)), and not only discontinuous, at the LO: this LO divergent behaviour produces a strong real-virtual mismatch and an ensuing logarithmic enhancement at the NLO (see Fig. 10  The unphysical shape (i.e., the cusp behaviour) of dσ/dM γγ at NLO persists at each subsequent perturbative order. Such an unphysical fixed-order behaviour can be removed by a proper allorder resummation of soft-gluon effects [60]. Resummation leads to a smooth behavior (Sudakov shoulder) [60] of both dσ/dM γγ and its slope from the peak region at M γγ ∼ 57 GeV to the 'approximate' threshold at M γγ ∼ 47 GeV.
In the context of fixed-order computations, the unphysical behaviour produces perturbative instabilities that can be reduced only by considering observables that are sensitive to a 'sufficiently- ‡ Note that such recoil is forbidden in the case of symmetric p T cuts with p H = p S (see a related discussion in Sect. 3.5).

smeared' region in the vicinity of M γγ ∼ M LO
dir . The degree of 'sufficient' smearing or insensitivity depends on various factors, such as the type of observable, the bin size and the kinematical cuts. We briefly comments on these factors. As for the dependence on the observable, the total cross section has little sensitivity to these instabilities, the differential cross section dσ/d cos θ * at central values of cosθ * has some sensitivity, and the differential cross section dσ/dM γγ in the region close to M LO dir is certainly sensitive. The bin size dependence is obvious at the qualitative level, but it is less obvious at the quantitative level. A sufficiently-large M γγ -bin at M γγ ∼ M LO dir removes the cusp behaviour, but it leads to a binned value of dσ/dM γγ that depends on the bin size and also on the average slope of dσ/dM γγ in the region close to the LO threshold (a larger value of the average slope increases the sensitivity to the bin size). In our discussion on the results in Fig. 3, we have argued that dσ N LO /dM γγ is possibly sensitive to the perturbative instabilities even if the M γγ -bin size is 2 GeV. The dependence on the p T cuts has been pointed out in our comments on the results in Fig. 6, where we have observed that the size of the perturbative instabilities can be affected by the average slope of dσ/dM γγ at M γγ ∼ M LO dir , which directly depends on the difference between M LO dir and M LO (i.e., on the values of the p T cuts p H and p S ). Certainly, the relevance of these perturbative instabilities also depends on the required theoretical accuracy of the QCD calculation. In Sect. 3.5 we present additional results and comments on this kind [60] of instabilities and related observables.

Diphoton production at the LHC and NNLO results
In this section we consider diphoton production in pp collisions at LHC energies and we present perturbative QCD results at the NNLO. We use smooth cone isolation since the NNLO calculation for standard cone isolation has not yet been performed. Within smooth isolation, we also present corresponding results at LO and NLO to directly comment on features of the perturbative QCD expansion.
In Ref. [35] we presented NNLO results for diphoton production at the LHC in diphoton kinematical configurations that are typically used in the context of Higgs boson searches and studies. In particular, we considered the region M γγ < 200 GeV and highly-asymmetric cuts on the photon transverse momenta (i.e., p hard T > 40 GeV, p sof t T > 25 GeV). In the following we consider different kinematical configurations and we discuss various aspects of the NNLO results. Diphoton production results at NNLO and comparisons with LHC data are also presented in Refs. [4,7,8,37]. NNLO results for diphoton production in pp collisions at √ s = 100 TeV are presented in Ref. [97] (see Sect. 8.3 therein).
In our computation the radius of the photon isolation cone is set at the value R = 0.4. We use the smooth isolation function χ(r; R) in Eq. (3) (the same form of the isolation function is used in the NNLO results reported in Refs. [4,7,8]) and the value of the power n is set to n = 1 for most of the results, although we comment on the n dependence of total cross sections and of some differential cross sections.
The QCD results are obtained by using the programs 2γNNLO and MATRIX. The theoretical setup of our calculation is the same as described at the beginning of Sect. 2.3 (in particular, we use the value α = 1/137 of the QED coupling constant and the MMHT 2014 sets [91] of PDFs). The only difference with respect to Sect. 2.3 regards the central value µ 0 of the renormalization (µ R ) and factorization (µ F ) scales. Unlike the case of Sect. 2.3 (where µ 0 = M γγ ), throughout this section we use the dynamical central value µ 0 = M γγ 2 + p 2 T γγ = M T γγ (M T γγ is the transverse mass of the diphoton system), which also depends on the transverse momentum p T γγ of the photon pair. We consider independent scale variations of µ R and µ F within the ranges 0.5 ≤ µ R /µ 0 ≤ 2 and 0.5 ≤ µ F /µ 0 ≤ 2 around the central value µ 0 . Practically, we obtain the results for nine scale configurations (we independently combine µ R /µ 0 = {0.5, 1, 2} and µ F /µ 0 = {0.5, 1, 2}) and we evaluate scale uncertainties by considering the maximum value and minimum value among these results. We have checked that, for most of the computed quantities (including total cross sections), the maximum and minimum values correspond to the scale configurations {µ R = µ 0 /2, µ F = 2µ 0 } and {µ R = 2µ 0 , µ F = µ 0 /2}, respectively.
The bulk of the diphoton cross section is produced at small values of p T γγ (p T γγ ≪ M γγ ). Therefore, to the purpose of computing the total cross section, the choice of the dynamical scale µ 0 = M T γγ leads to results that are basically similar to those obtained with the scale choice µ 0 ∼ M γγ . The dynamical scale µ 0 = M T γγ sizeably differs from the diphoton invariant mass M γγ only at high values of p T γγ (i.e., in the kinematical regions where M γγ ≪ p T γγ ). High values of p T γγ can be reached either in the highly unbalanced regime where p sof t T γ ≪ p hard T γ or in the highly boosted regime where the two photons have comparable transverse momenta (p sof t T γ ∼ p hard T γ ) and small values of both the azimuthal angle separation ∆Φ γγ and the rapidity separation ∆y γγ . In both these regimes, p T γγ is balanced by a recoiling high-p T jet, and the dynamical scale µ 0 = M T γγ ∼ p T γγ parametrically mimics the scale of the invariant mass M γγjet (p T γγ ∼ O(M γγjet )) of the diphoton+jet final-state system, which is the characteristic scale of the underlying hardscattering subprocesses.
In this paper we do not evaluate PDF uncertainties. NLO PDF uncertainties at the LHC are computed in the diphoton studies of Refs. [4,7,8] and combined with scale variation uncertainties (PDF uncertainties are found to be typically smaller than scale uncertainties). In Ref. [7] the PDF uncertainty on the NLO total cross section is explicitly quoted and it amounts to about ±5%. The PDF uncertainties that are computed in Ref. [8] are at the level of ±2%.
As main reference kinematical configuration we consider LHC collisions at the centre-of-mass energy √ s = 7 TeV and we use the kinematical acceptance cuts implemented by the ATLAS Collaboration in the analysis of Ref. [4]. We require p hard T ≥ 25 GeV and p sof t T ≥ 22 GeV, we restrict the rapidity of both photons to the regions |y γ | < 1.37 and 1.52 < |y γ | ≤ 2.37, and the minimum angular separation between the two photons is R min γγ = 0.4. These acceptance cuts coincide with the main reference cuts in Sect. 2.3, apart from the exclusion of the rapidity interval 1.37 < |y γ | < 1.52. The ATLAS data are selected by using standard cone isolation with the isolation parameters R = 0.4 and E T max = 4 GeV. We use smooth cone isolation with R = 0.4 and E T max = 4 GeV. All the results in this section refer to this configuration and to this value of E T max , unless otherwise explicitly stated.

Total cross sections
The results of the LO, NLO and NNLO total cross sections are reported in Table 2.
At the LO the value of the total cross section is σ LO = 9.293 pb +10.9 % −11.9 % (scale), and it does    Table 1), apart from an overall reduction (by roughly 10%) of the cross section values, which is due to the exclusion of the photon rapidity region 1.37 < |y γ | < 1.52 .
The value of the NNLO total cross section with n = 1 is σ N N LO = 39.50 pb +7.0 % −7.8 % (scale). This result is obtained by using MATRIX with the default setup [41]. The value of the NNLO cross section is obtained by the code through an r cut → 0 extrapolation of the r cut dependence for r cut > 0.15%. The systematic uncertainty of the extrapolated result, as explicitly reported in Table 2, is at the ±O(1%) level. Such uncertainty is much smaller than the NNLO perturbative uncertainties, and thus still acceptable to the purpose of the present paper. More accurate results could in principle be obtained with MATRIX by using the option switch qt accuracy=1, which lowers the value of the minimum r cut down to 0.05%. The MATRIX result is in agreement with the corresponding result obtained by using 2γNNLO within the systematical uncertainties.
We have computed the n dependence of the total cross section in the range 0.5 ≤ n ≤ 2 (see Table 2). At the NNLO (NLO) we find that the result with n = 1 increases by about 4% (3%) with n = 0.5 and decreases by about 5% (2%) with n = 2. At both NLO and NNLO we note that this n dependence of the total cross section is monotonic and in qualitative agreement with the physical expectation in Eq. (8). In particular, we point out that the n dependence of the cross section increases (especially at n = 2) in going from NLO to NNLO: the NNLO results are more sensitive to the value of n. We also note that variations of n in the interval 0.5 ≤ n ≤ 2 produce (at both NLO and NNLO) variations of the total cross section that are smaller than those produced by the scale dependence at fixed n. In view of this, we limit ourselves to using n = 1 (unless otherwise stated) for most of our subsequent studies.
Throughout this section we present explicit results on the contribution of the different initialstate partonic channels to various observables at the NNLO. To simplify the presentation we consider the partition of the total result in three contributions: the contribution of the gg partonic channel (and its partial component from the box contribution gg → γγ) ‡ , the contribution of the qg+qg partonic channel (it is dominated by its qg component and it is briefly labelled as qg channel) and the contribution of the remaining partonic channels (this contribution is dominated by the qq initial state and it is briefly labelled as qq channel). This decomposition in partonic channels has a mild scale dependence and we always present results at the central values µ R = µ F = µ 0 = M T γγ of the scales. About 9% of σ N N LO is due to the gg channel. Therefore, the contribution of the gg channel at NNLO is not sizeable (it is quantitatively comparable to the size of the scale dependence of σ N N LO ), despite the fact that the contribution of the box gg → γγ is approximately equal to one half of the LO total cross section. We also note that the total gg contribution at NNLO is partly smaller than its box component: the additional NNLO contributions from gg collisions turn out to be negative and have an absolute size that is approximately one quarter of the box contribution. The NNLO result for the total cross section is dominated by the qg and qq channels, which approximately equally contribute to σ N N LO . The NNLO cross section receives contributions of about 48% from the qq channel and of about 43% from the qg channel. The NLO total cross section is roughly 3 times larger than σ LO . The NNLO K factor, K N N LO = σ N N LO /σ N LO , at central values of the scales is K N N LO ≃ 1.4. We see that both NLO and NNLO corrections are sizeable. The large size of the QCD radiative corrections up to NNLO is justifiable and understandable following and extending the reasoning that we have presented in Sect. 2.3.1 (see Eqs. (9) and (11) and accompanying comments). At the LO only the qq channel contributes. At the NLO, σ N LO receives contribution also from the qg channel. The relative NLO correction from the qg channel is of the order of α S L qg /L qq and the large value of the PDF luminosity ratio L qg /L qq compensates the suppression factor produced by α S . Therefore, the presence of a new NLO partonic channel, the qg channel, with a large PDF luminosity implies that, at the quantitative level, there is no parametric hierarchy of O(α S ) between NLO corrections and LO result. As a consequence, the NLO result has to be quantitatively regarded as an effective lowestorder estimate of the cross section. The next-order corrections to this result are parametrically of O(α S ) (the contribution of the new gg channel at NNLO is not particularly sizeable) and they turn out to have a 'moderate' size. Indeed, the value K N N LO ≃ 1.4 has a size that is not much different from the typical (and expected) size of NLO K factors for various hard-scattering processes at hadron colliders. We also note that the scale dependence of the total cross section is partly reduced in going from NLO to NNLO (it is not so in going from LO to NLO). However, we remark that the values of σ N N LO and σ N LO do not overlap by including their corresponding scale dependence. This implies that the computed scale dependence of σ N N LO cannot be consistently regarded as a reliable estimate of uncalculated radiative corrections to the total cross section. The 'true' theoretical uncertainty of σ N N LO due to higher-order corrections is certainly larger than the NNLO scale dependence that we have computed.

Differential cross sections
We now move to consider kinematical distributions. In particular, we consider the differential cross sections dσ/dM γγ , dσ/d cos θ * , dσ/d∆Φ γγ and dσ/dp T γγ , which are also considered in Ref. [4], and we use the same kinematical bins as used for the corresponding experimental data in Tables 2-5 of Ref. [4].
The LO, NLO and NNLO results (including their scale dependence) for dσ/dM γγ are presented in the main panel of Fig. 11-left. In the lower panel of Fig. 11-left we present the NNLO K factor and the relative scale dependence at NLO.
The NNLO K factor is defined as where dσ N N LO (µ)/dx is the scale dependent NNLO result for the differential cross section dσ/dx with respect to the kinematical variable x, and dσ N LO (µ 0 )/dx is the corresponding NLO result at central scales (µ R = µ F = µ 0 = M T γγ ). The relative scale dependence at NLO is defined as where dσ N LO (µ)/dx is the scale dependent NLO result.
The comparisons between the numerical values of the quantities in Eqs. (37) and (38) gives a direct quantitative illustration of the degree of overlap (including scale dependencies) of the NNLO and NLO results. By inspection of Fig. 11-left we see that the two bands in the lower panel do not overlap. This lack of overlap has already been observed for the case of the NNLO and NLO total cross sections and the same qualitative features will be observed for the differential cross sections dσ/d cos θ * , dσ/d∆Φ γγ and dσ/dp T γγ that we present in the following.
As we have already remarked, this implies that the 'true' perturbative uncertainty of the NNLO result for these observables is larger than the corresponding NNLO scale dependence that we explicitly compute.
A minor comment on the results in Fig. 11 regards the first bin, where 0 ≤ M γγ ≤ 20 GeV. As discussed in Sect. 2.3.3 (see, in particular the comments that accompany Eq. (23)), the differential cross section dσ/dM γγ in Fig. 11 has the lower kinematical limit M γγ ∼ > 10 GeV. This implies that we do not actually compute dσ/dM γγ with M γγ → 0. It also implies that the rapid decrease of dσ/dM γγ from the second to the first bin is partly an artifact of the vanishing of dσ/dM γγ in the region where M γγ ∼ < 10 GeV (this region covers about one half of the first bin). Obviously the same kinematical artifact equally affects the experimental data in the first bin (the artifact has no effect on the comparison between data and theory and on the quantities in the lower panel of Fig. 11-left).
The main features of dσ/dM γγ and of the corresponding LO and NLO results are discussed in Sects. 2.3.2 and 2.3.3 (see, in particular, Fig. 3 and related comments). That discussion also includes some comments on our expectations about the NNLO results. As noticed below, those expectations are confirmed by the NNLO results in Fig. 11. Figure 11: The differential cross section dσ/dM γγ with the same photon kinematical cuts as in Table 2. The main panel in the left-hand side shows the LO (black dotted), NLO (red dashed) and NNLO (blue solid) results, with the corresponding scale dependence. The lower subpanel presents the NNLO K factor (including its scale dependence) and the relative scale dependence at NLO. The panel in the right-hand side shows the NNLO result at central scales and its decomposition in the contributions of different initial-state partonic channels: qq (green dot-dashed), qg (black dotted) and gg (blue dashed). The separate contribution of the box gg → γγ squared amplitude is also shown (light-blue dot-dot-dashed).
The presence of the (unphysical) LO threshold at M γγ = 50 GeV is responsible for the shape of dσ/dM γγ in Fig. 11. The two invariant-mass bins (40 GeV ≤ M γγ ≤ 50 GeV and 50 GeV ≤ M γγ ≤ 60 GeV) that are closest to the LO threshold have a relatively-large size, which does not offer enough resolution to examine the detailed shape of dσ/dM γγ at M γγ ∼ 50 GeV. Therefore, we simply comment on the high-mass (M γγ ∼ > 50 GeV) and low-mass (M γγ ∼ < 50 GeV) regions (comments on the region that is very close to M γγ = 50 GeV are postponed to Sect. 3.5).
In the high-mass region, the NLO and NNLO corrections to dσ/dM γγ have a size that is similar to that of the corresponding results for the total cross section. The value of K N N LO (M γγ ) (including its scale dependence) for dσ/dM γγ is very similar to the corresponding NNLO K factor for the total cross section. In particular, K N N LO (M γγ ) is remarkably independent of M γγ : we have K N N LO (M γγ ) ≃ 1.4 for 50 GeV ∼ < M γγ ∼ < 350 GeV; at higher values of M γγ , K N N LO (M γγ ) slightly decreases and K N N LO (M γγ ) ≃ 1.3 in the highest-mass bin (500 GeV ≤ M γγ ≤ 800 GeV) of Fig. 11. The NNLO contributions of the various initial-state partonic channels to dσ/dM γγ are shown in Fig. 11-right. The gg channel gives a little NNLO contribution to dσ/dM γγ (analogously to the case of the total cross section). In particular, the results in the high-mass region of Fig. 11 explicitly show that the total contribution of the gg channel is partly smaller than the sole component of the box contribution gg → γγ. In the region where 50 GeV ∼ < M γγ ∼ < 150 GeV (which gives the bulk of the total cross section), the qq and qg channels give comparable contributions to dσ/dM γγ (the qq contribution is slightly larger). At larger values of M γγ , the relative contribution of the qq channel increases.
In the low-mass region (M γγ ∼ < 50 GeV) the NNLO (NLO) result represents, at the formal level, an effective NLO (LO) perturbative prediction. In view of that and of the discussion in Sect. 2.3.3 (see, in particular, Eqs. (28) and (29) and accompanying comments), the NNLO corrections are expected to be large. Indeed, the NNLO corrections sizeably increase by decreasing M γγ . The NNLO K factor has the value K N N LO (M γγ ) ∼ 3.8 in the second bin (20 GeV ≤ M γγ ≤ 30 GeV), and it increases to K N N LO (M γγ ) ∼ 4.4 in the first bin. The scale dependence of the NNLO result also increases by decreasing M γγ , and it reaches the size of about ±20% in the first bin. As shown in Fig. 11-right, the contribution of the gg partonic channel to the NNLO result remains small also in the low-mass region (note that the box contribution gg → γγ is absent in this region). The major contribution to the NNLO result is due to the qg channel, because both the qg and qq channels already contribute at the lowest order in the low-mass region and the PDF luminosity L qg is larger than L qq . Figure 12: The differential cross section dσ/d cos θ * . The results are analogous to those in Fig. 11.
The NNLO results for dσ/d cos θ * (and the corresponding results at LO and NLO) are presented in Fig. 12. The main features of the shape of dσ/d cos θ * are discussed in detail in Sects. 2.3.2 and 2.3.3. By direct inspection of the NLO (and LO) results in Fig. 4 and Fig. 12-left, we note that the shape of dσ/d cos θ * in Fig. 12-left is slightly different in the central region (| cos θ * | ∼ < 0.5), since it tends to be more sharpened as cos θ * → 0. This shape distortion between the results in Fig. 4 and Fig. 12 is basically produced by the exclusion of the photon rapidity region 1.37 < |y γ | < 1.52 in the case of Fig. 12.
In the central region (| cos θ * | ∼ < 0.5) we have K N N LO (cos θ * ) ≃ 1.4 . This value of the NNLO K factor and its scale dependence are remarkably similar to those of the corresponding K factors for both the differential cross section dσ/dM γγ in the high-mass region and the total cross section. In this cos θ * region the qq and qg channels give comparable NNLO contribution to dσ/d cos θ * (the NNLO contribution of the gg channel is small, at the level of about 10%), as shown in Fig. 12-right.
All these features are perfectly consistent with the fact (as discussed in detail in Sect. 2.3.2; see, in particular, the second paragraph below Eq. (18)) that the central cos θ * region is kinematically strongly correlated (through the photon p T cuts) to the M γγ region that gives the bulk of the total cross section.
At large values of | cos θ * | (| cos θ * | ∼ > 0.5), the size of the NNLO corrections to dσ/d cos θ * increases by increasing | cos θ * | (see Fig. 12-left). In the second bin (0.84 < | cos θ * | < 0.92) we have K N N LO (cos θ * ) ≃ 1.7, and in the first bin (0.92 < | cos θ * | < 1) we have K N N LO (cos θ * ) ≃ 2.4 (the NNLO scale dependence also increases in the first bin). As discussed in Sect. 2.3.2, the large | cos θ * | region is kinematically more sensitive to high values of M γγ (values of M γγ higher than those that mostly contribute to the total cross section), and it is also relatively more sensitive to the angular distribution of the underlying partonic hard-scattering processes. In particular, this region can receive enhanced NNLO corrections from the qq and qq initiated processes in Eqs. (32) and (33). The increasing value of K N N LO (cos θ * ) at large | cos θ * | is consistent with our comments that accompany Eqs. (30)- (33) in Sect. 2.3.3 . As shown in Fig. 12-right, the relative NNLO contribution to dσ/d cos θ * of the qq channel increases in the region where | cos θ * | ∼ > 0.5. This increasing behaviour of the qq channel is consistent with that observed at high values of M γγ (see Fig. 11-right).
The NNLO (and NLO) results for the differential cross sections dσ/d∆Φ γγ and dσ/dp T γγ are presented in Fig. 13. These differential cross sections have some similarities and some differences.
The similarities certainly regard the Sudakov sensitive region, namely the region that is close to the exclusive boundary of the phase space (either ∆Φ γγ ∼ π or p T γγ ∼ 0). In this region the LO result (which is not shown in Fig. 13) is non-vanishing only in the most extreme bin (the bin that includes either ∆Φ γγ = π or p T γγ = 0). As already recalled in Sect. 2.3.3, fixed-order QCD computations of dσ/dp T γγ at small values of p T γγ are affected by large logarithmic contributions (powers of ln(M γγ /p T γγ ) ∼ ln(M LO /p T γγ )) that invalidate the physical predictivity of the fixedorder result. Since p T γγ → 0 implies ∆Φ γγ → π, ensuing large logarithms (ln(π − ∆Φ γγ )) appears in the fixed-order computation of dσ/d∆Φ γγ . These large logarithmic corrections produce a rapid change of the shape of the differential cross section order-by-order in QCD perturbation theory. Such a variation of shape is clearly visible in the results of Fig. 13-left for ∆Φ γγ ∼ > 2.8 and Fig. 13right for p T γγ ∼ < 15 GeV (see, in particular, the shape of K N N LO in the lower subpanels of Fig. 13). As a consequence of this shape variation, the NLO and NNLO results tend to overlap (leading to K N N LO ∼ 1) in a tiny region. This overlap and the much reduced scale dependence of the NNLO result in this region should not be regarded as a signal of perturbative convergence: in contrast, they are just a consequence and an artifact of the order-by-order perturbative instability of the shape of the differential cross section (see Ref. [98] and, in particular, Sect. 3 therein for a completely related discussion in the context of Z production). In the Sudakov sensitive region, reliable QCD predictions requires all-order resummation of the large logarithmic contributions [42,43].
The NNLO contribution of the various initial-state partonic channels to dσ/dp T γγ in the region where p T γγ < 40 GeV is presented in Fig. 14. Note that the partial contribution of the box gg → γγ amplitude is non-vanishing only in the lowest-p T γγ bin (p T γγ < 2 GeV). In this first bin the qq and gg channels give comparable and positive contributions to the NNLO differential cross section, while the contribution of the qg channel is large and negative (it is outside the range of the vertical scale in Fig. 14). If p T γγ > 2 GeV the gg channel always gives a minor contribution: in particular, Figure 13: The differential cross sections dσ/d∆Φ γγ (left) and dσ/dp T γγ (right). The NLO and NNLO results are analogous to those in Fig. 11. the contribution is negative in the entire p T γγ region (excluding the lowest-p T γγ bin) of Fig. 14, and it becomes positive at larger values of p T γγ . At very low values of p T γγ , the qq channel gives the largest NNLO contribution (this channel is responsible for the dominant logarithmic enhancement at small values of p T γγ ). By increasing the value of p T γγ , the qg channel tends to give the largest NNLO contribution because of its larger PDF luminosity. The contribution of the qg channel remains the largest also in the region where p T γγ > 40 GeV.
We do not explicitly report the decomposition of dσ N N LO /d∆Φ γγ in initial-state partonic channels, since it is qualitatively similar to the decomposition in Fig. 13 through the correspondence (low p T γγ ) ↔ (large ∆Φ γγ ) and (large p T γγ ) ↔ (small ∆Φ γγ ). In particular, the qq channel gives the largest NNLO contribution at ∆Φ γγ ∼ > 3, whereas the qg channel gives the largest contribution at moderate and small values of ∆Φ γγ . The contribution of the gg channel is always small and, in particular, it is negative at ∆Φ γγ ∼ > 2.7 .
Outside the Sudakov sensitive region the LO result vanishes and, therefore, the NNLO (NLO) results in Fig.13 are 'effective' NLO (LO) QCD predictions. We comment on the results for dσ/d∆Φ γγ and dσ/dp T γγ in turn.
In the region from moderate to small values of ∆Φ γγ (Fig. 13-left) the NNLO corrections monotonically increase by decreasing ∆Φ γγ . We have K N N LO (∆Φ γγ ) ∼ 2 at ∆Φ γγ ∼ 2.5, and K N N LO (∆Φ γγ ) ≃ 3.4 in the lowest-∆Φ γγ bin (0 < ∆Φ γγ < 0.5). The large size of the NNLO corrections at small ∆Φ γγ was pointed out in Ref. [99]. The scale dependence of the NNLO result throughout this region is at the level of about ±15%, and it is basically unchanged with respect to the scale dependence of the NLO result.
In the region from moderate to large values of p T γγ the size of the NNLO corrections does not have a monotonic dependence on p T γγ (Fig. 13-right). We have K N N LO (p T γγ ) ∼ 1.7 at Figure 14: The NNLO cross section dσ/dp T γγ at central scales in the low-p T γγ region of Fig. 13right. The NNLO result is decomposed in the contributions of different initial-state partonic channels: qq (green dot-dashed), qg (black dotted) and gg (blue dashed). The box gg → γγ squared amplitude only contributes in the first (lowest-p T γγ ) bin.
p T γγ ∼ 20 GeV and K N N LO (p T γγ ) ≃ 1.8 in the highest-p T γγ bin (250 GeV < p T γγ < 500 GeV); at these values of p T γγ the NNLO scale dependence is at the level of about ±12% (the corresponding NLO scale dependence is about ±18%). In the intermediate p T γγ region, the size of the NNLO corrections has a maximal value at p T γγ ∼ 50 GeV. In the bin with 50 GeV < p T γγ < 55 GeV we have K N N LO (p T γγ ) ≃ 2.4 with a scale dependence of about ±15% (it is basically the same scale dependence as at NLO, since in the region where 50 GeV ∼ < p T γγ ∼ < 100 GeV the NLO scale dependence tends to be minimal).
The NNLO results for dσ/d∆Φ γγ and dσ/dp T γγ deserve some overall comments.
We first note that the size of the NNLO corrections to dσ/d∆Φ γγ and dσ/dp T γγ is typically larger than that of the NNLO corrections to the total cross section, to dσ/dM γγ at high masses and to dσ/d cos θ * . This feature is consistent with the fact that we are dealing with 'effective' NLO (rather than NNLO) QCD predictions in the case of the ∆Φ γγ and p T γγ differential cross sections.
Then, we also observe a hierarchy of the size of the NNLO corrections for various differential cross sections that are computed at 'effective' NLO: K N N LO (M γγ ) at low M γγ is larger than K N N LO (∆Φ γγ ) at small ∆Φ γγ , which is in turn larger than K N N LO (p T γγ ) at moderate and large values of p T γγ . This hierarchy is in agreement with the expectations of our discussion in Sect. 2.3.3. In particular, it is qualitatively consistent with the effect of the NNLO processes in Eqs. (28) and (29): they are dynamically enhanced at small values of M γγ (see the accompanying comments to Eqs. (28) and (29)) and, because of kinematics, they are still enhanced but with a decreasing relevance at small ∆Φ γγ and, in turn, at large p T γγ .
Finally, we note that the values of K N N LO for dσ/d∆Φ γγ and dσ/dp T γγ within smooth isolation are larger than the corresponding differences between standard and smooth isolation results at NLO (see Sect. 2.3.3). This fact and the discussion throughout the paper is a strong indication of the presence of sizeable NNLO radiative corrections to diphoton production also within the standard cone isolation criterion.
We present some additional comments on dσ/dp T γγ in the region of relatively-large values of p T γγ . As we have noticed in Sect. 2.3.3, the NLO standard isolation results in Fig. 8-right have a shoulder-type behaviour at p T γγ ∼ 50 GeV. The p T γγ shoulder is also visible in the NNLO smooth isolation results of Fig. 13-right (see also Figs. 15-left and 18-right, which have a linear scale in p T γγ ), although it is less pronounced since the results refer to E T max = 4 GeV, which is smaller than the value E T max = 10 GeV of Fig. 8-right. As previously recalled, in the context of standard isolation the p T γγ shoulder was discussed in Ref. [95]. Here we present related comments in the context of the smooth isolation criterion. Figure 15: Left panel: the NNLO cross section dσ/dp T γγ (blue dashed) in Fig. 13-right and its partition into the contributions from the two complementary regions where ∆Φ γγ > π/2 (green dotted) and ∆Φ γγ < π/2 (red solid). Right panel: the fractional contribution to dσ/dp T γγ from the region where ∆Φ γγ < π/2 at NLO (blue dashed) and NNLO (red solid).
We follow Ref. [95] and we consider the partition of dσ/dp T γγ in the contributions of the two complementary kinematical regions with relatively-large (∆Φ γγ > π/2) and relatively-small (∆Φ γγ < π/2) values of ∆Φ γγ . These contributions are denoted by dσ (l) /dp T γγ (if ∆Φ γγ > π/2) and dσ (s) /dp T γγ (if ∆Φ γγ < π/2) in the following text. The complete NNLO result for dσ/dp T γγ and the corresponding NNLO results from the two complementary regions of ∆Φ γγ are presented in Fig. 15-left at values of p T γγ around the shoulder. The component dσ (l) /dp T γγ has an approximately-constant slope at p T γγ ∼ > 30 GeV, whereas the total contribution tends to flatten out at p T γγ ∼ 50 GeV and it remains systematically higher than dσ (l) /dp T γγ at larger values of p T γγ . These differences are obviously due to the component dσ (s) /dp T γγ , which exactly vanishes at p T γγ ∼ < 34 GeV and has a maximum value at p T γγ ∼ 50 GeV.
To comment on the behaviour of the results in Fig. 15-left we consider the effect of the kinematical constraint in Eq. (24) (see also its accompanying comments), which is a consequence of the photon p T cuts. The component dσ (l) /dp T γγ is insensitive to the constraint. Moving from large to smaller values of p T γγ , the size of this component increases since the production of photons with smaller transverse momenta is kinematically (because of energy conservation) and also dynamically favoured. At ∆Φ γγ ∼ 0, the constraint in Eq. (24) leads to p T γγ ∼ > 47 GeV (see Eq. (25)). Therefore, provided p T γγ ∼ > 47 GeV, the entire kinematical region of ∆Φ γγ is kinematically allowed and dσ (s) /dp T γγ is also insensitive to the constraint. It follows that at large values of p T γγ both components (and their total contribution) show a similar increasing behaviour as p T γγ decreases. If p T γγ ∼ < 47 GeV the constraint in Eq. (24) is instead effective on dσ (s) /dp T γγ and, consequently, on the total contribution to dσ/dp T γγ . Going from p T γγ ≃ 47 GeV to smaller values of p T γγ , the kinematical constraint forbids the region of small values of ∆Φ γγ : consequently, the increasing behaviour of dσ/dp T γγ flattens out since dσ (s) /dp T γγ is increasingly suppressed. In particular, Eq. (24) implies the restriction p T γγ ∼ > 34 GeV if ∆Φ γγ < π/2 and, therefore, dσ (s) /dp T γγ has a kinematical threshold at p T γγ ≃ 34 GeV.
From our discussion on the results in Fig. 15-left, it follows that the onset of the shoulder at p T γγ ∼ 50 GeV is kinematically driven (by the presence of the photon p T cuts). However, it also has a dynamical component, since the NNLO corrections enhance this effect (see the size of K N N LO (p T γγ ) at p T γγ ∼ > 50 GeV in Fig. 13-right).
To remark the NNLO enhancement we compute the ratio between dσ (s) /dp T γγ and dσ/dp T γγ , namely, the fractional contribution to dσ/dp T γγ from the region of relatively-small values of ∆Φ γγ . In Fig. 15-right we present the results for this ratio at both NNLO and NLO (the shape of dσ (s) /dp T γγ at NLO is qualitatively similar to that of the corresponding NNLO result in Fig. 15left). At both perturbative orders this ratio is sizeable and approximately independent of p T γγ in the region just above p T γγ ∼ 50 GeV: the ratio is about 0.42 at NLO and 0.54 at NNLO. Therefore, the region where ∆Φ γγ < π/2 contributes to roughly half of the complete differential cross section dσ/dp T γγ at p T γγ ∼ > 50 GeV, and this fact has a kinematical origin. The increase of about 30% of the ratio from the NLO to the NNLO result has instead a dynamical origin (incidentally, this increase implies that the NNLO K factor for dσ (s) /dp T γγ is larger than that for dσ/dp T γγ in Fig.13-right). It is due to the NNLO processes in Eqs. (28) and (29) that are dynamically enhanced at small values of M γγ and, consequently (see Eq. (23) and accompanying comments), partly enhanced also at very small values of ∆Φ γγ . There are dynamical similarities between these processes and the NLO fragmentation processes in Eqs. (26) and (27), which, according to Ref. [95], originate the p T γγ shoulder in the NLO standard isolation results. Therefore, our discussion on the p T γγ shoulder in the NNLO smooth isolation results is consistent with the observations in Ref. [95].
The partition of the phase space into the two regions with ∆Φ γγ > π/2 and ∆Φ γγ < π/2 was applied by the D0 Collaboration [6] to diphoton production data in proton-antiproton collisions at the Tevatron. The results in Ref. [6] show that the inclusion of the NNLO corrections considerably improves the description of the data in both phase space regions.

Dependence on isolation parameters
We add some illustrative results and comments on differential cross sections at NNLO and their dependence on isolation parameters.
We use dσ/d cos θ * and dσ/d∆Φ γγ as representative cross sections of typical features of the NNLO results. We consider the corresponding reference NNLO results of this section (Figs. 12 and 13), which are obtained with E T max = 4 GeV and the power n = 1 in the isolation function χ(r; R), and we perform variations of the isolation parameters E T max and n. Specifically, we use either n = 0.5 or n = 2 by keeping E T max = 4 GeV fixed, and we use n = 1 by increasing E T max to the value E T max = 10 GeV. The ensuing NNLO results are presented in Fig. 16. In the lower subpanels of Fig. 16 we present the ratio of the results with different isolation parameters with respect to the reference NNLO result (n = 1, E T max = 4 GeV). All these results are obtained with central values of the scales (µ R = µ F = µ 0 = M T γγ ). The scale dependence of the reference NNLO result is also reported in Fig. 16 for comparison with isolation parameters effects. We comment on the results in Fig. 16 by excluding those in the region where ∆Φ γγ ∼ > 2.8 (comments on this region are postponed). As a main overall comment we observe that the NNLO results have a mild dependence on the isolation parameters. Indeed, the effects of the variations of the isolation parameters are not larger than the scale dependence at fixed isolation parameters. Moreover, the qualitative dependence on both E T max and n is in agreement with physical expectation (see Eqs. (6) and (8)) throughout the entire kinematical regions of cos θ * and ∆Φ γγ : the value of the differential cross sections decreases by requiring the photons to be more isolated. We also recall (see Sect. 2.3.3) that in the region where ∆Φ γγ ∼ < π − R ≃ 2.7 the NLO cross section dσ/d∆Φ γγ (within smooth isolation) does not depend on the isolation parameters. Therefore, the E T max and n dependence observed in the NNLO results for dσ/d cos θ * (dσ/d∆Φ γγ ) is 'effectively' an NLO (an LO) QCD effect.
We note that the n dependence (at fixed E T max = 4 GeV) is practically the same throughout the entire kinematical region of either cos θ * or ∆Φ γγ , so that it does not produce any shape variations. The overall size of the n dependence of the differential cross sections is basically equal to that of the NNLO total cross section (see Sect. 3.1).
The E T max dependence in Fig. 16 is substantially larger than the n dependence, because E T max = 10 GeV is sizeably different from E T max = 4 GeV. Nonetheless, also the E T max dependence at NNLO is quite moderate (as already noticed in Sect. 2.3.1 for the case of the NLO total cross section). In the case of dσ N N LO /d cos θ * , the variation of E T max does not produce any significant shape variation: the differential cross section uniformly increases by approximately 10% in going from E T max = 4 GeV to E T max = 10 GeV. Consequently, an equal increase ap-plies to the NNLO total cross section. We note that the size of the E T max dependence is slightly larger at NNLO than at NLO (see Sect. 2.3.1); an analogous increased sensitivity to the value of n (at fixed E T max ) has been already noticed in Sect. 3.1 . In the case of dσ N N LO /d∆Φ γγ , the variation of E T max also produces a visible shape dependence. In going from E T max = 4 GeV to E T max = 10 GeV, dσ N N LO /d∆Φ γγ increases by approximately 8% at ∆Φ γγ ∼ 2, and by approximately 16% at ∆Φ γγ ∼ 0.5 . This shape dependent effect is consistent with our previous discussions in Sect. 3.2 and in Sect. 2.3.3. At small values of ∆Φ γγ the NNLO cross section dσ N N LO /d∆Φ γγ receives substantial contribution from configurations in which the photons are accompanied by partonic (hadronic) transverse energy inside the isolation cones (see Eqs. (28) and (29)): by increasing E T max these configurations are less suppressed by the isolation requirement and the cross section increases (a qualitatively similar E T max dependence at small ∆Φ γγ is observed in the NLO standard isolation results of Fig. 7).
On the basis of our discussion about the results in Fig. 16 and of our previous discussions about similarities between various NNLO differential cross sections, we can argue about the E T max dependence of dσ/dM γγ and dσ/dp T γγ at NNLO. In the high-M γγ region the E T max dependence of dσ/dM γγ is similar to that of dσ/d cos θ * . The E T max sensitivity of dσ/dM γγ at small values of M γγ is partly enhanced with respect to that of dσ/d∆Φ γγ at small ∆Φ γγ . The E T max sensitivity of dσ/dp T γγ at moderate and large values of p T γγ is partly reduced with respect to that of dσ/d∆Φ γγ at small ∆Φ γγ .
We come to comment about the results in the large-∆Φ γγ region (∆Φ γγ ∼ > 2.8) of Fig. 16. This is definitely part of the Sudakov sensitive region. We point out that here the dependence on the isolation parameters (both E T max and n) is amplified with respect to the rest of the kinematical regions. We are dealing with diphoton production close to the exclusive boundary of the phase space, where the photon pair is accompanied by radiation of little transverse energy, also outside the photon isolation cones. In such a configuration the effects of variations of the isolation requirements are relatively enhanced at fixed order: therefore, an enhanced sensitivity to isolation parameters can be expected in the results of Fig. 16-right. As recalled in Sect. 3.2, the NNLO result (and, consequently, the ensuing isolation parameter dependence) is not physically reliable in this Sudakov sensitive region. Nonetheless, the observed isolation parameter dependence at NNLO indicates that photon isolation effects have to be carefully considered (and examined) in the context of refined resummed calculations in the Sudakov sensitive region.

Comparison of NNLO results and data
In Ref. [4] the ATLAS data on diphoton production at the centre-of-mass energy √ s = 7 TeV were compared with the NNLO QCD results. As stated in the conclusion of Ref. [4], the NNLO results are able to match the data very closely within the uncertainties, except in limited regions. In the following we present some comments on the comparison between NNLO results and data.
As a preliminary comment, we note that the NNLO results in Ref. [4] and those in this paper are not exactly equal. The differences have various origins. Part of the differences are due to the setup of the NNLO calculations. The results of the NNLO calculation presented in Ref. [4] use the PDFs of the MSTW2008 set, the central scale µ 0 = M γγ and the scale dependence that is obtained by considering the scale configurations with µ R = µ F = M γγ /2 and with µ R = µ F = 2M γγ . Here we use different PDFs, a different central scale and independent scale variations of µ R and µ F . The other differences are due to the fact that the NNLO results of Ref. [4] were obtained by using the first version of the numerical program 2γNNLO, which had an implementation error that was subsequently corrected [35]. All these effects produce differences between the NNLO results used in Ref. [4] and those reported here. Nonetheless the differences are relatively small (in particular, they are smaller than the scale dependence at NNLO) and they do not have a major effect on the comparison between NNLO results and data.
The ATLAS Collaboration performs a quantitative estimate of underlying event, pile-up and hadronization effects. These effects, which are not included in fixed-order QCD calculations, are taken into account [4] by applying a correction factor to the fixed-order QCD results. The typical size of the correction factor is around 0.95 [4]. The bin-by-bin correction factors for various differential cross sections are available in the database of the Durham HepData Project [100]. We apply these bin-by-bin correction factors to our NNLO calculations (for simplicity, we simply rescale the NNLO central values by the bin-by-bin corrections, without including the uncertainties of the correction factors themselves), and we present the ensuing results in Figs. 17 and 18. A technical comment about Fig. 17 regards the differential cross section dσ/d cos θ * . Since the two photons of the diphoton pair are identical particles, dσ/d cos θ * is symmetric with respect to cos θ * ↔ − cos θ * . In our QCD results throughout this paper we have always computed the differential cross section with respect to | cos θ * | and then the results are presented in symmetrized form. In the ATLAS data, cos θ * refers to the cosine of the polar angle θ * of the harder photon (the photon with momentum p hard T γ ) and, consequently (because of experimental effects in the measurement, including statistical fluctuations), dσ/d cos θ * is not exactly symmetric under cos θ * ↔ − cos θ * (although it is symmetric within the experimental errors). The NNLO result in Fig. 17-right is not exactly symmetric because of the non-symmetric correction factors that have been applied.
A main general comment about the comparison between data and NNLO results regards the photon isolation criteria. The ATLAS data use standard isolation (though the experimental details of the actual isolation procedure are quite involved and differ from a plain implementation of standard isolation), whereas the NNLO results use smooth isolation with the same values of the isolation parameters R and E T max . In our opinion the comparison is meaningful despite the differences between the isolation criteria. This follows from several observations that we list. Owing to the physical constraint in Eq. (5), NNLO results for smooth isolation give a lower bound on the NNLO results for standard isolation. Smooth and standard isolation results are quantitatively quite similar starting from computations at 'effective' NLO accuracy (see Sects. 2.3.1 and 2.3.3), and such similarity is expected to remain valid at NNLO (see, e.g., the mild dependence on isolation parameters that we noticed in Sect. 3.3). Calculations within both isolation criteria are affected by sizeable higher-order corrections, and at 'effective' NLO accuracy these corrections tend to be larger than the differences between the NLO results for smooth and standard isolation (see Sects. 2.3.1-2.3.3). Therefore, neglecting NNLO effects can have more impact than using different isolation criteria. Reference [4] compares the ATLAS data to both NNLO results for smooth isolation and lower-order results for standard isolation: the differences between these fixed-order results confirm that NNLO corrections are sizeable and relevant. Obviously all these observations about smooth and standard isolation results are valid within the corresponding perturbative uncertainties.
We present more specific comments about total and differential cross sections.
The experimental value of the total cross section is [4] σ AT LAS = 44.0 +3.2 −4.2 pb. We compute the corresponding NNLO result from the bin-by-bin corrected differential cross section dσ/d cos θ * (Fig. 17-right), and we obtain σ N N LO AT LAS = 37.2 +3.2 −3.3 (scale) pb. The measured value and the NNLO result are consistent within the corresponding experimental and scale dependence uncertainties. Moreover, as remarked in Sect. 3.1, the NNLO scale dependence cannot be regarded as a consistent estimate of the perturbative uncertainty due to uncalculated higher-order terms. The 'true' perturbative uncertainty of σ N N LO is larger than the scale dependent effect that we have computed. A better uncertainty estimate can be obtained, for instance, by considering enlarged scale variations, such that the ensuing NNLO and NLO results have some degree of overlap. An alternative, similar and simpler procedure consists in comparing the NNLO and NLO results at central values of the scales (incidentally, we have checked that NLO results for smooth and standard isolation are quantitatively very similar) and using half of the difference between them to assign the perturbative uncertainty of the NNLO result. This procedure leads (see Table 2) to an NNLO uncertainty of about ±14%, which amounts to almost double the scale uncertainty of the NNLO result for the total cross section.
The ATLAS data [4] on differential cross sections are reported in Figs. 17 and 18 together with the corresponding NNLO results. In the lower subpanels of these figures we present the ratio between the data and the NNLO results at central scales, and the scale dependence of the NNLO results. At large values of ∆Φ γγ (e.g., ∆Φ γγ ∼ > 2.8) and small values of p T γγ (e.g., p T γγ ∼ < 20 GeV) the shape of the data (Fig. 18) is definitely different from that of the NNLO results. This is expected since, as discussed in Sect. 3.2, the shape of the NNLO results is not physically reliable in these scale dependence) for dσ/d∆Φ γγ (left) and dσ/dp T γγ (right). The NNLO results are corrected for hadronization and underlying event effects (see text).
Sudakov sensitive regions. All-order resummed calculations [42,43,8] have to be used here. Allorder resummation effects are implemented also in parton shower event generators, which can lead to a consistent description of the data [4,8].
The NNLO results in Fig. 18 (outside the Sudakov sensitive regions) and those in the low-mass region (M γγ < 50 GeV) of Fig. 17-left are perturbative results at 'effective' NLO accuracy, while those in the high-mass region and in the entire region of cos θ * (Fig. 17) are results at 'effective' NNLO accuracy. In 'effective' NNLO regions there is a good agreement between data and NNLO results within the corresponding experimental and scale uncertainties. In 'effective' NLO regions the data tend to be systematically higher than the NNLO results at central values of the scales. Note that the shape of the differences between data and NNLO results tends to qualitatively follow the shape of the NNLO K factors (see the lower subpanels in Figs. 17 and 18 and those in Figs. 11, 12 and 13). This feature is consistent with our previous remarks on the fact (see the differences between NLO and NNLO results in Figs. 11, 12 and 13) that the computed scale dependence at NNLO does underestimate the true perturbative uncertainties of the NNLO results. If we proceed to assign a perturbative uncertainty on the basis of the differences between NNLO and NLO results for differential cross sections (analogously to the procedure that we have previously mentioned for the total cross section in this subsection), the NNLO uncertainty is almost doubled with respect to the NNLO scale dependence in most of the regions of Figs. 17 and 18 and it is further increased in limited extreme regions (such as at very small values of M γγ , at very large values of | cos θ * | and at very low values of ∆Φ γγ ). Such an NNLO uncertainty makes the ATLAS data consistent with the NNLO results in both 'effective' NNLO and 'effective' NLO regions.
The tendency of the data to be systematically higher than the NNLO results at central values of the scales is not inconsistent with the expectation (see Eq. (5)) that the NNLO results with smooth isolation should be a lower bound on the corresponding results for standard isolation. One can try to reduce the differences between the smooth and standard isolation criteria by decreasing the isolation effects within smooth isolation. This can be done by using different smooth isolation parameters such as, for instance, a smaller value of the power n of the isolation function χ(r; R) or a larger value of E T max . As shown by the results in Sect. 3.3, such a procedure can reduce the systematic differences between ATLAS data and NNLO smooth isolation results. Nonetheless, this 'tuning' procedure does not affect the overall features of the comparison between data and NNLO results, because of the relatively-small dependence on the isolation parameters (see Sect. 3.3) and of the substantial NNLO theoretical uncertainties that we have previously discussed.
In this subsection we have considered the ATLAS diphoton data at √ s = 7 TeV and we have presented comments on the main features of the data/NNLO comparison. The same features and comments are equally valid for other LHC diphoton measurements [7,8] and related data/theory comparisons [7,8,37]. In particular, the inclusion of the NNLO radiative corrections greatly improves the description of the data [4,7,8] with respect to lower-order results. The effect of the NNLO corrections is from moderate to sizeable in different kinematical regions. This leads to a good or consistent (depending on kinematical regions) agreement with the data by taking into account the perturbative uncertainty of the NNLO results.

Asymmetric and symmetric photon p T cuts
Throughout the paper we have remarked that the presence of the photon p T cuts (p hard T γ ≥ p H , p sof t T γ ≥ p S ) has relevant effects on various features of diphoton production observables. All the numerical results presented so far (with the exception of those in Fig. 6) use the values p H = 25 GeV and p S = 22 GeV. In this subsection we present some results and related comments on effects that are due to variations of p H and p S . In particular, we consider symmetric p T cut configurations (p H = p S ) in addition to asymmetric configurations (p H > p S ).
We consider the reference kinematical configuration (and theoretical setup) used throughout this section, but we vary the photon p T cuts p H and p S . All the numerical results presented in this subsection refer to the central value of the scales (µ R = µ F = µ 0 = M 2 γγ + p 2 T γγ ).
We begin our quantitative study by fixing p S = 22 GeV and varying the value of p H . The dependence on p H is parametrized by ∆p T = p H −22 GeV. The numerical results for the total cross section at LO, NLO and NNLO in the range 3 GeV < ∆p T < 50 GeV are reported in Fig. 19-left. At each of these perturbative orders we note that the total cross section σ(∆p T ) monotonically decreases by increasing p H (i.e, the value of ∆p T ), in agreement with physical expectations. We also note that the size of the NNLO corrections, as given by the K factor K N N LO (∆p T ) = σ N N LO (∆p T )/σ N LO (∆p T ) (lower subpanel of Fig. 19-left), is very weakly dependent on ∆p T : it varies in the range between 1.4 and 1.5 .
This weak dependence on ∆p T of the NNLO radiative corrections has to be contrasted with a corresponding strong dependence of the NLO radiative corrections. In the results of Fig. 19left the NLO K factor K N LO (∆p T ) = σ N LO (∆p T )/σ LO (∆p T ) is much larger than K N N LO (∆p T ) and it sizeably increases by increasing ∆p T . We have K N LO ∼ 3 at ∆p T ∼ 3 GeV, K N LO ∼ 4 at ∆p T ∼ 10 GeV, and K N LO ∼ 5.4 at ∆p T ∼ 50 GeV. At the LO, values of p sof t T γ smaller than p H are Figure 19: Dependence on the photon p T cuts p H and p S (p hard T γ ≥ p H , p sof t T γ ≥ p S ) of the total cross section for diphoton production (the other kinematical cuts are the same as in Table 2). The NLO (red dashed) and NNLO (blue solid) results are separately presented in the regions of asymmetric (left panel) and nearly-symmetric (right panel) p T cuts. The LO result (black dotted) is presented only in the case of asymmetric p T cuts (left panel). The lower subpanels present the corresponding NLO (red dashed) and NNLO (blue solid) K factors.
kinematically forbidden, and they do not contribute to σ LO . At the NLO the kinematical subregion with p S < p sof t T γ < p H is allowed: by increasing its extension (i.e. by increasing ∆p T = p H −p S ) this subregion produces increasingly large NLO corrections. The p sof t T γ region where p S < p sof t T γ < p H is kinematically allowed at both NLO and NNLO, and the NNLO K factor turns out to be weakly dependent on the asymmetry of the photon p T cuts.
We now move to consider the total cross section in configurations with symmetric (p H = p S ) or nearly-symmetric (p H ∼ p S ) p T cuts. As a reference symmetric configuration we consider the case with p H = p S = 22 GeV. Nearly-symmetric configurations can be obtained by either increasing p H at fixed p S = 22 GeV (we define ∆p T = p H − 22 GeV > 0) or decreasing p S at fixed p H = 22 GeV (we define ∆p T = p S − 22 GeV < 0). The NLO and NNLO results for the total cross section as a function of ∆p T are reported in Fig. 19-right. More precisely, σ(∆p T ) is computed at nine values of ∆p T , namely ∆p T (GeV) = {−5, −3, −1, −0.5, 0, +0.5, +1, +3, +5}, and the results are reported as data points in Fig. 19-right (the continuous lines in Fig. 19-right are just a graphical interpolation between the data points).
According to physical expectations, the total cross section σ(∆p T ) should be a monotonically decreasing function of ∆p T . Indeed, the amount of physical diphoton events that contribute to σ(∆p T ) decreases (increases) by increasing p H (decreasing p S ). The results in Fig. 19-right (especially the NLO results) do not show such a physical behaviour. In particular, σ N LO (∆p T ) has a local maximum at ∆p T ∼ 1 GeV and a local minimum at ∆p T = 0. Considering the NLO results with ∆p T ≥ 0 in Fig. 19-right, the unphysical behaviour could be ascribed to either the local maximum (a cross section value that is unphysically large at some finite, though small, value of ∆p T ) or the local minimum (a cross section value that is unphysically small at ∆p T = 0). However, the NLO cross section does not show an evident unphysical (non-monotonic) behaviour at finite and negative values of ∆p T . Therefore, we can conclude that the pathological behaviour of the NLO results in Fig. 19-right is located in the region of symmetric (or nearly-symmetric) p T cuts (|∆p T | → 0).
The unphysical behaviour of the fixed-order results for the diphoton total cross section in the presence of symmetric p T cuts is expected, since a similar behaviour was first observed and discussed in Ref. [101] in the context of dijet photoproduction. Specifically, in the case of nearlysymmetric cuts the NLO cross section behaves as This implies that the local minimum of σ N LO at ∆p T = 0 is unphysical. This also implies that σ N LO in Fig. 19-right has a cusp at ∆p T = 0, with a ∆p T -slope that diverges to +∞ (−∞) if ∆p T tends to vanish from positive (negative) values of ∆p T . The NLO numerical results in Fig. 19-right are consistent with these features, although the cusp behaviour is not clearly visible since it is located in a very narrow region of small values of |∆p T |. The double-logarithmic behaviour in the right-hand side of Eq. (39) exactly follows from the same reasoning as used in Ref. [102] in the context of dijet production. We note that the behaviour in the right-hand side of Eq. (39) is different from the single-logarithmic behaviour in the corresponding expression of Eq. (2.9) of Ref. [101]. The NLO single-logarithmic contribution of Ref. [101] is produced by hard-collinear radiation from the initial state (see Eq. (2.8) therein). The dominant double-logarithmic term in Eq. (39), which originates from the NLO process in Eq. (36), is instead due to radiation that is both soft and collinear to the direction of the initial-state colliding partons (related comments on this effect are postponed to our discussion of the results in Fig. 20) The unphysical behaviour of the total cross section at p H ≃ p S persists at each subsequent perturbative orders and a physical (smooth and monotonic) dependence of σ(∆p T ) on ∆p T can be recovered only by a proper all-order resummation of soft-gluon effects [102]. Such a resummation is beyond the scope of the present paper and we limit ourselves to comment on the quantitative reliability of the fixed-order diphoton results that we have presented.
By direct inspection of the NLO and NNLO results in Fig. 19-right, we tend to conclude that the onset of the unphysical fixed-order behaviour due to nearly-symmetric p T cuts occurs in a region of small values of ∆p T , such as |∆p T | ∼ < 2 GeV. Therefore, in the case of asymmetric cuts with p H − p S ∼ > 3 GeV, we argue that the unphysical behaviour has little effect on the total cross section.
At smaller values of p H − p S we can try to quantify the effect of the unphysical behaviour. For instance, at each fixed order we can assume that a 'tentative' physical value of the cross section is in the range between the values of σ(∆p T ∼ −2 GeV) and σ(∆p T ∼ +2 GeV) at the corresponding order. Then, we can use the size of this range and the difference with respect to the 'unphysical' computed value of σ(∆p T ) at small |∆p T | to assign a systematic theoretical uncertainty to this computed value. From the results in Fig. 19-right, this procedure leads to effects of about 10%-15% at NLO and of several percent at NNLO. These quantitative effects can be regarded as a rough estimate of the uncertainty due to the unphysical soft-gluon effects at fixed orders. We note that such uncertainty quantitatively decreases by increasing the perturbative order. This decrease is expected since the unphysical behaviour is located in an increasingly smaller region of ∆p T by increasing the perturbative order (the effect is visible in the NLO and NNLO results of Fig. 19-right). We also note that such uncertainty is quantitatively similar to the scale dependence of the cross section at the corresponding order (see Table 2). More importantly for 'practical' purposes, such uncertainty is definitely smaller than the typical size (about 40%) of the NNLO corrections (the value of the NNLO K factor in the lower subpanel of Fig. 19-right increases very slightly, from 1.4 to 1.5, toward the region where ∆p T ∼ 0). This implies that even in the case of nearly-symmetric p T cuts the bulk of the NNLO corrections to the total cross section is due to hard-parton radiation rather than to unphysical soft-gluon effects. This also implies (as we have discussed throughout this section and, in particular, in Sect. 3.4) that the perturbative uncertainty of the NNLO result for the total cross section is dominated by the effect of the large NNLO radiative corrections.
From our discussion of the results in Fig. 19, we can make an overall comment about the total cross section at NNLO: the size of the NNLO corrections and of the NNLO theoretical uncertainties weakly depends on the amount ∆p T of the asymmetry of the photon p T cuts. Some validation of this overall comment can be found in features of data/theory comparisons. Indeed, LHC measurements of diphoton total cross sections have been performed in configurations with different values of ∆p T (∆p T = 3 GeV [4], 10 GeV [8] and 15 GeV [7]), and in all these configurations the comparison between data and NNLO results (see Sect. 3.4 and Refs. [4,7,8,37]) shows a similar degree of consistency.
We further study the effects produced by symmetric p T cuts by considering some differential cross sections. Analogously to the case of the total cross section, we consider the symmetric p T cut configuration with p H = p S = 22 GeV and we compute the M γγ -differential cross section dσ/dM γγ and the inclusive transverse-momentum spectra dσ/dp hard T γ and dσ/dp sof t T γ of the harder and softer photon. The numerical results at NLO and NNLO are presented in Fig. 20. Incidentally, we note that the main features of the results in Figs. 19 and 20 do not depend on the specific value p H = p S = 22 GeV of the symmetric p T cuts (we have explicitly checked this by considering values of p H = p S in the range between 20 GeV and 30 GeV).
In Fig. 20-left we report the results for dσ/dM γγ in the M γγ region that is close to the LO threshold at M γγ = M LO dir . For comparison, we present the results for two different configurations with asymmetric (p H = 25 GeV, p S = 22 GeV) and symmetric (p H = p S = 22 GeV) p T cuts. We remark that the two results are obtained by only varying the value of p H , while all the other kinematical cuts and parameters of the calculation are unchanged.
We first briefly comment on the case with asymmetric p T cuts. The NLO and NNLO results in Fig. 20-left exactly correspond to those in Figs. 11-left and 17-left, the only difference being the much smaller M γγ bin size, which is equal to 0.2 GeV. The behaviour of the LO and NLO results for dσ/dM γγ has been discussed in detail at the end of Sect. 2.3.3 (see Fig. 10). In particular, the LO result has a threshold at M γγ = M LO dir = 50 GeV and the NLO result has an upward double-side cusp at M γγ = M LO dir . The NLO cusp behaviour is due to soft-gluon radiation effects [60]. The 'unphysical' soft-gluon effects persist at the NNLO level and their dominant (at the formal level) contribution leads to a negative double-logarithmic enhancement (∝ − ln 2 (ǫ M )) of the NLO cusp behaviour in Eq.  and NNLO (solid) results for the differential cross sections dσ/dp hard T γ (blue) and dσ/dp sof t T γ (red) of the harder and softer photon in the configuration with symmetric p T cuts, p H = p S = 22 GeV.
The shape of dσ/dM γγ is quite different in the two configurations with asymmetric and symmetric p T cuts (Fig. 20-left). The LO result with symmetric p T cuts is not shown in Fig. 20-left. Its shape is exactly similar to that with asymmetric p T cuts (see Fig. 10 and Eq. (34)) since the value of p S does not matter at the LO and, in particular, dσ LO /dM γγ has its threshold at M γγ = M LO dir = 44 GeV (since p H = 22 GeV). In the case of symmetric p T cuts, the NLO and NNLO values of dσ/dM γγ are quite small below the LO threshold. Just above the LO threshold (M γγ ∼ > M LO dir ), the NLO result is (relatively) large and negative and the NNLO result is also (relatively) large but positive. This NLO behaviour is unphysical. In particular, the relatively-large and negative differential cross section dσ N LO /dM γγ at M γγ ∼ 44 GeV 'explains' the unphysical behaviour of the NLO total cross section at ∆p T = 0. In the case of symmetric (or nearly-symmetric) p T cuts dσ N LO /dM γγ becomes negative at M γγ ∼ M LO dir and, after integration over M γγ , this negative contribution is responsible for the decreasing behaviour (see Fig. 19-right) of the NLO total cross section at small values of ∆p T (0 ∼ < ∆p T ∼ < 1 GeV). We note that also dσ/dM γγ (and not only the total cross section) for symmetric p T cuts is physically expected to be larger than the corresponding differential cross section for asymmetric p T cuts. This physical expectation is not fulfilled by the NLO results of Despite their apparent shape difference, the fixed-order behaviour of dσ/dM γγ at M γγ ∼ M LO dir for asymmetric and symmetric p T cuts is produced by the same underlying mechanism [60], as we are going to discuss below. In this sense, the unphysical fixed-order behaviour of the total cross section at p H ≃ p S can also be regarded as a consequence (after integration) of unphysical soft-gluon effects for non-smooth differential distributions [60] (e.g., dσ/dM γγ and, also, dσ/dp hard T γ and dσ/dp sof t T γ as discussed below).
Analogously to the case of asymmetric p T cuts (see Eq. (35)), in the case of symmetric p T cuts we have examined the NLO shape of dσ/dM γγ at M γγ ∼ M LO dir in analytic form, and we find the dominant behaviour where b 0 and b (+) are positive constants (i.e., they do not depend on M γγ ) and the dots in the righthand side denote subdominant contributions in the limit M γγ → M LO dir . According to Eq. (40), dσ N LO /dM γγ is finite at M γγ = M LO dir . Its behaviour just above the LO threshold (M γγ > M LO dir ) is analogous to that in Eq. (35) and, in particular, the first derivative of dσ N LO /dM γγ with respect to M γγ (i.e., the slope of dσ N LO /dM γγ ) diverges to −∞. Therefore, dσ N LO /dM γγ has an upward cusp at M γγ = M LO dir . At variance with Eq. (35), the behaviour of the result in Eq. (40) is smooth just below the LO threshold (M γγ < M LO dir ) and, in particular, the slope of dσ N LO /dM γγ is finite (subdominant terms of O(ǫ M ) are neglected in the right-hand side of Eq. (40)) in this M γγ region. Therefore, dσ N LO /dM γγ has a single-side cusp at M γγ = M LO dir (rather than a double-side cusp as in the case of asymmetric p T cuts).
The NLO quantitative results in Fig. 20-left for the symmetric p T cut configuration are consistent with the analytic behaviour in Eq. (40). We also note ( Fig. 20-left) that the NLO value of dσ/dM γγ at M γγ = M LO dir = 44 GeV is very small and, in particular, this implies that b 0 in Eq. (40) is much smaller than a 0 in Eq. (35). This is not unexpected since in the case of symmetric p T cuts (p H = p S ) the LO threshold M LO dir and the 'approximate' threshold M LO for hard radiation (see Eq. (22)) actually coincides. As a consequence of the small value of b 0 , the upward singleside cusp drives dσ N LO /dM γγ to negative values in the region just above the LO threshold. By further increasing M γγ , the physical (positive) behaviour of dσ/dM γγ sets in and, consequently, dσ N LO /dM γγ has a local minimum (with a negative value) in the vicinity of M LO dir .
The NLO behaviours in Eqs. (35) and (40) (though they are partly different) are directly related. Indeed, they are both due to the non-smooth behaviour of dσ LO /dM γγ at M γγ = M LO dir that produces an unbalance between real and virtual soft-gluon effects at higher perturbative orders [60]. The mechanism that leads to both Eqs. (35) and (40) is analogous (see the discussion that accompanies Eq. (35) in Sect. 2.3.3) and the main important difference regards the role of the real soft-emission contribution in the process of Eq. (36) (owing to transverse-momentum conservation the one-loop virtual contribution to the process qq → γγ is independent of the value of p S and, consequently, it is insensitive to the difference between asymmetric and symmetric p T cuts). The radiated soft gluon produces a transverse-momentum unbalance (p hard T γ = p sof t T γ ) between the two photons and it can preferably lead to either a decrease of p sof t T γ or an increase of p hard T γ with respect to the LO configuration (p hard T γ = p sof t T γ = p H ) depending on the phase space that is available in the presence of the photon p T cuts. In the case of asymmetric p T cuts and M γγ ∼ M LO dir , the soft-gluon momentum recoil is 'absorbed' by the softer photon (i.e., the value of p sof t T γ decreases below its LO value p sof t T γ = p H , whereas p hard T γ ≃ p H ) and this produces diphoton events with M γγ < M LO dir (see the accompanying comments to Eq. (36)). In the case of symmetric p T cuts, the momentum p sof t T γ of the softer photon cannot decrease since it is always constrained to be larger than p H = p S (as in an LO configuration) and, consequently, the soft-gluon momentum recoil is necessarily 'absorbed' by the harder photon § . The momentum p hard T γ tends to increase § This observation is also relevant for our subsequent discussion of the results in Fig. 20-right. above its threshold value p H (p hard T γ > p H ), whereas p sof t T γ ≃ p H , and this leads to diphoton events with M γγ > M LO dir . Therefore, below the LO threshold (M γγ < M LO dir ) there are no dominant real and virtual soft-gluon effects (here dσ N LO /dM γγ is smooth since it only receives contributions from hard radiation), whereas just above the LO threshold (M γγ > M LO dir ) virtual soft-gluon effects dominate since real soft-gluon radiation tends to produce diphoton events that are sufficiently far from M γγ = M LO dir . This real-virtual kinematical mismatch produces the NLO double-logarithmic enhancement (see Eq. (40)) of the non-smooth behaviour of dσ/dM γγ at M γγ ∼ > M LO dir .
In the case of symmetric p T cuts, the single-side cusp behaviour of dσ/dM γγ at M γγ ∼ M LO dir occurs at each perturbative order. Just above the LO threshold the slope of dσ/dM γγ alternatively diverges to +∞ or −∞ at subsequent perturbative orders [60]. This non-smooth order-by-order behaviour is removed by all-order-resummation of soft-gluon effects [60]. After resummation, dσ/dM γγ rapidly increases just above the LO threshold but it has a finite slope at M γγ ∼ M LO dir . In the context of fixed-order calculations, the unphysical shape of dσ/dM γγ at M γγ ∼ M LO dir is localized in a mass region whose size tends to decrease by increasing the perturbative order. This is consistent with the results in Fig. 20-left. At the NNLO the behaviour of dσ/dM γγ is qualitatively consistent with physical expectations, with the sole exception of a narrow region very close to M γγ = 44 GeV, where the divergent (to +∞) slope of dσ N N LO /dM γγ produces an excess of diphoton events that is eventually also responsible for the (quantitatively small) unphysical (increasing) behaviour of the NNLO total cross section as ∆p T → 0 ( Fig. 19-right).
We comment on the results in Fig. 20-right for the transverse-momentum spectra dσ/dp sof t Fig. 19-right. Indeed, the total cross section can be obtained by integrating dσ/dp hard sections are affected in different ways. The results in Fig. 20 show that dσ/dM γγ and dσ/dp hard T γ have evident unphysical behaviour in limited regions of M γγ and p hard T γ , respectively. In the case of dσ/d cos θ * we expect 'unphysical' effects of normalization (and, possibly, shape) in the region of central values of cos θ * , which is mostly sensitive to the M γγ region that is relatively close to M LO dir . The case of the differential cross section dσ/dp sof t T γ in the region where p sof t T γ > p H is somehow 'special', since its fixed-order results do not show evident signs of unphysical behaviour (Fig. 20-right). However, the integral of dσ/dp sof t T γ over the region with p sof t T γ ≥ p min T exactly corresponds to the total cross section in a configuration of symmetric p T cuts with p H = p S = p min T . Since the fixed-order result of such total cross section is unphysical, it turns out that the overall normalization (independently of the detailed shape) of dσ/dp sof t T γ is sensitive to soft-gluon perturbative instabilities at all values of p sof t T γ with p sof t T γ > p H .
As discussed in Sect. 2.2 the NNLO results presented in this paper are affected by a systematic uncertainty related to the finite value of the parameter (q T cut or r cut ) that is used in the numerical implementation of the q T subtraction method. An estimate of these uncertainties for total cross sections is explicitly reported in the NNLO results presented in Table 2. At fixed value of q T cut (r cut ), the systematic uncertainties on differential cross sections tend to be larger than the corresponding uncertainties on fiducial (total) cross sections. Such uncertainties can be enhanced in the presence of soft-gluon instabilities (as in the case of observables that we have discussed in this subsection), since the instabilities are smeared by the finite value of q T cut (r cut ). More detailed studies are needed to assess the numerical precision of the NNLO results in these situations. Such studies, however, cannot improve the physical predictivity of the NNLO result, since the latter is affected by sizeable theoretical uncertainties produced by the same (unphysical) soft-gluon effects.

Summary
In this paper we have considered diphoton production in hadronic collisions at LHC energies and we have presented a study of QCD radiative corrections at NLO and NNLO. At NLO we have performed a thorough analysis of photon isolation and its perturbative QCD effects, comparing results obtained by using smooth cone and standard cone isolation. While the former facilitates theoretical calculations by removing the fragmentation component (and the ensuing effects of the poorly-known non-perturbative fragmentation functions of the photon), the latter is the isolation criterion that is typically used in experimental analyses. We have then extended our study to NNLO, where only smooth cone isolation results can be presently obtained.
Our main results can be summarised as follows.
• We have shown that lowest-order results for diphoton production are affected by large radiative corrections at higher orders. The large radiative corrections are due to partonic subprocesses with high parton multiplicity in the final state. As a consequence, the radiative corrections are typically positive and they enhance diphoton production rates and related kinematical distributions. Since the final-state partons can be produced outside the photon isolation cones, it also follows that the large radiative corrections have a relatively-mild dependence on the photon isolation prescription. In particular, the quantitative differences between smooth and standard isolation tend to decrease upon the inclusion of radiative corrections.
• We have presented a detailed comparison of fiducial cross sections and differential distributions obtained within standard and smooth cone isolation at NLO. In the case of tight isolation (E T max ∼ few GeV) the comparison shows that the two isolation procedures lead to results that are consistent within the corresponding scale uncertainties. In the case of moderate isolation (E T max ∼ 10 GeV), the same features hold true for diphoton observables that are 'effectively' computed at NLO accuracy. In 'effective' LO regions (such as at small values of M γγ or ∆Φ γγ ) the two isolation procedures lead to NLO differences that are much smaller than the size of the NNLO corrections computed within smooth isolation. We thus conclude that the smooth cone isolation provides a consistent theoretical framework to compute radiative corrections up to NNLO and that, if photon isolation is sufficiently tight, the ensuing predictions can be reliably compared with experimental measurements carried out by using standard cone isolation with the same values of the isolation parameters (E T max and cone isolation radius R).
• An alternative approximation scheme of standard cone isolation consists in complementing the NLO calculation of the direct component with the LO computation of the fragmentation component. We have shown that such a scheme features an unphysical NLO dependence on E T max and, therefore, it is not recommended.
• The NNLO computation of fiducial and differential cross sections shows that the NNLO corrections are rather large in the phase space regions where the calculation is an 'effective' NNLO prediction. The impact of the NNLO corrections can become huge in phase space regions where the calculation is 'effectively' an NLO prediction, as for the cases of small values of M γγ , low values of ∆Φ γγ and relatively-large values of p T γγ . We have also observed that NLO and NNLO uncertainties obtained through scale variations do not overlap. This is mainly due to the significant contribution of the qg initial-state partonic channel, which, although parametrically suppressed by one power of α S with respect to the LO qq contribution, is enhanced by the large L qg luminosity. As a consequence, the true theoretical uncertainty is larger than the one obtained by performing customary scale variations. A more reliable estimate of the perturbative uncertainties can be obtained by properly taking into account the differences between the NLO and NNLO results. For instance, the perturbative uncertainty at NNLO (NLO) can be defined as the half difference between the central scale predictions at NNLO (NLO) and NLO (LO).
• The comparison of the NNLO calculation to the experimental results shows a clear improvement with respect to NLO in the description of the LHC data. The data tend to overshoot the NNLO predictions but, if the perturbative uncertainty is properly taken into account, the LHC data are consistent with the NNLO results in both 'effective' NNLO and 'effective' NLO regions.
• We have discussed and remarked how diphoton rates and kinematical distributions are strongly affected by the selection cuts that are typically applied on the transverse momenta of the photons. These selection cuts also lead to unphysical thresholds and ensuing perturbative instabilities in the fixed-order computations of the M γγ distribution, of the transverse-momentum spectra of the photons and of related observables. We have discussed the behaviour of the NLO and NNLO results in these threshold regions, by presenting the logarithmic structure of the perturbative instabilities in analytic form. The effect of the perturbative instabilities is tamed by considering sufficiently smeared observables. Such fixed-order perturbative instabilities can be eliminated only through a proper all-order resummation of the logarithmically-enhanced contributions.