NNLO QCD corrections to three-photon production at the LHC

: We compute the NNLO QCD corrections to three-photon production at the LHC. This is the (cid:12)rst NNLO QCD calculation for a 2 ! 3 process. Our calculation is exact, except for the scale-independent part of the two-loop (cid:12)nite remainder which is included in the leading color approximation. We estimate the size of the missing two-loop corrections and (cid:12)nd them to be phenomenologically negligible. We compare our predictions with available 8 TeV measurement from the ATLAS collaboration. We (cid:12)nd that the inclusion of the NNLO corrections eliminates the existing signi(cid:12)cant discrepancy with respect to NLO QCD predictions, paving the way for precision phenomenology in this process. plays in the apparent perturbative convergence of these two processes. To quantify this in (cid:12)gure 6 we show the composition of the (cid:12)ducial cross-section at LO, NLO and NNLO organized by initial-state partonic reactions. We show the results for three di(cid:11)erent H T -based central scales; the results for the corresponding M T -based scales are very similar.


Introduction
Over the last decade, Next-to-Next-to Leading Order (NNLO) QCD calculations for hadron collider processes have sustained tremendous progress. Owing to the development of many independent approaches  almost all non-loop-induced 2 → 1 and 2 → 2 processes have now been computed, typically in more than one computational approach. Such massive theoretical progress has led to the creation of public codes and has started to produce valuable and solid LHC phenomenology on a massive scale.
The computation of 2 → 3 hadron collider processes represents a natural step beyond the current state of the art in NNLO calculations. Since many of the available computational approaches are generic, in principle they should be able to handle the problem of double real radiation in 2 → 3 processes. The calculation of the so-called real-virtual correction to such processes should also be possible, in principle, due to the existence of numerically stable libraries for one-loop amplitudes. The only ingredient for such calculations which is not readily available are the two-loop five-point amplitudes. Thanks to the development of various new methods [32][33][34][35][36][37][38][39][40][41][42][43][44][45], first results for selected helicities, color structures or kinematics (typically Euclidean) have started to appear. This includes 5-point amplitudes computed in QCD [46][47][48][49][50][51][52][53][54][55][56], in pure Yang-Mills [57][58][59], in supersymmetric theories [60][61][62][63] and in gravity [64][65][66]. In this work we calculate the two-loop amplitude qq → γγγ which is the first time a 5-point two-loop QCD amplitude is derived explicitly, in analytic form, in the physical region. We discuss this result at length in section 2.2 below.

JHEP02(2020)057
The goal of the present work is to demonstrate the feasibility of existing calculational technologies to deal with 2 → 3 hadron collider processes. We have decided to apply this first-ever NNLO 2 → 3 calculation to the process pp → γγγ + X. Our motivation for choosing this process is twofold. First, the production of colorless final states has always occupied a special place among hadron collider processes and has been the pioneering work for both 2 → 1 and 2 → 2 processes. In addition, the calculation of the two-loop amplitudes is more feasible due to the smallest number of scales involved.
Second, the process pp → γγγ + X is of direct phenomenological interest. The crosssection for three isolated photons at the LHC 8 TeV was measured in detail by the ATLAS collaboration [67] (see also the earlier measurement [68]) and was found to be significantly above the NLO QCD prediction in a wide kinematic region. Since at NLO the theory error is completely dominated by missing higher-order terms this process represents a prime case for an NNLO QCD calculation. Indeed, we find that with the help of our calculation this discrepancy can be addressed (see section 3).
The paper is organized as follows. In section 2 we explain our calculation with emphasis on the derivation of the two-loop amplitude. In section 3 we present our predictions for the fiducial and differential cross-section. In section 3.4 we discuss the important question of perturbative convergence in this process. Our conclusions are summarized in section 4.

The calculation
In this work we follow the STRIPPER approach [11][12][13] previously applied at NNLO in QCD to top-pair [69][70][71][72][73] and inclusive jet [74] production at the LHC. The framework is implemented in a fully-differential partonic Monte Carlo program which can calculate any infrared-safe partonic observable. The technical details about our implementation can be found in ref. [74].
The complete calculation converges very well in terms of phase-space integration. Not counting the CPU time needed to evaluate the two-loop finite remainder (see section 2.2.2 for details), it took only about 2k CPU hours to complete. The slowest contribution (about 1k CPU hours) is the real-virtual finite contribution due to the slow evaluation of the 6point one-loop amplitude with OpenLoops 1. That contribution, however, converges fast in terms of required phase-space points.
The ingredients needed for the present calculation are tree-level amplitudes as well as the finite remainders of one-loop and two-loop amplitudes. Their calculation is described in detail in section 2.1 and section 2.2. Here we only point out that all required one-loop amplitudes are included exactly, with full color dependence. The finite remainder of the two-loop amplitude qq → γγγ is included in the leading color approximation, additionally excluding diagrams with closed fermion loops. The justification for this approximation is given in section 2.2 below.
The infrared subtraction operator (sometimes called "Z"-operator) is given in ref. [13]; its leading color approximation can be found in ref. [51]. We work in a theory with 5 massless active quark flavors and renormalize the amplitudes accordingly. No loops with massive fermions are included in our calculation. Their effect in the context of diphoton production has been discussed in ref. [75].

Tree-level and one-loop amplitudes
All tree-level diagrams are computed with the help of the library avhlib [76,77]. For the derivation of the two-loop finite remainder the one-loop amplitude qq → γγγ is needed to order ε 2 (where d = 4−2ε is the space-time dimension). We have computed it following the standard Feynman diagram plus Integration-by-Parts (IBP) identities [78,79] approach. All required master integrals expanded to that order in ε are available in electronic form in ref. [53]. The finite remainders for all one-loop amplitudes are obtained from the library OpenLoops [80,81], while the one-loop squared qq → γγγ contribution is taken from the library Recola [82].
Unlike the case of diphoton production, the gluon-initiated one-loop amplitude gg → γγγ vanishes and thus does not contribute to the process studied in this paper. Since the gg-flux is sizable, the vanishing of this contribution is of phenomenological significance and we discuss it in more detail in section 3.

The two-loop amplitude for qq → γγγ
An important novelty in this work is the calculation of the two-loop amplitude for the process qq → γγγ. Although our calculation is restricted to the leading color approximation, this is the first time a two-loop five-point amplitude is put in a form that can be used in a phenomenological application. For this reason we describe it in detail in this section.

Structure of the two-loop amplitude
We need the two-loop amplitude |M (2) (qq → γγγ) , multiplied by the Born one |M (0) (qq → γγγ) , and summed/averaged over helicities and color. Its color decomposition reads: where N c = 3 is the number of colors. In this work we simplify the calculation by utilizing the following approximation: i.e. we neglect the non-planar contribution M (np) as well as all contributions M (f) with a fermion loop (both planar and non-planar). The non-planar contribution M (np) is suppressed by a factor of 1/N 2 c relative to the leading color one. It is thus expected to be numerically subdominant. The non-planar contribution cannot be computed at present since the required IBP solutions (topologies B 1 and B 2 in the notation of ref. [38]) are not yet fully known.
The contribution M (f) contains all diagrams with one closed fermion loop. Both planar and non-planar diagrams contribute to it. It cannot be currently derived since the required non-planar IBPs (topology B 2 from ref. [38]) are not yet known. The term M (f) is suppressed with respect to the leading color one by a single power 1/N c . At the same time JHEP02(2020)057 some diagrams 1 are enhanced by the number of massless fermion flavors n f = 5. Therefore, although in the strict N c → ∞ limit M (f) is suppressed relative to the leading color contribution, its numerical value may not necessarily be sub-dominant with respect to eq. (2.2). For this reason, to be conservative, one should assume that it is comparable numerically to the leading color contribution. As we show in section 3 below (see in particular figure 6), the numerical impact of the leading color approximation eq. (2.2) to the differential crosssection is itself negligible, at the percent level, which aposteriori justifies the approximation M (f) ≈ 0. In the future, once the corresponding contribution M (f) becomes known (same for the term M (np) ), we can easily update our cross-section predictions.

Calculation of the two-loop amplitude
To compute the two-loop amplitude we use a standard Feynman diagram-based approach. The diagrams are generated with the help of a private software. After multiplying with the Born amplitude and then computing the traces of spin tensors and color factors, the resulting scalar integrals are mapped to master integrals using the IBP results of ref. [38].
The last step is the inserting of the results for the required master integrals. To that end we utilize the results of ref. [39] where a set of integrals has been explicitly computed in terms of the so-called pentagon functions f ij . This set of integrals can be algebraically related to the set of master integrals in ref. [38] with the help of the IBP solution derived in that latter reference.
At this point the bare amplitude can be computed numerically using the routines for the numerical evaluation of the set of integrals provided with ref. [39]. We do not follow this approach here for two reasons. First, we would like to provide an explicit analytic result in terms of basic functions, like the set of pentagon functions. Second, the complete results involves not just the 61 master integrals but also many integrals that are obtained from them by crossings of external legs. In practice we have more than 70 set of crossings that need to be applied to the set of master integrals. While not every master integral will need to be crossed for all crossings, the complete set of integrals, accounting for all crossings, far exceeds the dimension of functions needed to describe the amplitude. For this reason such an approach would not be minimal and could lead to more severe loss of precision during the numerical evaluation.
For the above reasons, we use the explicit representation of master integrals in terms of pentagon functions f ij [39] and have applied the momentum crossings directly to those functions. To minimize the set of functions, we have derived various functional identities between those functions with different arguments. As a result, we have derived an explicit expression for the squared amplitude eq. (2.2) as a polynomial in transcendental constants and f ij functions with various arguments. The coefficients of this polynomial are rational functions of the kinematic invariants. They have been factorized and simplified, in some cases using the finite-field reconstruction package FiniteFlow [83].
Besides the usual ζ(2) and ζ(3) a new set of constants, collectively called bc4, appear at weight 4 [39]. Their treatments requires special attention. These constants are associated JHEP02(2020)057 with the master integrals at weight 4 and, despite being called "constants", in general take different numerical values in the various physical regions. We have accounted for this possibility in the process of applying momentum crossings to the master integrals. Many of these constants take the same value in the various physical regions. We have utilized the numerical values which are included in the numerical code accompanying ref. [39]. Similarly, the analytic continuation of the pentagon functions f ij across the various physical regions is performed automatically by the numerical library of ref. [39]. To check the correctness of our manipulations, we have compared in each physical region the numerical predictions for each master integral, constructed by us as described above, with the numerical value for the master returned directly by the library of ref. [39] and have found complete agreement. We have also checked many integrals against the numerical program pySecDec [84], finding agreement in all cases.
In summary, we have expressed the complete analytic result for the bare amplitude in a basis of about 1800 transcendental terms involving ζ(2), ζ(3) and f ij functions plus about 100 terms involving bc4 "constants" of weight 4.
Most of the rational coefficients are small (i.e. kB size) but some exceed 1MB. The loss of numerical precision due to cancellation between the various terms is thus of particular concern. To minimize such cancellation we have evaluated all rational coefficients with exact rational arithmetic. Specifically, we rationalize each phase-space point by preserving the accuracy of the original floating point number, and then use its rational form to compute the rational coefficients as exact rational numbers. This is implemented with the help of the CLN library [85]. The evaluation is much slower than the evaluation in double precision, however the overall timing is negligible compared to the evaluation of the slowest pentagon functions of weight 4. We have performed various tests for the depth of numerical cancellations and have found them to be under control in all test cases.
The numerical evaluation of the functions f ij is performed with the help of the C++ library provided with ref. [39]. The time it takes to evaluate these functions depends strongly on their weight. All functions through weight 3 are standard polylogarithms and can easily be computed with full double precision in negligible time. The functions of weight 4 are the slowest and can take up to several minutes per phase-space point. Their precision is less than full double precision due to conflicting requirements of precision and speed as well as the numerical stability of the integration routines used for their calculation. With the help of extensive experimentation we have found that computing them with at least 7 significant digits is sufficient for our purposes. To test the depth of numerical cancellations we have also computed the weight 4 functions requiring 5 significant digits. This results in a finite remainder with, typically, 2 significant digits.
It takes about 10-50 min, depending on the phase space point, to compute the finite remainder in a single phase-space point on a single CPU (i.e. without any parallelization). The average time is about 17 min when a relative precision of 10 −7 for the weight 4 functions is requested. In general, several hundred thousand points are required in order to integrate the three-photon phase space over the required bins. Such a calculation requires significant, cluster-size computer resources. While a one-off evaluation is possible it poses non-trivial problems, especially if re-evaluation of the amplitude is needed (for example for a different JHEP02(2020)057 setup or collider energy). To minimize this computational effort we have utilized a two-fold strategy.
First, we have produced an optimized set of phase-space points which have been generated according to the Born cross-section. Such an approach has already been used in refs. [86][87][88][89][90] and it allows us to obtain a good quality double-virtual contribution with a reduced number of events. In practice, we have computed 30k events.
Second, we have utilized an approach to the implementation of two-loop finite remainder where the above mentioned 30k phase space points have been used to construct a (four-dimensional) interpolating function for the real part of the finite remainder. Constructing multi-dimensional interpolating functions is a hard problem. In our case we have used the purposely developed library GPTree [91] which uses advances in machine learning to optimize the interpolation tables and to produce an estimate of the interpolation error. The output of the GPTree library is a C++ library which is portable and very easy to link to a C++ code and to use. It has the advantage that if more phase-space points are computed in the future the interpolation tables can be refitted and thus be further improved. This library will be made public in a future publication. We have also found it very useful as an additional monitor for the appearance of numerical instabilities.
In view of the phenomenological application of the two-loop amplitude eq. (2.1), see section 3, we make explicit the scale dependence of its finite reminder H (2) : where s 12 is the partonic c.m. energy squared. Since the coefficients c n can be determined exactly from the tree-level and one-loop amplitudes, throughout this work the scale dependence of the two-loop finite reminder is included with full color dependence. Therefore, the approximation eq. (2.2) is applied only to the first term in the r.h.s. of eq. (2.3). As can be seen in figure 6 below, the numerical impact on the NNLO cross-section of the scale-independent part of the two-loop finite remainder H (2) (s 12 ) is rather small, at the percent level. The explicit expressions for the amplitude, as well as further details about its evaluation, will be given in a subsequent publication.

LHC setup
Our calculational setup follows the 8 TeV ATLAS measurement [67]. The definition of histograms and experimental data is taken from the corresponding HEPData entry [92]. Our event selection is based on the following phase-space cuts: • E T cut for the three photons: E T,γ 1 > 27 GeV, E T,γ 2 > 22 GeV and E T,γ 3 > 15 GeV, where γ 1 represents the hardest photon while γ 3 the softest one.

JHEP02(2020)057
• Photon separation: the angular distance ∆R between any two photons is required to • A minimum three-photon invariant mass is required: m γγγ > 50 GeV.
• Following Frixione [93], we impose smooth photon isolation. Specifically, for any angular distance ∆R away from each photon, subject to ∆R ≤ ∆R 0 , we require where ∆R 0 = 0.4 and E max T = 10 GeV. The energy E iso T (∆R) is defined as The sum in the above equation is over all final-state partons i, and E T,i and ∆R i,γ are parton i's transverse energy and angular distance with respect to the photon.
Our calculation uses the NNPDF31 nnlo as 0118 pdf set [94]. We have not computed the pdf error; it was estimated in ref. [67] and found to be below the (NLO) scale variation.
In this work we have utilized two different forms for the dynamic factorization and renormalization scales: Our default central scale choice is µ 0 = H T /4, which follows from the findings of ref. [72]. In fact, in the following we have studied the choices µ 0 = H T /n, with n = 1, 2, 4 as well as the alternative choices µ 0 = M T /n, with n = 1, 2, 4, that are based on the transverse mass of the three-photon system. The M T -based scale was used in the latest diphoton production study [95]. We find that the differences between calculations with central scales M T /n and H T /n, with n = 1, 2, 4, are relatively small. Scale variation of the factorization and renormalization scales is derived with a standard seven-point variation around the central scale µ 0 .

Fiducial cross-section
In figure 1 we compare ATLAS data [67] with the predictions for the fiducial cross-section as defined in section 3.1. We compare predictions based on 6 different renormalisation and factorization scales, in LO, NLO and NNLO QCD. In all cases we observe large shifts from LO to NLO and from NLO to NNLO which are much larger than the scale variations at, respectively, LO and NLO. Specifically, for our default scale µ 0 = H T /4 we have an NLO/LO correction of about 2.8 while the NNLO/NLO correction is about 1.6. We discuss this important feature in section 3.4 below. The scales µ 0 = H T /4 and µ 0 = H T /2 both agree with data, especially the H T /4 one. The scale µ 0 = H T is only just outside the measurement's uncertainty band. For simplicity in this work we did not compute the full scale variation around the scale H T (same for M T ) which is why only the central value is shown for these two scales. We do not expect the scale variation around the two scales will be much different that the pattern already emerging from figure 1.
In general we observe that the scale variation increases when going from LO to NNLO and that all scales are consistent at a given order within their scale uncertainties. For a proper interpretation of the reliability of the theoretical predictions it is therefore imperative to understand the issue of perturbative convergence. We devote section 3.4 to this issue but here we only say in advance that we believe the NNLO predictions are probably the first order for which the theory prediction, with its associated scale variation, is reliable.
To summarize, based on the above discussion we conclude that our default scale choice is in perfect agreement with the experimental measurement

Differential distributions
A very large number of differential distributions have been measured by the ATLAS collaboration in ref. [67]. In this work we have computed the theory predictions in NNLO QCD for all of them. We start by showing in figure 2 the predictions for the p T distributions of the three individual photons: the hardest one γ 1 (left), γ 2 (center) and the softest one γ 3 (right). The plots for all other differential distributions shown next follow the same pattern: in figure 3 we show the three ∆Φ angles between the three pairs of photons, in figure 4 we show the three rapidity differences |∆η| between the three pairs of photons and, finally, in figure 5 we show the invariant mass distributions between the three pairs of photons as well as the invariant mass of the three-photon system.
Overall, a very consistent picture arises from all differential distributions, both in relation to the properties of the theory predictions as well as in relation to their agreement with data.
The most notable feature evident in all differential distributions are the large jumps from LO to NLO and from NLO to NNLO. The difference between orders is much larger than the corresponding scale variations at LO and NLO which, in principle, raises the question of the validity of perturbative convergence in this process. This behavior closely JHEP02(2020)057    resembles the one already discussed for the fiducial cross-section. At this point we will only mention that we believe the NNLO QCD predictions is likely already a reliable prediction which can be confidently compared to data. We leave the detailed discussion of this point to section 3.4.
The second notable feature is the overall good agreement between NNLO QCD predictions based on a scale H T /4 with data. While in most distributions there are bins that do not agree with the NNLO prediction, the overall shape as well as normalization of all distributions is clearly correctly described at NNLO. In fact in some of the bins where deviations is observed could be due to larger statistical fluctuations in data. An improved future measurement will clearly be very useful to clarify this. Interestingly, it is the distributions ∆Φ and ∆η that are described best and in fact in those two we observe perfect agreement between NNLO and data for all pairs of photons.
The relative MC error on the differential NNLO predictions shown here is below 3 percent. The theoretical predictions for all distributions, based on our default scale H T /4, are available in electronic form as supplementary material attached to this paper.
In summary, we would like to stress that in this calculation we have only accounted for the QCD corrections through NNLO. Other theoretical contributions should at this point also be revisited. These include electroweak corrections and (refining the study of) JHEP02(2020)057 effects due to photon isolation. Effects due to pdfs appear to be subdominant to the scale variation at NNLO but this should also be cross-checked in a more complete study. The issue of the "best" scale choice can always be debated and at this level of precision seems to be a dominant source of theoretical uncertainty. For a detailed phenomenological study the MC error of the predictions shown here can be improved further. Finally, for completeness, one would like to have the complete NNLO prediction by including the contributions to the two-loop finite remainder neglected in this work, although we expect them to be phenomenologically insignificant. JHEP02(2020)057

Discussion of perturbative convergence
As in diphoton production [75,95,96], the inclusive production of three photons exhibits behavior that at first glance is inconsistent with perturbative convergence. Indeed, as emphasized above, we observe large jumps from LO to NLO and from NLO to NNLO. These jumps are much larger than the corresponding scale variation bands at LO and NLO. This behavior is evident in all differential distributions as well as in the fiducial cross-section. Specifically, we recall that for our default scale choice the fiducial cross sections at NLO exceed the LO one by a factor of 2.8 while the NNLO/NLO K-factor is about 1.6. This behavior is very similar to the one encountered in diphoton production. Various arguments have been given in the past for the appearance of such large Kfactors in diphoton production. Two of those arguments are the presence of asymmetric cuts imposed on the two photons as well as the sizable loop-induced gg → γγ contribution. While these arguments have their merit, it is easy to see that they are not the drivers behind the behavior we are trying to understand in diphoton (as well as three-photon) production. For example, the asymmetric cuts should not play an appreciable role for three-photon production because the Born state is naturally asymmetric. Similarly, while the loop-induced reaction is very large relative to the LO diphoton cross-section, its relative contribution at NNLO is not that sizable, only of the order of 10% [95]. While such a contribution is important it is not large enough to be the driver behind the large K-factors observed in both processes. In fact, this issue can be cleanly understood in three-photon production process where the corresponding loop-induced amplitude gg → γγγ vanishes.
The above analysis of the gg-driven correction brings a very important point, namely, the role the initial-state flux plays in the apparent perturbative convergence of these two processes. To quantify this in figure 6 we show the composition of the fiducial cross-section at LO, NLO and NNLO organized by initial-state partonic reactions. We show the results for three different H T -based central scales; the results for the corresponding M T -based scales are very similar.
What we observe in figure 6 is very illuminating. First, we note that the gg flux does contribute (due to double real emissions and collinear subtractions) although its effects is marginal, in the few percent range, depending on the choice of scale. Clearly, despite the fact the gg flux is very large it nevertheless has negligible effect on the cross-section simply because the corresponding partonic cross-sections are very small. The large gluon pdf does have a substantial impact on the three-photon cross-section but this happens through the qg reaction. As also emphasized in ref. [95] for the case of diphoton production, the qg reaction starts to contribute only at NLO. This leads to a very unique interplay between purely partonic contributions, including their radiative corrections, and partonic fluxes. Specifically, the qq contribution receives a sizable but not huge NLO radiative correction. At NLO this contribution is now dwarfed by the newly generated qg correction which at this point is only LO. At NNLO the qq result gets another significant yet moderate correction, and the qg reaction also receives sizable but reasonable perturbative correction. At NNLO two new channels open -the gg and qq ones, the latter being much more significant than the former. As also concluded in ref. [95] for the case of diphoton production, NNLO is the first order where all large partonic reactions have already been included together with higher-order corrections to the largest ones; therefore one can reasonably expect that from this point on the yet-higher order N 3 LO corrections to be derived at some point in the future are likely to start showing a more convergent behavior.
Before closing this section we would like to emphasize that the pattern of scale dependence observed when going from LO through NNLO should not be viewed as anomalous. The fact that scale dependence increases towards NNLO is simply due to the fact that the scale variation at the lower orders is artificially small and that, as explained in this section, at each new order through NNLO new large partonic reactions enter the process thus increasing the overall scale dependance. The arguments given here imply that starting at N 3 LO the scale variation should start to decrease. This will be very interesting to check in the future. In summary, in our view, the above arguments imply that the scale dependence of the NNLO prediction is likely not artificially small.

Conclusions
In this work we calculate the NNLO QCD correction to three-photon production at the LHC. Our calculation is complete except for the scale-independent part of the two-loop finite remainder which is included in the leading color approximation. We estimate the effect of the missing two-loop contributions. We expect they are phenomenologically insignificant.
Our calculation is the first NNLO calculation for a 2 → 3 scattering process. Although the production of colorless final states is not as complicated as a generic 2 → 3 reaction, we believe that our calculation clearly demonstrates that current computational technology is capable of dealing with the complicated structure of infrared singularities in multi-final state processes. In particular, based on our experience with dijet production (which was JHEP02(2020)057 computed within the same STRIPPER framework as the present calculation) we think that the NNLO computation of three jets at the LHC is feasible.
An important part of our work is the calculation of the two-loop amplitude qq → γγγ in the leading color approximation. We have expressed the amplitude in a fully analytic form, defined directly in the physical region. To the best of our knowledge this is the first time a 5-point two-loop amplitude has been expressed in a form readily available for phenomenological applications. To this end we had to go beyond simply producing an analytic result that implements the many momentum crossings inherent to processes with many particles in the final state together with the relevant analytic continuations; we have extensively investigated the question of numerical stability and have been able to evaluate the amplitude numerically in about 30k phase-space points with sufficient numerical precision. We would like to stress that this problem is highly non-trivial due to the large size of the amplitude and the large number of independent transcendental functions that appear in it.
The evaluation of the two-loop amplitude is expensive in terms of CPU time. We have investigated two possibilities to mitigate this problem: one involves specially generated phase-space points that accelerate the convergence of the phase-space integration, while the other involves the construction of a four-dimensional interpolating function which internally utilizes machine-learning techniques. We find that both approaches lead to compatible predictions within the corresponding Monte Carlo errors.
We observe that the structure of higher-order corrections in inclusive three-photon production is very interesting and resembles closely the one known from diphoton production. We find very large higher-order corrections; the NLO prediction for the fiducial cross-section is larger than the LO one by a factor of 2.8 while the NNLO exceeds the NLO one by a factor of 1.6. We have presented detailed analysis of the anatomy of the higher-order corrections in this process and have concluded that the NNLO prediction is likely to be reliable.
Finally, we have compared our predictions with the high-quality LHC data available from the ATLAS Collaboration. We find that the sometimes huge discrepancies between QCD predictions and data noted previously at NLO are absent at NNLO and that the NNLO prediction agrees well with data for all distributions. This result clearly demonstrates how indispensable higher-order corrections are to quantitative phenomenological LHC analyses.

JHEP02(2020)057
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.