Diphoton production at hadron colliders: transverse-momentum resummation at next-to-next-to-leading logarithmic accuracy

We consider the transverse-momentum (qT) distribution of a diphoton pair produced in hadron collisions. At small values of qT , we resum the logarithmically-enhanced perturbative QCD contributions up to next-to-next-to-leading logarithmic accuracy. At intermediate and large values of qT, we consistently combine resummation with the known next-to-leading order perturbative result. All perturbative terms up to order \alpha_S^2 are included in our computation which, after integration over qT, reproduces the known next-to-next-to-leading order result for the diphoton pair production total cross section. We present a comparison with LHC data and an estimate of the perturbative accuracy of the theoretical calculation by performing the corresponding variation of scales. In general we observe that the effect of the resummation is not only to recover the predictivity of the calculation at small transverse momentum, but also to improve substantially the agreement with the experimental data.


Introduction
Diphoton production at hadronic colliders is a very relevant process, both from the point of view of testing the Standard Model (SM) predictions [1,2,3,4,5,6] as for new physics searches.
While a general understanding of QCD processes at hadronic collisions poses a serious challenge due the complicated environment, direct or prompt photons provide an ideal probe since they constitute a theoretically and experimentally clean final state. From the theory side, because they do not have QCD interactions with other final state particles and, from the experimental side, because their energies and momenta can be measured with high precision in modern electromagnetic calorimeters.
Besides purely QCD-related considerations, diphoton final states have played a crucial role in the recent discovery of a new boson at the LHC [7,8], whose properties are compatible with those of the SM Higgs. They are also important in many new physics scenarios [9,10], in particular in the search for extra-dimensions [11] or supersymmetry [12].
In this paper we are interested in the process pp → γγX, and in particular in the transversemomentum (q T ) spectrum of the diphoton pair. The lowest-order process (O(α 0 S )) occurs via the quark annihilation subprocess qq → γγ. The QCD corrections at the first order in the strong coupling α S involve quark annihilation and a new partonic channel, via the subprocess qg → γγq. First order corrections have been computed and implemented in several fully-differential Monte Carlo codes [13,14,15,16]. At the second order in the strong coupling α S the gg channel starts to contribute, and the large gluon-gluon luminosity makes this channel potentially sizeable.
The amplitudes needed to evaluate the corrections at the second order in the strong coupling α S , for diphoton production, have been presented in [17,18,19,20], and first put together in a complete and consistent O(α 2 S ) calculation in the 2γNNLO code [21]. The next-order gluonic corrections to the box contribution (which are part of the N 3 LO QCD corrections to diphoton production) were also computed in ref. [14] and found to have a moderate quantitative effect.
The calculation of the q T spectrum poses an additional challenge with respect to more inclusive calculations, such as the total cross section. In the large-q T region (q T ∼ M γγ ), where the transverse momentum is of the order of the diphoton invariant mass M γγ , calculations based on the truncation of the perturbative series at a fixed order in α S are theoretically justified. In this region, the QCD radiative corrections are known up to the next-to-leading order (NLO), including the corresponding partonic scattering amplitudes with X = 2 partons (at the tree level [18]) and the partonic scattering amplitudes with X = 1 parton (up to the one-loop level [19]). We remind the reader that at least one additional parton is needed in order to have q T = 0 for the diphoton pair. The q T spectrun of the diphoton pair has been calculated in fully-differential Monte Carlo codes at LO [13,14,15,16] and at NLO [21,22,23]. Recently, first calculations for diphoton production in association with two [24,25,26] and three [26] jets at NLO became available.
The bulk of the diphoton events is produced in the small-q T region (q T ≪ M γγ ), where the convergence of the fixed-order expansion is spoiled by the presence of large logarithmic terms, α n S ln m (M 2 γγ /q 2 T ). In order to obtain reliable predictions these logarithmically-enhanced terms have to be systematically resummed to all perturbative orders [27]- [38]. The resummed calculation, valid at small values of q T , and the fixed-order one at large q T have then to be consistently matched to obtain a pQCD prediction for the entire range of transverse momenta.
We use the transverse-momentum resummation formalism proposed in Refs. [37,39,40] (see also [41] for processes initiated by gg annihilation). The formalism is valid for a generic process in which a high-mass system of non strongly-interacting particles is produced in hadron-hadron collisions. The method has so far been applied to the production of the Standard Model (SM) Higgs boson [39,40,42,43,44], Higgs boson production in bottom quark annihilation [45], Higgs boson production via gluon fusion in the MSSM [46], single vector bosons at NLL+LO [47] and at NNLL+NLO [48], W W [49,50] and ZZ [51] pairs, slepton pairs [52], and DY lepton pairs in polarized collisions [53].
Finally, note that besides the direct photon production from the hard subprocess, photons can also arise from the fragmentation of QCD partons. The computation of fragmentation subprocesses requires (the poorly known) non-perturbative information, in the form of parton fragmentation functions of the photon (the complete single-and double-fragmentation contributions are implemented in DIPHOX [13] for diphoton production at the first order in α S ). However, the effect of the fragmentation contributions is sizeably reduced by the photon isolation criteria that are necessarily applied in hadron collider experiments to suppress the very large irreducible background (e.g., photons that are faked by jets or produced by hadron decays). Two such criteria are the socalled "standard" cone isolation and the "smooth" cone isolation proposed by Frixione [54]. The standard cone isolation is easily implemented in experiments, but it only suppresses a fraction of the fragmentation contribution. By contrast, the smooth cone isolation (formally) eliminates the entire fragmentation contribution. For all of the results presented in this paper we rely on the smooth isolation prescription, which, for the parameters used in the experimental analysis reproduces the standard result within a 1% accuracy [55].
The paper is organized as follows. In Sect. 2 we briefly review the resummation formalism of Refs. [37,39,40]. In Sect. 3 we present numerical results and we comment on their comparison with the LHC data [6]. We also study the scale dependence of our results with the purpose of estimating the corresponding perturbative uncertainty. In Sect. 4 we summarize our results.

Transverse-momentum resummation
We briefly recall the main points of the transverse-momentum resummation formalism of Refs. [37,39,40], referring to the original papers for the full details. The formalism is general, as long as the measured final state is composed of non strongly-interacting particles (transverse-momentum resummation for strongly-interacting final states, such as heavy-quark production, has been developed in Refs. [56,57]). Here we specialize to the case of diphoton production only for ease of reading. The inclusive hard-scattering process considered is where h 1 and h 2 are the colliding hadrons with momenta p 1 and p 2 , γγ is the diphoton pair with invariant mass M γγ , transverse momentum q T and rapidity y, and X is an arbitrary and undetected final state.
The corresponding fully differential cross section, in q T , M γγ and y, which we denote for simplic-ity (since our focus is on the q T distribution) by dσ γγ /dq 2 T , can be written using the factorization formula as (up to power-suppressed corrections), where the f a/h (x, µ 2 F ) (a = q,q, g) are the parton densities of the hadron h at the factorization scale µ F , α S ≡ α S (µ 2 R ), dσ γγ ab /dq 2 T is the pQCD partonic cross section, s (ŝ = x 1 x 2 s) is the square of the hadronic (partonic) centre-of-mass energy, and µ R is the renormalization scale.
In the region where q T ∼ M γγ the QCD perturbative series is controlled by a small expansion parameter, α S (M γγ ), and a fixed-order calculation of the partonic cross section is theoretically justified. In this region, the QCD radiative corrections are known up to next-to-leading order (NLO) [17,18,19,20].
In the small-q T region (q T ≪ M γγ ), the convergence of the fixed-order perturbative expansion is spoiled by the presence of powers of large logarithmic terms, α n S ln m (M 2 γγ /q 2 T ). To obtain reliable predictions these terms have to be resummed to all orders.
To perform the resummation, we start by decomposing the partonic cross section as The first term on the right-hand side contains all the logarithmically-enhanced contributions, which have to be resummed to all orders in α S , while the second term is free of such contributions and can thus be evaluated at fixed order in perturbation theory. Using the Fourier transformation between the conjugate variables q T and b (b is the impact parameter), the resummed component dσ where J 0 (x) is the 0th-order Bessel function. The form factor W γγ is best expressed in terms of its double Mellin moments W γγ N 1 N 2 , taken with respect to the variables z 1 , z 2 at fixed M γγ , with the resummation structure of W γγ N 1 N 2 can be organized in an exponential form † were we have defined the logarithmic expansion parameter L ≡ ln(µ 2 res b 2 /b 2 0 ), and b 0 = 2e −γ E (γ E = 0.5772... is the Euler number). The scale µ res (µ res ∼ M γγ ), which appears on the right-hand side † For the sake of simplicity we consider here only the case of the diagonal terms in the flavour space of the partonic indices a, b. For a detailed discussion, we refer to Ref. [39,40].
of Eq. (6), is the resummation scale [39]. Although W γγ N 1 N 2 (i.e., the product H γγ N 1 N 2 × exp{G N 1 N 2 }) does not depend on µ res when evaluated to all perturbative orders, its explicit dependence on µ res appears when W γγ N 1 N 2 is computed by truncation of the resummed expression at some level of logarithmic accuracy (see Eq. (7) below). Variations of µ res around M γγ can thus be used to estimate the size of yet uncalculated higher-order logarithmic contributions.
The form factor exp{G N 1 N 2 } is universal ‡ and contains all the terms that order-by-order in α S are logarithmically divergent as b → ∞ (or, equivalently, q T → 0). The resummed logarithmic expansion of the exponent G N 1 N 2 is defined as follows: where the term L g (1) collects the leading logarithmic (LL) O(α p+n s L n+1 ) contributions, the function g (2) N 1 N 2 includes the next-to-leading leading logarithmic (NLL) O(α p+n s L n ) contributions [32], g N 1 N 2 controls the NNLL O(α p+n s L n−1 ) terms [33,34,58,59] and so forth; p is the number of powers of α s in the LO (Born) process. In Eq. (7), α S L is formally of order 1, so there is an explicit O(α S ) suppression between different logarithmic orders. The explicit form of the functions g (1) , g (2) N 1 N 2 and g (3) N 1 N 2 can be found in Ref. [39]. The process dependent function H γγ N 1 N 2 does not depend on the impact parameter b and it includes all the perturbative terms that behave as constants as b → ∞. It can thus be expanded in powers of α S : where σ (0) γγ is the partonic cross section at the Born level. Since the formalism applies to non strongly-interacting final states, in general the Born cross-section can only correspond to a qq or gg initial state. In the specific case of the diphoton production, both channels contribute, but at different orders in α s : the qq subprocess initiates as a pure QED process (O(α s ) 0 ), while the gg one requires a fermion loop, starting at O(α s ) 2 .
In the present work, we keep contributions up to an uniform order in α s (and all orders in α S L), namely up to α n S L n−1 . For the qq channel, this requires the inclusion of the H coefficients of Eq. (8) up to order 2: the first-order coefficients H γγ(1) N 1 N 2 are known since a long time [58], while the second-order coefficients H γγ(2) N 1 N 2 were computed only relatively recently [21,38]. For the gg channel, it is sufficient to include the leading H contribution (that is, the Born cross-section) and the appropriate G in the exponential of Eq. (7). Since it does not require any additional numerical effort, we decided, in all the plots presented in the paper, to include all the terms up to g (3) N 1 N 2 in the exponential G factor also for this channel. In this way, we technically include some terms which are of higher order in α s with respect to those in the qq channel; however we checked that those terms result in a negligible numerical effect (at 1% accuracy), that is, the difference produced by including the higher order terms is within the error bands obtained by the scale variations, which verifies the stability of the calculation. ‡ The form factor does not depend on the final state; all the hard-scattering processes initiated by qq (gg) annihilation have the same form factor.
Within a straightforward ('naive') implementation of Eq. (6), the resummation of the large logarithmic contributions would affect not only the small-q T region, but also the region of large values of q T . This can easily be understood by observing that the logarithmic expansion parameter L diverges also when b → 0. To reduce the impact of unjustified higher-order contributions in the large-q T region, the logarithmic variable L in Eq. (6) is actually replaced byL ≡ ln (µ 2 res b 2 /b 2 0 + 1) [39,42]. This (unitarity related) replacement has an additional and relevant consequence: after inclusion of the finite component (see Eq. (9)), we exactly recover the fixed-order perturbative value of the total cross section upon integration of the q T distribution over q T (i.e., the resummed terms give a vanishing contribution upon integration over q T ).
We now turn to consider the finite component of the transverse-momentum cross section (see Eq. (3)). Since dσ (fin.) γγ does not contain large logarithmic terms in the small-q T region, it can be evaluated by truncation of the perturbative series at a given fixed order. In practice, the finite component is computed starting from the usual fixed-order perturbative truncation of the partonic cross section and subtracting the expansion of the resummed part at the same perturbative order.
Introducing the subscript f.o. to denote the perturbative truncation of the various terms, we have: This matching procedure between resummed and finite contributions guarantees to achieve uniform theoretical accuracy over the region from small to intermediate values of transverse momenta.
At large values of q T , the resummation (and matching) procedure is eventually superseded by the customary fixed-order calculations (their theoretical accuracy in the large-q T region cannot be improved by resummation of the logarithmic terms that dominate in the small-q T region).
In summary, the inclusion of the functions g (1) , g (1) N 1 N 2 in the resummed component, together with the evaluation of the finite component at LO (i.e. at O(α S )), allows us to perform the resummation at NLL+LO accuracy. This is the theoretical accuracy used in previous studies [16,60,61,62] of the diphoton q T distribution. Including also the functions g Using the H γγ (2) N 1 N 2 coefficient [21,38], we are thus able to present the complete result for the diphoton q T -distribution up to NNLL+NLO accuracy. We point out that the NNLL+NLO (NLL+LO) result includes the full NNLO (NLO) perturbative contribution in the small-q T region. In particular, the NNLO (NLO) result for the total cross section is exactly recovered upon integration over q T of the differential cross section dσ γγ /dq T at NNLL+NLO (NLL+LO) accuracy.
We conclude this section with some comments on the numerical implementation of our calculation. Within our formalism, the resummation factor W γγ N 1 N 2 (b, M γγ ) is directly defined, at fixed M γγ , in the space of the conjugate variables b and N 1 , N 2 . To obtain the hadronic cross section, we have to perform inverse integral transformations: the Bessel transformation in Eq. (4) and an inverse Mellin transformation, implemented following the prescription introduced in Ref. [63]. These integrals are carried out numerically. The Mellin inversion requires the numerical evaluation of some basic N-moment functions that appear in the expression of the the second-order coefficients H γγ (2) N 1 N 2 [21,38]: this evaluation has to be performed for complex values of N 1 and N 2 ; to evaluate some of the needed special function at complex value, we use the numerical routines of Ref. [64]. We recall [39] that the resummed form factor exp{G . We avoid this singularity by introducing a smooth effective cut-off at small b, which is shown to have a negligible effect in the final result.
It is known that at small values of q T , the perturbative QCD approach has to be supplemented with non-perturbative contributions, since they become relevant as q T decreases. A discussion on non-perturbative effects on the q T distribution is presented in Ref. [39], and related quantitative results are shown in Sect. 3.

Numerical results for diphoton production at the LHC
In this section we consider diphoton production in pp collisions at LHC energies ( √ s = 7 TeV).
We present our resummed results at NNLL+NLO accuracy, and compare them with NLL+LO predictions and with available LHC data [6]. Since the present formulation of the q T resummation formalism, is restricted to the production of colourless systems F , it does not treat parton fragmentation subprocesses (here F includes one or two coloured partons that fragment); therefore, we concentrate on the direct production of diphotons, and we rely on the smooth cone isolation criterion proposed by Frixione [54] (see also Ref. [65,66]) which is defined by requesting with a suitable choice for the function χ(r). This function has to vanish smoothly when its argument goes to zero (χ(r) → 0 , if r → 0 ), and it has to verify 0 < χ(r) < 1, if 0 < r < R . One possible choice is where n is typically chosen as n = 1. This condition implies that, closer to the photon, less hadronic activity is allowed inside the cone. At r = 0, when the parton and the photon are exactly collinear, the energy deposited inside the cone is required to be exactly equal to zero, and the fragmentation component (which is a purely collinear phenomenon in perturbative QCD) vanishes completely. Since no region of the phase space is forbidden, the cancellation of soft gluon effects takes place as in ordinary infrared-safe cross sections. That is the main advantage of this criterion: it eliminates all the fragmentation component in an infrared-safe way. By contrast, it can not be implemented within the usual experimental conditions; the standard way of implementing isolation in experiments is to use the prescription of Eq. (11) with a constant χ(r) = 1. In any case, from a purely pragmatic point of view, it has been recently shown [55] that if the isolation parameters are tight enough (e.g., E T max < 6 GeV, R = 0.4), the standard and the smooth cone isolation prescription coincide at the 1% level, which is well within the theoretical uncertainty of our predictions.
The acceptance criteria used in this analysis ( √ s = 7 TeV) are those implemented by the ATLAS collaboration analysis [6]; in all the numerical results presented in this paper, we require p harder T ≥ 25 GeV, p softer T ≥ 22 GeV, and we restrict the rapidity of both photons to satisfy |y γ | < 1.37 and 1.52 < |y γ | ≤ 2.37. The isolation parameters are set to the values E T max = 4 GeV, n = 1 and R = 0.4, and the minimum angular separation between the two photons is R γγ = 0.4.
We use the Martin-Stirling-Thorne-Watt (MSTW) 2008 [67] sets of parton distributions, with densities and α S evaluated at each corresponding order (i.e., we use (n + 1)-loop α S at N n LO, with n = 0, 1, 2), and we consider N f = 5 massless quarks/antiquarks and gluons in the initial state. The default renormalization (µ R ) and factorization (µ F ) scales are set to the value of the invariant mass of the diphoton system, µ R = µ F = M γγ , while the default resummation scale (µ res ) is set to µ res = M γγ /2. The QED coupling constant α is fixed to α = 1/137. Non-perturbative (NP) effects are expected to be important at very small q T . In this paper we follow the strategy of Ref. [39], implementing them by multiplying the b-space form factor W γγ of Eq. (4) by a 'NP factor' which consists of a gaussian smearing of the form where a denotes the initial state channel, a = F for qq and a = A for gg (as usual, C F = (N 2 c − 1)/(2N c ) and C A = N c ). In order to asses the importance of the NP contributions, we vary g N P in the interval from g N P = 0 GeV 2 (no NP contributions) to g N P = 2 GeV 2 , corresponding to moderate NP effects [39].
An additional and potentially important source of theoretical uncertainty arises from an ambiguity in the definition of the photon momenta in the resummation formalism. In fact, in the main resummation formula (4), which is used to define both the resummed and (via the subtraction Eq. (9)) the finite contributions to the partonic cross-section, the diphoton pair total transverse momentum q T is not associated with the recoil of any extra physical particles in the final state.
After q T resummation the angular distributions of the photons are still provided by the Born level functions (σ (0) γγ , H γγ N 1 N 2 ), which appear as multiplicative factors in front of the Sudakov form factor of Eq. (7). At this point there are two strategies to follow, which differ by corrections that are of O(q T /M γγ ) order-by-order in the perturbative expansion [69] (after having matched the resummed calculation with a complete N k LO calculation, these corrections start to contribute at the N k+1 LO).  12)) is also shown.
One of them is to just use the Born phase space for all the angular distributions in front of the Sudakov form factor of Eq. (7). In this case the transverse momentum q T is neglected in the calculation of the photon momenta, while the momenta q i , (i = 1, 2) of the colliding partons are given by: where x i (i = 1, 2) are the parton momentum fractions and P i (i = 1, 2) are the momenta of the colliding hadrons. The momenta q i respects the Born level kinematics where q γi , (i = 1, 2) are the momenta of the two photons in the final state. In this case the two photons are always in the back-to-back configuration, and therefore the kinematic effects of the transverse-momentum recoil, are not included in the momentum of each single photon. We call this approach the 2-body phase space; all the differential distributions which use this method have a (2-body) label.
A more elaborate, and arguably more physical, approach is to consider the effects of the transverse-momentum recoil in the two-photon final state. Therefore these kinematic effects have to be 'absorbed' by the incoming parton momenta q 1 and q 2 , in order to respect the Born level kinematics and the LO kinematics of Eq. (13). There are different consistent implementations of this approach, all of them giving equivalent (up to higher order corrections) results § . The implementation that we use in the following is essentially equivalent to define the separate photon momenta according to the Born level angular distribution computed in the Collins-Soper frame [70] of the diphoton system (in practice, this is the same procedure used in Ref. [16]). We call this method the 3-body phase space, and distributions obtained using this method are labelled (3-body). The modifications introduced by this phase space parametrization in the momentum fractions x 1 , x 2 of the incoming partons are neglected because they produce negligible effects (O(q T /M γγ )) in the cross-section. Notice that it is not the formal 3-body phase space (the 3-body effect is produced by the Lorentz transformation with finite q T of the Collins-Soper frame). We have checked that if the formal 3-body phase space is used ¶ , the changes in the cross section (comparing the 2-body and 3-body cases) are of order O(q T /M γγ ). The ambiguity in the treatment of the photon momenta is considered as an additional source of theoretical uncertainty.  (3)), is also shown for comparison (dotted line). We observe that the LO result diverges to +∞ as q T → 0, as expected.
The finite component is regular over the full q T range, it smoothly vanishes as q T → 0 and gives an important contribution to the NLL+LO result in the low-q T region. That is mostly originated by the qg channel, which starts at NLO and provides a subleading correction in terms of logs (single logarithmic terms) but contributes considerably to the cross-section due to the huge partonic luminosity compared to the formally leading qq channel. The resummation of the § The details of the implementation of these kinematic effects are discussed in a forthcoming paper [69]. ¶ Which violates the Born level kinematics of Eq. (14).
small-q T logarithms leads to a well-behaved distribution: it vanishes as q T → 0 and approaches the corresponding LO result at large values of q T .
The results in the right panel of Fig. 1 are systematically at one order higher: the q T spectrum at NNLL+NLO accuracy (solid line) is compared with the NLO result (dashed line) and with the NLO finite component of the spectrum (dotted line). The NLO result diverges to −∞ as q T → 0 and, at small values of q T , it has an unphysical peak that is produced by the compensation of negative leading and positive subleading logarithmic contributions. The contribution of the NLO finite component to the NNLL+NLO result is of the order of the 50% at the peak and becomes more important as q T increases. A similar quantitative behaviour is observed by considering the contribution of the NLO finite component to the NLO result. At large values of q T the contribution of the NLO finite component tends to the NLO result. This behaviour indicates that the logarithmic terms are no longer dominant and that the resummed calculation cannot improve upon the predictivity of the fixed-order expansion.
We also observe that the position of the peak in the NNLL+NLO q T distribution is slightly harder than the corresponding NLL+LO q T distribution. This effect is (in part) due to the large transverse-momentum dependence of the fixed order corrections.
As discussed in Sect. 2, the resummed calculation depends on the factorization and renormalization scales and on the resummation scale µ res . Our convention to compute factorization and renormalization scale uncertainties is to consider independent variations of µ F and µ R by a factor of two around the central values µ F = µ R = M γγ in independent way in order to maximise them: (µ F = 2 M γγ , µ R = M γγ /2, µ res = M γγ /2) and (µ R = 2 M γγ , µ F = M γγ /2, µ res = M γγ /2). The uncertainty due to the resummation scale variation is assessed separately by varying it between µ res = M γγ /4 and µ res = M γγ at fixed µ F and µ R . Figure 4: The q T spectrum of diphoton pairs at the LHC. The NNLL+NLO result is compared with the NLL+LO result, for the window 0 GeV < q T < 40 GeV (left panel) and the full spectra (right panel). We use the 3-body parametrization in the resummed cross-section, and set g N P = 2 GeV 2 . The bands are obtained by varying µ R and µ F as explained in the text.
In Fig. 2 we compare the impact of the various sources of theoretical error for the NNLL+NLO  Figure 5: Comparison of the theoretical prediction for the q T spectrum of diphoton pairs at the LHC to the experimental data. The NNLL+NLO result is compared with the ATLAS data of Ref. [6], for the window 0 GeV < q T < 40 GeV (left panel) and the full spectra (right panel). In the right panel the NNL+LO distribution at central scale (dotted line) is also shown in order to compare it with the data and the NNLL+NLO result. In the theoretical curves we use the 2body parametrization for the resummed cross-section and set g N P = 2 GeV 2 ; the bands (solid and dashed lines) are obtained by varying µ R and µ F as explained in the text.
predictions: the 2-body and 3-body parametrizations of the photon phase space, the effect of the variation of the non-perturbative contribution, and the variation of the factorization and renormalization scales. In the left panel of Fig. 2 we use only the central scale (µ F = µ R = M γγ , µ res = M γγ /2), and we note that the 3-body method results in a slightly larger cross section (about 10%) around the peak (q T ∼ 5 GeV) than the 2-body one. For q T > 20 GeV all the contributions in the left panel coincide, which is consistent with the fact that at these values the resummed component starts to vanishes and the fixed order result dominates the cross section, as we can anticipate from Fig. 1. Also we note that if we use a NP parameter, the peak of the distribution is located at larger values of q T (q T ∼ 5 GeV) and the shape of the distribution is slightly different (for values 0 GeV< q T < 22 GeV) from the case in which the NP parameter is not implemented. These differences, which are stronger for q T < 10 GeV, have their origin in the resummed component, which is the only contribution that depends on g N P . In the right panel of Fig. 2 we show the comparison between the variation of the scales in the 3-body approach with the central scale results of the 2-body and 3-body frameworks (with and without non perturbative parameter, respectively). Evidently, the uncertainty due to the ambiguity of the parametrization of the photon momenta turns out to be subdominant with respect to the one arising from the variation of the scales which provides, by itself, a reasonable estimate of theoretical uncertainties.
In Fig. 3 we show the NNLL+NLO transverse momentum distribution for three different implementations of the µ res parameter (µ res = M γγ /4; M γγ /2; M γγ ) at fixed µ F = µ R = M γγ . The 2-body phase space and a non-perturbative parameter g N P = 2 GeV 2 were used. We notice the small impact of the variation of the µ res scale in the cross section (at per-cent level). In the left panel of Fig. 3 we present the transverse momentum distribution for values of q T within the We vary the scales from (µ F = 2M γγ ; µ R = M γγ /2) to (µ R = 2M γγ ; µ F = M γγ /2) at fixed µ res = M γγ /2.  Figure 6: Comparison of the theoretical prediction for the q T spectrum of diphoton pairs at the LHC with the experimental data. The NNLL+NLO result is compared with the ATLAS data of Ref. [6], for the window 0 GeV < q T < 40 GeV (left panel) and the full spectra (right panel). In the right panel the NNL+LO distribution at central scale (dotted line) is also shown in order to compare it with the data and the NNLL+NLO result. In the theoretical curves we use the 3-body parametrization for the resummed cross-section and set g N P = 2 GeV 2 . The bands (solid and dashed lines) are obtained by varying µ R and µ F as explained in the text.
interval 0 GeV< q T < 40 GeV, and in the right panel of Fig. 3 the full spectra. We also notice that the strongest effect of the variation of the µ res scale appears in the last bin of right panel of Fig. 3. This is expected since the resummation scale effectively sets the value of transverse momentum at which the logarithms are dominant. A choice of a very large resummation scale affects the distribution at larger transverse momentum and might in general result in a mismatch with the fixed order prediction due to the artificial introduction of unphysically large logarithmic contributions in that region. Similar results are obtained if the 3-body phase space is used instead of the 2-body one.
In Fig. 4 we compare the variation of the scales of the NNLL+NLO and NLL+LO predictions (3-body phase space), for the interval 0 GeV< q T < 40 GeV (left panel) and the full spectra (right panel). We notice that the dependence on the scales is not reduced when going from NLL+LO to NNLL+NLO. This is mostly because at NNLL+NLO a new channel (gg) opens, in which the box contribution (effectively "LO" but formally O(α 2 S )) ruins the reduction of the scale dependence usually expected when adding second order corrections for the qq channel and first order corrections for the qg channel. This effect has the same origin that the reported behaviour of the diphoton production at NNLO of Ref. [21], when the variation of the scales is implemented. Since NNLL+NLO is the first order at which all partonic channels contribute, it is possible to argue that this is the first order at which estimates of theoretical uncertainties through scale variations can be considered as reliable. The same results are obtained if the 2-body phase space is used instead of the 3-body one.
The peak observed in the right panel of Fig. 4 is the so called Guillet shoulder [68], which is a real radiation effect and has its origin in the fixed order contribution. It appears stronger in the NNLL+NLO q T distribution than in the NLL+LO, due to the larger size of the real contributions In both panels, we compare the resummed prediction using the 3-body parametrization in the resummed cross-section (g N P = 2 GeV 2 ) with the fixed order prediction.
at NLO.
In the small q T region (q T < 4 GeV) the real radiation effects are no longer dominant in the q T distribution. The finite and resummed component (which vanish as q T → 0) are of the same order for q T < 4 GeV. In absence of the strong real radiation effects the contributions are almost completely Born like. This is the main reason why the NNLL+NLO bands overlaps with those at the previous order.
In Figs. 5 and 6 we compare the LHC data ( √ s = 7 TeV) from ATLAS [6] with our resummed theoretical predictions (at NNLL+NLO and NLL+LO) using the 2-body and 3-body approaches, respectively. In both cases we use a NP parameter different from zero (g N P = 2 GeV 2 ) and we estimate the theoretical uncertainty by the variation of the µ R and µ F scales. In the left panels we show the q T distribution in the window (0 GeV< q T < 40 GeV), while in the right ones we show the full spectra in logarithmic scale.
We observe in general an excellent agreement between the resummed NNLL+NLO prediction and the experimental data, that is accurately described within the theoretical uncertainty bands in the whole kinematic range. Also we observe that the NLL+LO result is not enough to describe the phenomenology of the transverse-momentum distribution of the LHC data (Figs. 5 and 6, right panel). By direct comparison to the fixed order prediction, we notice that the effect of resummation is not only to recover the predictivity of the calculation at small transverse momentum, but also to improve substantially the agreement with LHC data [6].
While the resummation performed in this work reaches NNLL accuracy formally only for the diphoton transverse momentum distribution, its predictions can be extended to other observables as well, since at least the leading logarithmic contributions have a common origin from soft and collinear emission. Note that the 3-body approach is more suitable for a consistent implementation of these leading logarithmic effects on the observables. In Fig. 7 we show results on more exclusive observables: the q T distributions of the harder (left-hand plot) and softer (right-hand plot) photon at the central scale. We compare the results obtained with the 2-body and 3-body phase space approaches, using a non-perturbative parameter g N P = 2 GeV 2 with the fixed order result at NNLO for diphoton production [21].
The 2-body phase space transverse-momentum distribution at NNLL+NLO provides the same result than the NNLO fixed order cross section for diphoton production [21]. This is consistent with two following related facts: i) the single photon momentum does not carry any information about the recoil due to the transverse momentum q T in the 2-body approach; ii) the NNLO (NLO) result for the total cross section is exactly recovered upon integration over q T of the differential cross section dσ γγ /dq T at NNLL+NLO (NLL+LO) accuracy.
The effects of resummation are only present in these more exclusive observables if the recoil due to q T are absorbed by the photons in the final state * * , which is equivalent to the implementation of a 3-body like phase space. In this way, because the single photon momentum depends on q T , the integral over q T of the NNLL+NLO (NLL+LO) distribution does not recover the NNLO (NLO) result as we can observe in Fig. 7. Here the last integral over the single photon momentum of the NNLL+NLO (NLL+LO) distribution, is required in order to recover the NNLO (NLO) result for the total cross section.
We comment on the q T distribution of the softer photon in the region around the back-to-back threshold q S T γ ∼ 25 GeV (see Fig. 7 right panel). The NLO fixed order result has a step-like behaviour, and this necessarily produces [71] integrable logarithmic singularities at each subsequent perturbative order. The peak of the NLO fixed order distribution at q S T γ ∼ 25 GeV is an artifact of these perturbative instabilities. The instability is cured by all-order perturbative resummation, which leads to a smooth q T distribution with a shoulder-like behaviour [71] in the vicinity of the back-to-back threshold.
In Fig. 8 we present the results of the cross section as a function of the azimuthal angle ∆Φ γγ . In the left panel we compare the fixed order (NLO), finite (NLO) and full (NNLL+NLO) ∆Φ γγ distributions. The fixed order component dominates the cross section over the whole ∆Φ γγ range. However, as could be expected, the effect of resummation is stronger for kinematic configurations near the ∆Φ γγ ∼ π which correspond to q T ∼ 0 GeV. As in the case of the fixed order q T distribution, the ∆Φ γγ fixed order differential cross section is not well-behaved near the back-toback configuration: it actually diverges as ∆Φ γγ → π (q T → 0). The finite contribution (Eq. (9)) is well-behaved near the back-to-back configuration, and the full result (NNLL+NLO) improves the description in the region near ∆Φ γγ ∼ 0.
In the right panel of Fig. 8 we compare our theoretical prediction at NNLL+NLO level of accuracy with the LHC data [6] using the variation of the µ R and µ F scales to estimate the theoretical uncertainty. We observe that the transverse momentum resummation provides a better description of the data with respect to the fixed order result. In both panels of Fig. 8 we used the 3-body approach to describe the diphoton phase-space. In fact, in this kind of observables (also see Fig. 7 pp -> γγ+X ATLAS √s=7TeV Figure 8: The ∆Φ γγ distribution of diphoton pairs at the LHC. In the left panel, we show the fixed order prediction, the complete resummed prediction and just the resummed contribution at the central scale. In the right, the full NNLL+NLO result (using the extreme values for µ R and µ F to estimate the theoretical uncertainty) is compared with the ATLAS data of Ref. [6]. section.

Summary
In this paper we performed the transverse momentum resummation for diphoton production at NNLL accuracy in hadron collisions. At small values of q T , the calculation includes the resummation of all logarithmically-enhanced perturbative QCD contributions, up to next-to-next-to-leading logarithmic accuracy; at intermediate and large values of q T , it combines the resummation with the fixed next-to-leading order perturbative result. The combination is performed in such a way as to exactly reproduce the known next-to-next-to-leading order result for the total cross section; in the end, the calculation consistently includes all perturbative terms up to formal order α 2 S . The theoretical uncertainty was estimated by varying the various scales (renormalization, factorization and resummation) introduced by the formalism as well as the parametrization of the diphoton phase-space. The result was compared to experimental data, showing good agreement between theory and experiment over the whole q T range. With respect to the fixed-order calculation, the present implementation provides a better description of the data and recovers the correct physical behaviour in the small q T region, with the spectrum smoothly going to zero. The same set-up also allows the calculation of more exclusive observable distributions; the q T spectrum of the individual photons and the ∆Φ γγ distribution are given as examples.