Triphoton production at hadron colliders in NNLO QCD

We present next-to-next-to-leading-order (NNLO) QCD corrections to the production of three isolated photons in hadronic collisions at the fully differential level. We employ qT subtraction within MATRIX and an efficient implementation of analytic two-loop amplitudes to achieve the first on-the-fly calculation for this process at NNLO accuracy. Numerical results are presented for proton-proton collisions at energies ranging from 7 TeV to 100 TeV. We find full agreement with the 8 TeV results of arXiv:1911.00479 and confirm that NNLO corrections are indispensable to describe ATLAS 8 TeV data. We demonstrate the significance of NNLO corrections for future precision studies of triphoton production at higher collision energies.

is important to constrain anomalous Higgs couplings in rare Higgs boson decays [39][40][41] or in the rare Higgs boson production process in association with a photon [42] with the Higgs boson decaying into a pair of photons. Moreover, triphoton production is relevant as a background to the associated production of a photon with a beyond-the-SM (BSM) particle that decays into a photon pair, see Refs. [43][44][45] for instance.
On the theoretical side 2 → 3 reactions are the current edge for NNLO QCD calculations, limited mostly by the complicated computation of two-loop corrections to five-point functions. However, for the progression of precision phenomenology at the LHC it is indispensable to go beyond the current state-of-the-art for 2 → 3 processes, which is next-to-leading order (NLO) QCD accuracy. Triphoton production is the only 2 → 3 process for which NNLO corrections have been calculated [1]. The NLO cross section for the production of three isolated photons had been calculated already some time ago [46] using smooth-cone isolation, and also considering the fragmentation contribution at leading order (LO) in Ref. [47].
In principle, there are two mechanisms that are relevant for the production of isolated photons: the direct production in the hard process, which can be described perturbatively, and the production through fragmentation of a quark or a gluon, which is non-perturbative. Since the latter production mechanism relies on the experimental determination of fragmentation functions with relatively large uncertainties, we will exploit smooth-cone isolation as suggested by Frixione in Ref. [48] to completely remove the fragmentation component in an IR-safe manner. This substantially simplifies theoretical calculations of processes with isolated photons beyond the LO. Although the finite granularity of the calorimeter prevents the usage of smooth-cone isolation in experimental analyses so that they have to rely on an isolation with a fixed cone, using certain smooth-cone parameters mimics fixed-cone isolation criteria, see e.g. Ref. [49].
In this letter, we present a new calculation of the fully differential NNLO cross section for the production of three isolated photons. We exploit the Matrix 1 framework [50], using its fully general implementation of the q T -subtraction formalism [51], and demonstrate that this NNLO approach is suitable to cope with 2 → 3 colour singlet processes. We achieve an efficient implementation of the two-loop helicity amplitudes for qq → γγγ based on the analytic calculation of Ref. [52]. This is the first time that a five-point two-loop amplitude is presented which is sufficiently fast to be calculated directly in the physical region during phase space integration. Therefore, at variance with Ref. [1] our calculation allows us to obtain NNLO corrections on the fly. Its implementation in the Matrix framework for the first time enables a fully flexible calculation of a 2 → 3 process at NNLO accuracy. 2 We study numerical results for fiducial cross sections and distributions in proton-proton collisions for centre-of-mass energies ranging from 7 TeV to 100 TeV.
We consider the production of three isolated photons, i.e. the process where X indicates inclusiveness over additional radiation. At LO, the triphoton production cross section is of O(α 3 ), where α is the electroweak (EW) coupling, and it is qq initiated as shown in Figure 1 (a). At variance with diphoton production, the loop-induced gg-initiated contribution shown in Figure 1 (b), which would enter at O(α 3 α 2 S ), with α S being the strong coupling, vanishes in 1 Matrix is the abbreviation of Munich Automates qT subtraction and Resummation to Integrate X-sections by M. Grazzini, S. Kallweit, M. Wiesemann. The program is available under http://matrix.hepforge.org. 2 The implementation of NNLO corrections to triphoton production will be made publicly available with the next release of Matrix. A preliminary version of the code is available from the authors upon request. the case of triphoton production due to the charge-conjugation symmetry of QED⊗QCD. The first non-vanishing contributions of loop-induced type, see Figure 1 (c) and (d), are gg-or qg-initiated, and they enter at O(α 3 α 3 S ), i.e. beyond the nominal accuracy of our calculation. Nevertheless, since those contributions are separately finite and represent the first non-vanishing order of the loop-induced production mode, we have calculated them and found their effect to be below 1% of the integrated NNLO cross section, and only slightly larger when considered differentially in phase space. Thus, we can safely neglect the loop-induced contributions in what follows.
For our calculation we employ the Matrix framework [50]. All tree-level and one-loop amplitudes are evaluated with OpenLoops [53][54][55]. At the two-loop level we have performed a novel implementation based on the analytic expressions in Ref. [52], which is highly efficient and evaluates the 2-loop hard function within few seconds for each phase space point, as discussed in more detail below. NNLO accuracy is achieved by a fully general implementation of the q T -subtraction formalism [51] within Matrix. The NLO parts therein (for γγγ and γγγ+1-jet) are calculated by Munich 3 , which uses the Catani-Seymour dipole subtraction method [58,59]. The Matrix framework features NNLO QCD corrections to a large number of colour singlet processes at hadron colliders. 4 It has already been used to obtain several state-of-the-art NNLO QCD predictions [19, 20, 23, 24, 26-30, 32, 33] 5 , and for massive diboson processes it was extended to combine NNLO QCD with NLO EW corrections [68] and with NLO QCD corrections to the loop-induced gluon fusion process [69,70]. Through the recently implemented interface [71,72] to the code RadISH [73,74] this framework now also includes the resummation of transverse observables such as the transverse momentum of the produced colour singlet final state.
The q T -subtraction formalism is employed for the very first time for a colour singlet process of the given complexity. Not only is triphoton production a 2 → 3 process, it also involves the isolation of all of the three photons. Already processes with a single isolated photon, such as Zγ production, feature rather large power corrections in the transverse momentum of the colour singlet system, as shown in Ref. [50]. Due to the interplay between the smooth-cone isolation criteria and the slicing cutoff (introduced as r cut in Ref. [50]) of the q T -subtraction approach, the production of three isolated photons could be subject to relatively large systematic uncertainties at NNLO. To deal pp → γγγ @ 8 TeV, µ 0 = H T /4 Figure 2: Dependence of the NNLO cross section for pp → γγγ + X on the slicing parameter r cut (red points with numerical error bars), the extrapolated cross section for r cut → 0 (orange, solid) and comparison to the results from Ref. [1] (blue, dashed).
with this issue, the r cut slicing parameter is fully monitored and controlled by Matrix, including a completely automated cutoff extrapolation r cut → 0 performed with every run [50]. Figure 2 shows the NNLO cross section for triphoton production, within the fiducial cuts specified in Table 1, as a function of r cut (which denotes a lower cut on the dimensionless quantity r = p T,γγγ /m γγγ ) and its extrapolation to r cut = 0 with an estimate of the respective uncertainties. We find that q T subtraction is fully capable of dealing with 2 → 3 colour singlet processes, even when three isolated photons are involved, and that we can control the numerical integration and the systematic uncertainties induced by the r cut dependence of the cross section at the few permille level, which fully suffices for any phenomenological application. 6 Before presenting phenomenological results we comment in more detail on our calculation of the double-virtual contribution. We exploit the compact analytic expressions for the finite remainders of two-loop and one-loop helicity amplitudes reported in Ref. [52]. They were obtained within the C++ framework Caravel [75], based on the numerical unitarity method and analytic reconstruction techniques [76][77][78][79][80]. We use the analytic expressions for the finite remainders to implement the numerical evaluation of the two-loop hard function, defined in Eqs. (12) and (62) of Ref. [81], for the qq → γγγ process within Matrix. 7 The finite remainders are expressed in terms of a judiciously constructed basis of multivariate transcendental functions [82]. Their numerical evaluation relies on the library PentagonFunctions++ provided in Ref. [82]. The evaluation of the two-loop hard function on average takes only a few seconds per phase space point, and the numerical error is far below the error of the phase space integration. Thus, the two-loop amplitude can be evaluated 6 For the differential distributions presented in this letter we have performed a bin-wise extrapolation r cut → 0 as used for instance in Refs. [30,61] and found extrapolation effects to depend rather mildly on the region of phase space. Therefore, rescaling distributions computed at a fixed r cut value with the integrated cross section in the r cut → 0 limit yields a suitable approximation. 7 To account for different definitions of the infrared subtraction operators, the necessary finite shifts are applied.
on the fly during phase space integration, so that there is no need to construct an interpolating function from a pre-generated grid in phase space for the two-loop contribution, as it was done in Ref. [1]. This implicates not only a vast improvement in terms of flexibility and applicability of the resulting NNLO calculation, but it also eliminates any interpolation uncertainties from the ensuing results.
The hard function at two-loop level for the process qq → γγγ can be written as where N c is the number of colours, N f the number of light quark flavours, C F = (N 2 c − 1)/(2N c ) the quadratic Casimir operator of the fundamental representation, and Q f the ratio of the electric charge of the quark with flavour f to the electric charge of the quarks in the initial state, and T F = 1/2. The functions H (2,N f ) and H (2,Ñ f ) describe the contributions from two types of Feynman diagrams with closed quark loops. The former captures the diagrams with no photons coupling to the quark loop, and the latter those with two of the photons attached to the quark loop. 8 The functions H (2,1) and H (2,Ñ f ) require the evaluation of non-planar Feynman diagrams that are beyond reach with current computational techniques. Thus, while our predictions include the full colour dependence in all other contributions, we compute H (2) in the approximation where we keep only the leading term N 2 c H (2,0) in the formal N c → ∞ limit. 9 The contribution from H (2,1) is suppressed by a factor of 1/N 2 c and can therefore be safely neglected. In a naive counting, the terms H (2,N f ) and H (2,Ñ f ) are expected to be of similar size as the leading-colour term H (2,0) , but for other processes, such as diphoton [83] or dijet [84] production for instance, those contributions are numerically subleading with respect to the leading-colour term. In order to gauge the quality of our approximation, we consider the directly related diphoton process as a proxy. The two-loop hard function for this process is of the same form as given in Eq. (1), and all subleading corrections are known [83]. For brevity we consider only the dominant subprocess uū → γγ here. Figure 3 shows the impact of each individual correction beyond leading colour. All three corrections are of the order of 10%, adding up to at most 30%. In the case of triphoton production, only the functions H (2,0) and H (2,N f ) are known [52]. Since not all two-loop amplitudes with closed quark loops have been calculated, for consistency we do not include H (2,N f ) in our calculation. Nevertheless, we calculated the H (2,N f ) contribution separately and found that its size is about −18% of the H (2,0) contribution, and largely independent of phase space region and centre-of-mass energy. This is in line with our assumption that N f -related corrections are numerically subleading. We have determined the contribution of H (2) in our approximation to be about 3% for the integrated NNLO cross section at 8 TeV, 10 while it reaches up to 6% in the differential distributions considered in Ref. [37]. We confirm this statement also for high-energy tails that are not resolved in this 8 The contribution from diagrams with one or three photons attached to the quark loop vanishes for the same reason as the amplitude gg → γγγ. 9 The same approximation was used in Ref. [1]. 10 Ref.
[1] quotes a similar relative size of the finite part of the double-virtual corrections in their subtraction scheme for the integrated cross section. They further argue that even in a conservative uncertainty estimate of the approximation at hand for the two-loop amplitude, where one assumes it to be wrong by at most 100%, the ensuing uncertainty is still below the scale uncertainties, which should be sufficient for phenomenological applications. In fact, we reckon that such uncertainty estimate might be even too conservative as we argue in the main text.   1), and the full result. The lower frame shows the relative differences with respect to H (2,0) . measurement, and we find slightly larger effects only for the TeV ranges of some invariant-mass distributions. The overall impact of H (2) on the fiducial cross section continuously decreases with increasing collision energy, by roughly a factor of 2 when going to 100 TeV, and the same trend is observed for differential observables as well. We therefore conclude that the effect of our approximation is negligible for phenomenological applications.
We present predictions for proton-proton collisions at centre-of-mass energies ( √ s) ranging between 7 TeV and 100 TeV. We employ the complex-mass scheme [85] and use the G µ scheme for the EW parameters, i.e. we evaluate the EW coupling as α = √ 2 G F |µ 2 W (1 − µ 2 W /µ 2 Z )| /π and the mixing angle as cos θ 2 while a customary 7-point variation is used to estimate the uncertainties due to missing higher-order corrections. Thus, the renormalization and factorization scales are varied around µ 0 by a factor of two with the constraint 0.5 ≤ µ R /µ F ≤ 2, and the band between the minimum and maximum value of the cross section estimates the uncertainty. The scale setting in Eq. (2) was the preferred one in Ref. [1], which allows us to directly compare to their results.
fiducial setup for pp → γγγ + X; used in the ATLAS 8 TeV analysis of Ref. [37] p   We study integrated cross sections and distributions in the fiducial region. The cuts are summarized in Table 1 and correspond to the fiducial phase space definition of the ATLAS 8 TeV analysis [37]. Those cuts involve different transverse-momentum thresholds for the three photons as well as a requirement on their pseudorapidities. Furthermore, we impose a separation of each pair of photons in ∆R = ∆φ 2 + ∆η 2 and a lower bound on the invariant mass of the three-photon system. Finally, the photons are required to be isolated, which is achieved by means of Frixione's smooth-cone isolation [48] with a fixed value (instead of one relative to the transverse momentum of the respective photon) of the threshold E ref T in the smooth-cone condition as given in Eq. (3) of Ref. [50].
We start by discussing fiducial rates in proton-proton collisions as a function of the machine energy shown in Figure 4. The corresponding numbers and K-factors are quoted in Table 2 at the LHC  energies 7 TeV, 8 TeV, 13 TeV and 14 TeV, and for two potential future colliders, the HE-LHC at 27 TeV and the FCC-hh at 100 TeV. We recall that producing these results is possible due to our fast and efficient calculation of triphoton production within Matrix, which allows us to obtain predictions for any setup and any centre-of-mass energy within a couple of days on a medium-sized cluster. We find full agreement with the results in Ref. [1], which are restricted to 8 TeV, with a quoted cross section of 67.5 +11% −8% fb at 8 TeV. We have also reproduced all differential distributions in the identical setup as in Ref. [1], and we found full agreement with the results in the ancillary files of Ref. [ Table 2: Predictions for fiducial pp → γγγ + X cross sections at different centre-of-mass energies; the numbers in brackets are integration errors, while at NNLO they also include systematic uncertainties from r cut dependence, see Ref. [50]; the percentages correspond to scale uncertainties; were subject to a minor bug, which we could reproduce by treating them as overflow bins, i.e. by including all events beyond the upper bound into the last bin, while keeping its normalization unchanged. 11 The results in Figure 4 and Table 2 show that higher-order corrections to triphoton production are very sizeable and absolutely crucial to obtain a reliable prediction. Already at LHC energies the NLO cross section is about a factor of three larger than the LO one, and the corrections further increase with the collision energy. Even NNLO corrections are larger than 50% at the LHC and reach almost a factor of two at the FCC-hh. Thus, NNLO accuracy is indispensable for an accurate prediction of the triphoton cross section. Indeed, the large discrepancy of several standard deviations observed in the ATLAS 8 TeV measurement of Ref. [37] when compared to NLO-accurate predictions is completely eliminated by the NNLO corrections, with the NNLO cross section being fully consistent with the measured 8 TeV cross section of 72.6 +6.5 −6.5 (stat) +9.2 −9.2 (syst) fb. Due to these sizeable corrections, one may question the validity of scale variations to estimate uncertainties from missing higher orders. Clearly, neither the scale uncertainties at LO nor the ones at NLO capture the corrections of the next order. In addition, contrary to what one usually expects, the uncertainties do not (or hardly) decrease upon inclusion of both NLO and NNLO corrections, especially at lower collider energies, where they even slightly increase in some cases. The fact that scale uncertainties at LO and NLO do not cover higher-order corrections can be largely attributed to a well-known observation that opening of new partonic channels induces sizeable corrections in fixed-order predictions, especially if the opened channels are enhanced by large initial-state flux. 12 Indeed, other processes with similarly large NLO and NNLO corrections, such as single Higgs production, receive a rather small contribution at N 3 LO, well within the uncertainties estimated from NNLO scale variations. This argument is also supported by an observation made in Ref. [88] about the pp → W bb process which suffers from large NLO K-factors. Considering the same process with a large enough number of additional jets in the final state, such that all partonic channels are included already at LO, one observes that the K-factors shrink significantly, and the NLO predictions end up being within the LO scale uncertainties. For triphoton production all partonic initial states are first included at NNLO. Therefore, corrections beyond NNLO are expected to be much smaller and within the estimated uncertainty. Also the significant decrease from relative NLO to NNLO corrections can be seen as a heuristic hint towards perturbative convergence. We note that it would be interesting to investigate the impact of scale choices and photon isolation criteria on the theoretical predictions and the respective uncertainties in the future, as done recently in the case of diphoton production [22]. A more conservative uncertainty estimate that is reliable also at LO and NLO can be obtained by following the approach suggested in Ref. [89]. Unfortunately, this approach has been formulated only for a fixed setting of the factorization and renormalization scales so far and thus cannot be directly applied to our predictions.
We continue our discussion by comparing differential distributions at LO, NLO and NNLO to ATLAS data at 8 TeV [37] in Figures 5-7. We show the invariant-mass spectrum of the triphoton dσ/dp T,ɣ1 [fb/GeV] system (m γγγ ) and of photon pairs (m γ i γ j ) in Figure 5, the transverse-momentum distributions (p T,γ i ) of the three photons in Figure 6, the pseudorapidity differences between two photons (∆η γ i ,γ j with i, j ∈ {1, 2, 3} and i = j) in Figure 7 (left), and the difference in the azimuthal angle between two photons (∆φ γ i ,γ j ) in Figure 7 (right). The agreement between the NNLO predictions and data is truly remarkable. With only few exceptions, in particular the last bin of the transverse-momentum distributions of both the hardest (p T,γ 1 ) and the second-hardest (p T,γ 2 ) photon in Figure 6, all data points agree with the NNLO predictions within one standard deviation. Apart from the pseudorapidity differences, where the NNLO corrections are essentially flat, they induce large effects on the shapes of all other distributions. Indeed, besides the corrected normalization these shape distortions at NNLO are absolutely crucial to describe the measured distributions well. In several cases the NNLO/NLO K-factor becomes as large as a factor of two, and it can even reach a factor of three and more for the ∆φ γ i ,γ j distributions in Figure 7. Assuming no substantial BSM effects, the excellent agreement with available data suggests that NNLO scale bands indeed do not substantially underestimate uncertainties due to missing higher-order corrections. At the same  Figure 7: Same as Figure 5, but for the difference in η and φ for each photon pair.
time, it is clear that the LO results including their uncertainties are completely insufficient to yield any meaningful prediction, and even NLO predictions can hardly be trusted when showing large discrepancies with data of several standard deviations.
We now move to studying differential distributions and the impact of NNLO corrections for higher centre-of-mass energies. Unfortunately, there has not yet been any 13 TeV measurement at the LHC, which is why we focus on theoretical predictions only. We obtain differential results for all collision energies considered in Table 2. In comparison to the results at √ s = 8 TeV, we do not find any substantial differences in the shapes of the K-factors. The corrections generally increase with energy, which we have already inferred from the results in Table 2, and, as expected, the jet activity increases. Since our results are not NNLO accurate any longer when requiring a jet, we focus on distributions that involve only the kinematics of the colour singlet final state. We note, however, that some observables intrinsically require jet activity in certain phase space regions, which is reflected by vanishing LO predictions. All relevant features will be discussed using the 13 TeV results as a reference.
In Figure 8 we present various differential distributions at 13 TeV, and we use a much finer binning to better resolve certain features. The upper plots of Figure 8 show the azimuthal difference between the hardest and the second-hardest photon (∆φ γ 1 ,γ 2 ) as well as between the second-hardest and the third-hardest photon (∆φ γ 2 ,γ 3 ). The hierarchy of the p T -ordered photons induces significant differences between the two cases. For LO kinematics, γ 1 and γ 2 need to recoil against each other since γ 3 does not carry sufficient energy to provide the recoil when the momenta of those two harder photons align. Correspondingly, γ 2 and γ 3 cannot be produced in back-to-back configurations at LO since these two photons need to recoil against the hardest photon γ 1 . As a consequence, the LO cross section vanishes for ∆φ γ 1 ,γ 2 < 2π/3 and ∆φ γ 2 ,γ 3 > 2π/3, respectively. Those phase space regions are filled only upon inclusion of real QCD radiation through higher-order corrections, which is required to overcome the kinematic constraints at LO. Accordingly, the NLO (NNLO) predictions in these regimes are effectively only LO (NLO) accurate, which is reflected by the increased size of both corrections and uncertainty bands. We find that back-to-back configurations of γ 1 and γ 2 are still preferred at higher orders, whereas the distribution of the azimuthal separation between γ 2 and γ 3 becomes much more uniform when adding higher-order corrections.
In the central plots of Figure 8 we show the invariant-mass and transverse-momentum distributions of the three-photon system. The invariant-mass distribution peaks around 100 GeV. Below the peak the distribution falls off steeply with a lower bound imposed by the phase space selection cut m γγγ ≥ 50 GeV. In that low m γγγ region radiative corrections increase quite strongly. By contrast, higher-order corrections become successively smaller in the tail of the m γγγ distribution, which is an important region for new-physics searches through small deviations from SM predictions. Around m γγγ = 1 TeV NLO and NNLO predictions become almost compatible within uncertainties. Also in the tail of the triphoton transverse-momentum spectrum NNLO corrections become successively smaller. We note that this observable vanishes for p T,γγγ > 0 at LO, and it diverges for p T,γγγ → 0 at any order in QCD perturbation theory. Only a suitable resummation of logarithmic contributions (cf. Ref. [71]) warrants a physical prediction at small p T,γγγ . The increase of NLO (NNLO) uncertainties in the p T,γγγ spectrum is again caused by the fact that the effective accuracy is decreased to LO (NLO) for this observable.
We conclude our analysis by studying the transverse-momentum spectrum of the hardest photon (p T,γ 1 ) and of the sum of the second-and third-hardest photons (p T,γ 2 γ 3 ) in the lower plots of Figure 8. The distribution in p T,γ 1 shows a rather similar pattern as the m γγγ distribution in terms  of the NNLO corrections. The distribution is cut at p T,γ 1 = 27 GeV in our fiducial setup, which induces an interesting behaviour in the p T,γ 2 γ 3 spectrum around that value at NLO and NNLO. The reason is that p T,γ 2 γ 3 = p T,γ 1 for LO kinematics, so that this distribution is not filled below 27 GeV at LO. At higher orders the spectrum develops a perturbative instability at this threshold caused by an incomplete cancellation of virtual and real contributions from soft gluons, which is logarithmically divergent, but integrable [90]. This instability slightly decreases when going from NLO to NNLO, but only a resummation of the relevant logarithmic contributions would yield a stable result. Practically, choosing a wide bin around the instability would alleviate this unphysical behaviour significantly. Note that such perturbative instabilities are present also for the p T,γ 1 γ 2 and p T,γ 1 γ 3 spectra, at the fiducial cut values imposed on p T,γ 3 and p T,γ 2 , respectively.
We stress again that LO results clearly fail to provide any reasonable prediction, and even NLO results strongly underestimate the cross section. Thus, one should bear in mind that care must be taken when determining new-physics contributions through higher-dimensional operators, see Refs. [38,40,41] for instance. Relying on LO cross sections for that matter might be insufficient because of the vastly inappropriate modelling of the cross section and distributions at this order. Even when relative BSM effects are computed by taking the ratio to the SM prediction at the same order, those effects are hardly trustworthy when based on a calculation at LO in QCD.
To summarize, we have presented a new calculation of the NNLO QCD corrections to the hadronic production of three isolated photons. This is the very first 2 → 3 process known at this accuracy, and we found our results to be fully compatible with the earlier calculation of Ref. [1]. Our calculation shows that even for highly non-trivial colour singlet processes q T subtraction is suitable to obtain NNLO predictions with systematic uncertainties at the few permille level, and that even complex five-point two-loop amplitudes can be implemented so efficiently that they do not pose a numerical bottleneck and can be evaluated directly during phase space integration. We have presented numerical results at various pp collision energies. Radiative corrections to this process are enormous: The NLO/LO K-factor is a factor of three (at 7 TeV) to five (at 100 TeV), while NNLO corrections turn out to be still as large as +50% (at 7 TeV) to +100% (at 100 TeV). NNLO is the first order to yield reliable predictions for triphoton production. Indeed, the comparison to 8 TeV data indicates substantial discrepancies with NLO results, while being in excellent agreement with the NNLO predictions. This is true not only for the fiducial rate, but also for differential distributions in the fiducial phase space. Our study of distributions at 13 TeV substantiates that NNLO accuracy is mandatory for precision phenomenology of this process in the future. We reckon that our results and our calculation 13 will be very useful for both measurements of the production of three isolated photons at the LHC and related BSM searches.