The Diphoton $q_T$ spectrum at N$^3$LL$^\prime$+NNLO

We present a $q_T$-resummed calculation of diphoton production at order N$^3$LL$^\prime$+NNLO. To reach the primed level of accuracy we have implemented the recently published three-loop $\mathcal{O}(\alpha_s^3)$ virtual corrections in the $q\bar{q}$ channel and the three-loop transverse momentum dependent beam functions and combined them with the existing infrastructure of CuTe-MCFM, a code performing resummation at order N$^3$LL. While the primed predictions are parametrically not more accurate, one typically observes from lower orders and other processes that they are the dominant effect of the next order. We include in both the $q\bar{q}$ and loop-induced $gg$ channel the hard contributions consistently together at order $\alpha_s^3$ and find that the resummed $q\bar{q}$ channel without matching stabilizes indeed. Due to large matching corrections and large contributions and uncertainties from the $gg$ channel, the overall improvements are small though. We furthermore study the effect of hybrid-cone photon isolation and hard-scale choice on our fully matched results to describe the ATLAS 8 TeV data and find that the hybrid-cone isolation destroys agreement at small $q_T$.

The production of prompt isolated photon pairs at hadron colliders is a test of QCD with a clean experimental signature and constitutes as a background to the h → γγ decay. Differential measurements at the LHC are available at 7 TeV , both by ATLAS and CMS [1-3], at 8 TeV by ATLAS [4], and very recently at 13 TeV also by ATLAS [5]. Further, many models of new physics predict resonant large mass diphoton decays for which events with diphoton invariant masses of up to ∼ 2 TeV [6, 7] are used to constrain them.
Photons can either be produced directly or through the fragmentation of QCD partons. Photons produced through fragmentation require treatment of their singular collinear splitting. These singularities can be renormalized into non-perturbative fragmentation functions, but which are available only at NLO so far with large uncertainties [8,9]. In higher-order calculations the collinear singularity is typically removed with a smoothcone isolation prescription [10], which eliminates the collinear singularity, but does not prevent the usual cancellation of soft singularities. Quite some attention has been paid to the smooth-cone photon isolation procedure in recent years due to the uncertainty associated with the difference to experimental photon isolation [11][12][13].
NNLO calculations [14,15] for diphoton production have shown large perturbative corrections with uncertainties that significantly underestimate the difference going from NLO to NNLO, estimated by using typical scale-variation procedures. This is in part due to a new loopinduced gg channel entering at NNLO [15][16][17][18], but even true for just the Born-level induced qq channel. This has theorists led to deviate from the usual scale-variation prescription and suggest taking (half) the difference between NLO and NNLO as an estimate of the NNLO uncertainty [13]. This in turn mandates the calculation of higher-order predictions until uncertainties obtained from scale-variation stabilize. The inclusion of the recently published three-loop qq → γγ hard function [19] in the calculation presented in this paper contributes to these improvements.
Resummed predictions at small q T are available at N 3 LL [20] and at N 2 LL [21], both matched to NNLO 0 fixedorder and have been computed previously at lower order matched to NLO 0 [22][23][24][25]. 1 Recently, a matching to parton shower has been considered [26]. NLO electroweak effects have been studied in ref. [27], which are found to be less than one percent for the q T distribution below 200 GeV , so they are not relevant for the results presented in this paper.
Uncertainties in the measured diphoton transverse momentum (q T ) and φ * [28] distributions at 7 TeV and 8 TeV are about 10%. This is about the uncertainty 1 With the subscript 0 we denote that the order is with respect to the Born-level topology γγ and not with respect to γγ+jet, i.e. at large qT . Latter would be denoted with a 1 subscript. Therefore N 3 LL+NNLO0 is equivalent to N 3 LL+NLO1. This is often ambiguous in the literature. In the following we will omit the 0 subscript. estimated by scale variation from fixed-order NNLO (O(α 2 s )) and N 3 LL q T -resummed calculations. While the uncertainties seem large, the differences between theory and measurement are found to be larger. Tensions start at q T ∼ 15 GeV and rise to differences of 30% between central values for q T larger than 50 GeV [20].
Before moving on to a discussion about the photon isolation, we need to set up the relevant notation. The smooth-cone photon isolation prescription in this paper [10] restricts the transverse hadronic energy E had T around photons to be where E iso T is an isolation cone energy that can either be fixed or dependent on the photon transverse momentum, R s is the isolation cone radius and n is a parameter. In this paper we also consider a simple hybrid-cone isolation that takes the smooth-cone isolation within an inner radius R s and a fixed-cone prescription in the outer cone with radius R o : Recently it has been argued that a bulk of the datatheory tensions are an artifact of two effects [11]. The first one is regarding the hard renormalization scale choice of m γγ that has often been used for predictions. Since the diphoton pair is not produced resonantly, m γγ is not the definite obvious choice for capturing the hard process kinematic scale. An alternative studied is to take the arithmetic mean of the photon transverse momenta q γ T . The second effect is due to the photon isolation, which likely needs to be assigned larger uncertainties than previously thought. The suggestion is to use the hybrid-cone isolation, which allows for a better matching to the experimentally used fixed-cone prescription by adjusting the inner cone radius. An isolation uncertainty can then be obtained by varying the inner cone radius by some amount.
In ref. [11] it is further argued that isolation scheme and hard scale have compensating effects, at least for distributions sensitive to the photon separation ∆R γγ like the m γγ distribution. The authors show that µ = m γγ plus smooth-cone isolation and µ = q γ T plus hybrid-cone isolation should go together, respectively, as there are compensating effects in the low and large mass region. Indeed they show that the agreement between data and theory is best for the combination of q γ T plus hybridcone isolation, but the experimental uncertainties in the large and small mass regions are also largest, leaving an unclear picture.
Every isolation prescription introduces an unphysical discontinuity due to the presence of step functions in the measurement function [29,30]. At NNLO the discontinuity translates into a Sudakov singularity that has been studied in more detail in ref. [11]: In the case of the smooth-cone isolation the affected observable is not of experimental relevance. But the hybrid-cone isolation places this singularity directly into the q γγ T distribution around E iso T and its only cure is a sufficiently large experimental binning. Effectively, an overall better agreement in some distributions like m γγ and at large q T is traded for a worse agreement at small q T , φ * [28] or a T [31] as well as for azimuthal photon separations ∆Φ ∼ π. While these are all regions that need to be addressed by the resummation of large q T /Q logarithms, we show that the discontinuity effects are amplified and that the agreement with data is considerably worsened: For the q T distribution around E iso T and equivalently for Φ * , the agreement of central values within a few percent is turned into disagreement of 50-60%, see our results in the following. Indeed in ref. [11] the authors anticipated problems with slicing subtractions. The q T resummation (at leading power) in that sense acts like a slicing procedure, adding only Born-topology corrections on top of the fixed-order prediction that exhibits the Sudakov singularity. There is no compensating mechanism from the fixed-order expansion of the resummed result, resulting in the large unphysical matching corrections. The hybrid-cone isolation and natural scale choice also can unfortunately not help to address the large difference between α s and α 2 s results, but they raise valid concerns about previous assumptions.
Exactly this current situation makes a calculation of higher-order (α 3 s ) effects necessary. They can hopefully unambiguously stabilize the perturbative series to allow for truncation uncertainties that can be trusted by finding overlapping bands between different orders in both the large and small q T regions. A first step in that direction at large q T is the very recent NNLO calculation of γγ+jet [32]. This calculation predicts positive corrections at the order of 10% (without the loop-induced gg channel) below 100 GeV , which likely fills the currently seen gaps between NLO large-q T predictions and data, see e.g. fig. 18 in ref. [20]. Also, since using the hybrid-cone isolation scheme with large-q T NLO predictions leads to agreement with data within uncertainties, see ref. [11] and our plots in the following, it will be interesting to investigate how the additional 10% effects from a large-q T NNLO prediction behave in the presence of this isolation scheme.
In the present study we demonstrate the effect of the α 3 s hard functions, which are three-loop for the qq [19] channel and two-loop for the gg channel [17] (implemented in MCFM in refs. [15,16]), and incorporate them with the recently published three-loop beam functions [33][34][35]. These ingredients are commonly referred to as "constant" pieces and including them to a higher order constitutes the primed accuracy, i.e. designated as N 3 LL + NNLO 0 . The inclusion of the primed contributions is typically a dominant effect of the next order and also largely responsible for stabilization of truncation uncertainties. This can be seen by comparing for example resummed spectra at NLL with N 2 LL and N 2 LL with N 3 LL. It is particularly true when matching corrections are small, such as in Drell-Yan production. This was recently observed in refs. [36,37], and we also observe it in our N 3 LL implementation for Drell-Yan.
But while for Drell-Yan production q T resummation works up to 40-50 GeV with almost negligible (1-2%) matching corrections [20,38], the situation for photon processes is different due to the photon isolation [39]. With fiducial cuts the photon isolation prescription induces large linear power corrections [20]. Furthermore, the typical minimum q γ,1 T and q γ,2 T cuts on the two photons completely invalidate q T resummation above ∼ q γ,1 T + q γ,2 T (= 70 GeV for the ATLAS 8 TeV study) and the matching corrections quickly grow towards that point. With these effects taken together, the matching corrections from fixed-order are 50-75% over the whole range of applicable q T . This means that more than half of the cross-section at small q T comes from the terms of the fixed-order prediction and are not described by the higher-order q T -logarithms. Ideally one would like to resum the linear power corrections O(q T /Q), but the isolation makes this difficult.
While several effects diminish stabilizing effects from the inclusion of the three-loop qq hard function, the present calculation allows for a first inclusion of the three-loop virtual corrections in a physical calculation without the complication of a full N 3 LO calculation, and therefore shows directly the impact of including these corrections. We also treat for the first time the qq and gg loop-induced α 3 s hard functions fully consistently together in the resummation. Both channels have to be added separately together, of course, which means to reach the α 3 s accuracy for the "constant" part, we add the N 3 LL resummed qq channel to the N 2 LL resummed gg loop-induced channel.
Implementation. We extend the existing framework CuTe-MCFM [20] which implements N 3 LL q T resummation in the SCET formulation of refs. [40,41] matched to fixed-order NNLO 0 [15,16]. This framework achieves an accuracy of α 2 s in improved perturbation theory at small and large q T . To upgrade to N 3 LL accuracy we have implemented the one-, two-and three-loop MSrenormalized virtual amplitudes from ref. [19] and restored the renormalization-scale dependence by solving the associated RGE [42] to order α 3 s . For the numerical evaluation of harmonic polylogarithms up to weight six in the hard function we use the hplog library [43]. 2 After that we find full agreement to machine precision with the existing one-and two-loop results in MCFM. The implemented MS-renormalized amplitudes constitute the hard function. We also implemented the three-loop beam functions [33][34][35] and find that the double logarithmic L ⊥ ∼ log(x 2 T µ 2 )-dependence, where x T is the fourierconjugate of q T , is as predicted by associated renormalization group equations [40]. For the resummation we employ an improved power counting L ⊥ ∼ 1/ √ α s (relevant at small q T ), factor out the double-logarithmic L ⊥ -dependence of the beam functions and exponentiate it through associated RGE [40]. In addition to the previously published version of CuTe-MCFM [20] we have also made the resummation scale uncertainties more robust by additionally varying the rapidity scale following ref. [44].
Results. Before showing differential results, we first discuss fiducial total cross-sections. At fixed order we can calculate these up to NNLO. For the q T -resummed predictions we can simply integrate over q T to obtain a total cross-section which includes higher-order logarithmic corrections. Ideally this is within the scaleuncertainties of the fixed-order result. Since the bulk of the cross-section comes from small q T , the resummation can, in principle, improve the prediction and uncertainties. With N 3 LL q T resummation we can give a consistent prediction including both the qq-and gginitiated hard functions at order α 3 s .
In table 1 we present the cross-sections for the qq and gg hard-function initiated processes at fixed-order and by integrating resummed cross-sections over q T . Our resummation is matched to fixed-order predictions using a transition function as detailed in ref. [20]. The matching uncertainty from varying the transition function is about one percent at higher orders, i.e. small compared to the other uncertainties. We neglect matching corrections below 1 GeV , which has an effect smaller than the numerical precision quoted. For the fixed-order crosssections the numbers in front of µ R and µ F denote a variation of the renormalization and factorization scale by a factor of 2 and 1/2, respectively. We only quote the maximum of upwards and downwards variation and take these as symmetric uncertainties. We also vary combinations, but their uncertainties are found to be smaller or similar. For the resummed predictions we correlate the renormalization scale µ R with the hard scale and correlate the resummation scale with the factorization scale µ F as in ref. [20]. To obtain robust uncertainties we additionally vary the rapidity-RGE scale, but find that these uncertainties are smaller than the resummation scale uncertainties for this process. In most studies all of these uncertainties are combined by taking the envelope. Here we just quote both numbers.
At NLO and NNLO we find indeed that the integrated resummed predictions agree within uncertainties with Table 1: Cross-sections at 8 TeV for qq and gg hardfunction initiated processes with µ R = m γγ and smooth-cone isolation.
the fixed-order result. At order α 2 s (NNLO) the uncertainties of the resummed prediction are smaller. The primed prediction taking into account the α 3 s constant terms ( N 3 LL + NNLO 0 ) shows no significant changes compared to the α 2 s prediction, neither in cross-section nor in uncertainties. To understand that underlying this is still a stabilization of the cross-section we have to look at the results differentially next. Taking together the qq and gg contributions, the prediction and measurement agree at best marginally within mutual uncertainties, which is a known issue.
As suggested in ref. [11] a more natural choice for the renormalization (hard) and factorization scale is to take the arithmetic average of the photon momenta µ R = q γ T , so we consider this in addition to our standard choice of µ R = (m γγ ) 2 + (q γγ T ) 2 . With this scale we present results at the highest orders in table 2. The cross-section increases noticeably by more than the size of the previous scale variation uncertainty, and is now in much better agreement with the measurement.
Additionally we compare with a hybrid-cone implementation where the inner smooth-cone radius is reduced to R s = 0.1 and the outer radius is R o = 0.4. The value of R s = 0.1 has been suggested in ref. [11] to best match the ATLAS 8 TeV measurement and the fragmentation calculation. We therefore present results with both µ R = q γ T and an inner cone R s = 0.1 in table 3. This time the central value already overshoots the measured Table 2: Cross-sections at 8 TeV for qq and gg hardfunction initiated processes with µ R = q γ T and smooth-cone isolation. Table 3: Cross-sections at 8 TeV for qq and gg hardfunction initiated processes with µ R = q γ T and hybrid-cone isolation with inner cone angle R 0 = 0.1.
central value, but is also within uncertainties. To visualize the results of tables 1 to 3 and directly compare with the measurement, we show selected combinations in fig. 1, where the qq and loop-induced gg channel are combined consistently at the same order of α s .
The resummed q T spectrum without matching. Differentially, we first discuss the resummed q T distributions without matching. We also focus only on the qq channel first, where we now include the α 3 s hard function at N 3 LL . This is shown in fig. 2.
The bottom plot shows the ratio to the highest-order prediction N 3 LL , while the upper plot shows the absolute distributions. The N 2 LL result takes into account just the α s hard function, while N 2 LL and N 3 LL take into account the α 2 s hard function but differ in the order of RGE solution that resums the large logarithms. Only N 3 LL is fully consistent to order α 2 s in improved perturbation theory.
The displayed uncertainties are obtained from the envelope of a variation of hard scale, resummation scale and rapidity scale. Since we use the envelope, the largest un- certainty determines the band, which in our case is from the resummation scale. Below 4 GeV the uncertainties turn constant for the following reason. In any resummation formalism a cutoff at small q T is necessary since, for example, otherwise α s would be evaluated at scales where non-perturbative effects become significant, i.e. where α s becomes large. We choose to set a minimum scale of 2 GeV , which consequently leads to frozen out uncertainties below 4 GeV when a downwards scale variation becomes ineffective. Without such a cutoff the uncertainties would become arbitrarily large and would not represent realistic perturbative uncertainties, in addition to numerical problems.
Overall the uncertainties decrease going from N 2 LL to higher orders. At the smallest q T the uncertainties for N 2 LL are about 20% and reduce to 12% going towards N 3 LL, but change little between N 3 LL and N 3 LL , likely an effect due to the same order in RGE running. For larger q T we compare N 2 LL, N 2 LL and N 3 LL and find that the primed accuracy, taking into account the higherorder hard and beam functions but solving the RGEs to a lower order, is responsible for the bulk of corrections. This is reflected by the good agreement between the orange and blue lines. Consequently one expects that the inclusion of the α 3 s hard and beam functions is responsible for the bulk of corrections within a consistent N 4 LL calculation. Indeed the highest-order prediction N 3 LL is between both lower order predictions and noticeably decreases uncertainties above 10 GeV .
These finding indicate a stabilization of the qq channel, but at order α 2 s the loop-induced gg → γγ channel enters that is enhanced due to the large gluon luminosity at the LHC at small momentum fractions. In fig. 3 we include this channel at the respective orders, i.e. the α 2 s hard function at N 3 LL and the α 3 s hard function at N 3 LL . The gg channel is a substantial contribution with huge uncertainties. At low q T the uncertainties are still at the order of 10%. Towards larger q T the gg channel contributes "only" half of the cross-section, but the uncertainties are so large that the qq channel uncertainties of 1-2% seen in fig. 2 blow up to 10% in the sum of both channels. This is not unexpected since the gg channel, despite being of order α 3 s is only NLO, respectively N 2 LL+NLO accurate. To increase the precision in the intermediate low q T region of about 10-50 GeV where resummation remains relevant, threeloop α 4 s corrections to the gg channel will therefore be important.
Fully matched results. We now move on to show fully matched results and directly compare with the 8 TeV ATLAS measurement. Our transition function is a function of x = q 2 T /m 2 γγ with a parameter x max that determines the transition region, see ref. [20] for a detailed description. For the following plots we use x max = 0.1 that performs the transition mostly in the region between 20 and 50 GeV . Since m γγ is not sharply peaked as in resonant boson production, there is a tail of larger m γγ , for which the transition is later, such that we need to choose x max relatively small to prevent reaching ∼ 70 GeV where the resummation breaks down due to the given photon cuts. We estimate the matching uncertainty by varying the transition function to use x max = 0.2. This shifts the transition to be between 30 and 70 GeV . The resulting difference is small compared to our presented uncertainty bands obtained by scale variation.
So far the stabilization of the resummed qq channel has been somewhat overshadowed by large gg channel uncertainties. We now show the fully matched result in fig. 4 with ATLAS binning. The matching corrections at α 2 s for diphoton production are sizable about 50%, as can also be seen by comparing with fig. 3. The top plot shows the absolute predictions for N 3 LL+NNLO and N 3 LL +NNLO, while the lower plot shows the ratio to the higher-order prediction without the gg channel. The higher-order corrections from the qq channel are small as we have seen, but the α 3 s corrections on the gg channel have a noticeable impact (at large q T we include the α 3 s matching corrections in this channel). In both cases the uncertainties are large and transition into fixed-order uncertainties of about 15% at large q T . Overall both predictions show uniform uncertainties of 10-15%.
To decrease uncertainties noticeably we will first have to include α 3 s matching-corrections also in the qq channel, which at low q T make up about 50% of the crosssection. Second, the gg channel has to be included at α 4 s since it contributes an equal amount to the total uncertainty.
We finally show the fully matched results in comparison with the ATLAS measurement in fig. 5. The top plot shows the ratio of the ATLAS measurement to our highest-order prediction as in fig. 4. We furthermore included a prediction where the hard scale is chosen as q γ T as suggested in ref. [11]. This more natural scale choice closes the uncertainty gap, and prediction and measurement have now overlapping uncertainty bands.
The bottom plot includes additionally a prediction with the hybrid-cone isolation using an inner-cone radius of R s = 0.1, where the previous agreement at small q T is now destroyed. As already shown in ref. [11], the fixedorder predictions with hybrid-cone isolation at small Φ * are worse than the smooth-cone isolation results. 3 The authors suggest that "The regions where agreement is notably worse are those in the neighborhood of the Sudakov singularities [...], and hence where poor agreement is expected in the absence of resummation". While this is true, the distortion due to the hybrid-cone isolation of the q T distribution cannot be cured by the present q T resummation, but would require some other resummation.
In the limit of inner cone approaching outer cone R s → R o the smooth-cone isolation is restored by definition. In this limit the resummed results agree well with the data. For successively smaller R s the cross-section coming from fixed-cone isolation grows. Since it is always larger than the contribution from smooth-cone isolation, the agreement at large q T also grows. But the increase in cross-section is unfortunately not just at large q T : the smaller R s is taken, the larger the ridge at q T = E iso T becomes (see e.g. fig. 8 in ref. [11]). This ridge effect can be smoothened out to some extend by choosing a Figure 5: Fully matched q T spectra at N 3 LL +NNLO in comparison with the ATLAS measurement. Top plot: Ratio to N 3 LL with smooth-cone isolation and µ R = m γγ in comparison with data and µ R = q γ T . Bottom plot: Similar, but in comparison with prediction using hybrid-cone isolation and µ R = q γ T .
non-constant E iso T , but would also have to be matched by the experimental definition.
We conclude that with present theory frameworks the hybrid-cone isolation is not the answer to a better modeling of photon isolation, especially with increased experimental precision. Distortion effects in the q T distribution due to hybrid-cone isolation are only exacerbated by the q T resummation. Better agreement at large q T is traded with drastic disagreement at small q T , where excellent agreement with resummation is achieved using the smooth-cone isolation. While it is possible to shift the Sudakov singularity in phase space, we believe that at this point the program for fragmentation functions will have to be revived.
For practical comparison with data, the lesson to be learned is likely to just take the more natural scale choice µ R = q γ T , which brings theory and data into better agreement, and include the NNLO γγ+jet corrections [32]. We also show comparison plots for the Φ * distribution in fig. 6 with similar observations at small Φ * , since it is directly correlated to small q T , but with worse agreement at large Φ * . Figure 6: Fully matched φ * spectra at N 3 LL +NNLO in comparison with the ATLAS measurement. Top plot: Ratio to N 3 LL with smooth-cone isolation and µ R = m γγ in comparison with data and µ R = q γ T . Bottom plot: Similar, but in comparison with prediction using hybrid-cone isolation and µ R = q γ T .
Conclusions. We have upgraded previous diphoton predictions at small q T accurate at the level of α 2 s in improved perturbation theory to include the α 3 s "constant" pieces. These include the recently published three-loop qq hard function [19] and the previously implemented two-loop gg hard function [17] together with the three-loop transverse momentum dependent beam functions [33][34][35]. This constitutes an overall primed resummation accuracy of N 3 LL + NNLO 0 . The resummation itself of the qq channel is noticeably stabilized with remaining uncertainties of a few percent at intermediate q T between 10 GeV and 50 GeV . But this is diminished by the large uncertainties from the gg channel at α 3 s and large matching corrections.
With that, we have eliminated the α 3 s qq hard part as a source of uncertainty, and can limit the remaining sources of higher-order uncertainty and contributions for the q T distribution: The dominating uncertainties will be reduced by a matching to γγ+jet at NNLO, which has 10% effects due to the qq channel below 100 GeV [32], and by α 4 s corrections to the gg loop-induced channel, which are mostly relevant for q T 200 GeV.
We have shown that the more natural hard and renormalization scale q γ T with smooth-cone isolation alone restores agreement with data for the q T distribution and for not too large Φ * . On the other hand the hybrid-cone isolation (with fixed E iso T ) introduces a discontinuity in the q T distribution that destroys the otherwise agreement with data at small q T . This is because the isolation is only a power-suppressed effect in the resummation. The leading-power resummation acts at the level of the Born-topology without isolation effects and cannot compensate the Sudakov singularity induced by the hybrid-cone isolation at fixed-order [11]. The upgrades presented in this paper will be included in the upcoming release of CuTe-MCFM 4 .
In the future we plan to match with a γγ+jet calculation at NNLO to take into account the full α 3 s matching corrections at as small q T as possible. If small enough q T can be achieved numerically, one could also extract fixed-order N 3 LO cross-sections using q T slicing, but this will require cutoffs lower than 1 GeV , which is already numerically quite expensive at NNLO [15]. The implementation of the three-loop beam functions also facilitates an implementation of a single-boson N 3 LL + NNLO 1 implementation in MCFM.
Neumann is supported by the United States Department of Energy under Grant Contract DE-SC0012704. This work was supported by resources provided by the Scientific Data and Computing Center (SDCC), a component of the Computational Science Initiative (CSI) at Brookhaven National Laboratory (BNL). [5] ATLAS collaboration, Measurement of the production cross section of pairs of isolated photons in pp collisions at 13 TeV with the ATLAS detector, 2107.09330.