Three photon production to NLO+PS accuracy at the LHC

In this paper, we present the next-to-leading order predictions for three photon production in the standard model, matched to the parton shower using the MC@NLO formalism. We have studied the role of the parton shower on various observables and we show a selection of results for the 14\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$14$$\end{document} TeV Large Hadron Collider.


Introduction
A wealth of data from the Large Hadron Collider (LHC) and the Tevatron involving large number of leptons, gauge bosons and hadrons in the final state not only provides an ample opportunity to test the predictions of the standard model (SM), but it also constrains various physics scenarios in the beyond standard model (BSM). Signatures of BSM are often plagued by the large SM background and hence careful study of wide variety of SM processes has been under way [1][2][3]. Precise predictions for such SM processes are important as the quantum corrections are often comparable to the BSM effects. In addition, they are essential to reduce the theoretical uncertainties of the leading order (LO) predictions, which arise from the missing higher order quantum corrections through the renormalisation and factorisation scales. This necessitates the calculation of the next-to-leading order (NLO) quantum effects through Quantum Chromodynamics (QCD) radiative corrections to these SM observables at the hadron colliders. Presently, the phenomenological results for almost all physical processes of interest at the LHC are available at this accuracy due to the tremendous advancement in automating the calculation of the virtual and real emission contributions. However, the situation becomes more involved a e-mail: mandal@hri.res.in b e-mail: prakash.mathews@saha.ac.in c e-mail: ravindra@imsc.res.in d e-mail: satyajit.seth@saha.ac.in in the case of next-to-next-to leading order (NNLO) calculation due to a lot of technical difficulties. In order to improve the theoretical predictions in a consistent way, it is customary to take into account all the higher order contributions, which are important in the complementary kinematic regions of phase space corresponding to the phase space relevant for the fixed order evaluation. However, to correctly treat these effects one meets with a lot of technical problems. In practice, it can be approximated via the parton shower (PS) algorithm, which not only gives a reasonable estimate of these effects in the collinear kinematic regions of the phase space, but also provides a very realistic final state configuration. In other words, parton level predictions have to be gone through as regards such showering of multi partons and recombination of these partons into hadrons through a hadronisation mechanism in order to compare them against the experimental data. Such predictions require careful matching of results at various orders to avoid double counting. Thus, NLO SM results supplemented with parton showering can provide a more reliable as well as realistic predictions that can serve in testing various BSM scenarios. Till now, there are mainly two different algorithms in existence incorporating the matching of a NLO calculation to parton showers, namely MC@NLO [4] and POWHEG [5,6]. We shall adopt the MC@NLO algorithm here, which has already been implemented and completely automated in aMC@NLO [7,8].
In this article, we revisit the three photon production process at the LHC at NLO in QCD and present study the consequences of matching it with the parton shower. Triple-photon production provides a background to techni-pion production in association with a photon, where the techni-pion decays into a photon pair [9]. This process has already been studied at LO [10,11], as well as at NLO level [12] in QCD. We extend the analysis including the effect of parton shower to get a realistic estimate of various kinematical distributions. We quantify the improvement in the predictions at small trans-verse momentum regions of the final state particles and the stabilisation of the cross section against the variation of the factorisation and renormalisation scales. This paper is organized as follows: in Sect. 2, we have described the details of the calculation, mainly the virtual as well as the real emission contribution. The numerical results of the fixed order calculation together with the NLO+PS accurate results have been discussed in Sect. 3 and finally, we conclude in Sect. 4.

Calculational details
LO (O(α 3 )) contributions to the production of three photons at the LHC come from quark anti-quark annihilation processes. At NLO O(α 3 α s ) in QCD, we encounter virtual as well as real emission contributions resulting from an additional parton, namely quark or anti-quark or gluon. Virtual amplitudes are already at O(α 3/2 α s ), hence only the interference of them with the LO Born amplitudes will contribute to the NLO level. The real emission processes at NLO level come from two types of processes, namely gluon emissions from the LO processes and scattering of a quark (anti-quark) and a gluon producing three photons along with a quark (antiquark). The ultra-violet (UV) divergences, coming from the virtual contributions, and the infra-red (IR) divergences, originating from the virtual as well as real emission contributions, need to be removed through the addition of proper counter terms. The resulting IR-safe parton level cross section up to NLO can be written as The first term is the Born contribution; d P S 3γ is the phase space measure of the three photon final states and S({ p} 1,m ) is the observable function which depends on the kinematic variables through the momenta of the external particles i.e., p 1 , p 2 , . . . , p m . The second term corresponds to virtual corrections to the Born process. They are often divergent when the loop momentum becomes very large and these UV divergences are first regularised and then renormalised using the counter terms given in the third term. The fourth term represents the real emission contributions at the NLO level coming from parton emissions by the initial and/or final state partons. Due to the massless quarks, anti-quarks and gluons participating in the hard processes, both virtual and real emission contributions encounter soft and collinear divergences. The divergences coming from soft gluons and from collinear partons in the final state of the real emission processes get cancelled with those coming from the virtual processes. The remaining collinear divergences from the initial states are removed by adding mass counter terms given in the last term of Eq. (1). The details of obtaining UV renormalised virtual contributions are discussed in the next section. The real emission contributions and the corresponding mass counter terms are obtained with the help of MadFKS [13], a set of routines available in the aMC@NLO [7,8], which, along with our in-house FORTRAN routines for calculating virtual contributions, can provide results on an event-by-event basis in terms of four momenta of all the particles involved in the scattering process and we use them to obtain the observables that we require to study. In the following sub-sections, we sketch a systematic outline of the complete computational procedure.

Virtual contribution
The virtual contribution comes from the interference between the Born diagrams and the one loop corrected virtual diagrams. The number of virtual diagrams to order α 3/2 α s for the three photon production is 48. Up to permutations of the final state photons, we find 1 pentagon diagram, 2 box diagrams, 3 triangle diagrams and 2 bubble diagrams. We have used QGRAF [14] to generate both LO and NLO amplitudes. It generates the symbolic description of the Feynman diagrams in terms of propagators and vertices. We have written a FORM [15] code, which translates the output of QGRAF into a suitable format, that can be used for further symbolic manipulations. We have supplied Feynman rules, identities for Dirac gamma matrices and equations of motion through this code, and we have performed various simplifications at the amplitude level. The loop integrals are regulated using dimensional regularisation. Both Lorentz contractions and Dirac gamma matrix simplifications are done in n = 4 + ε space-time dimensions. Both UV and IR divergences appear as poles in ε and they have been calculated using the MS scheme. Writing the virtual contribution in the following way: col spin where M (0) is the Born amplitude and M V,( ) s' are the distinct topologies of virtual diagrams. We compute only one particular topology. Then the permutations of photon momenta and their polarisations gave us the remaining contributions.
The reduction of tensor integrals to scalar ones in n dimensions is done using the standard procedureà la Passarino-Veltman [16][17][18]. The tensor integrals that appear at one loop level are of the form where One can decompose the above tensor integral in terms of the scalar coefficients as follows: where the square bracket implies the non-equivalent symmetrisation by giving the full set of non-equivalent permutations. We have written a FORM code for the purpose of doing this tensor reduction and expressed the virtual contributions in terms of these scalar coefficients. As described in [19], these coefficients are related to the scalar integrals in different space-time dimensions in the following way: where is a generalized scalar integral in shifted space-time dimension. These integrals in the shifted dimensions can be expressed in terms of integrals in n dimensions using the dimensional recurrence relations discussed in [20,21]. In this approach, inverse Gram determinants that result from the recurrence relations, often spoil the numerical stability of the integral. There exist a handful of solutions to this problem in the literature [22][23][24][25][26][27][28][29][30][31][32]. Recently, an elegant approach has been put forward in [33], where the authors have found signed minor algebraic relation, which avoids the appearance of inverse Gram determinants and thereby introducing a set of higher dimensional scalar integrals to cope with the small Gram determinants. These higher dimensional scalar integrals have been evaluated numerically after employing a series expansion in the small Gram region. This whole algorithm has been implemented in the numerical package, named PJFry [34,35], which we use to evalu-ate numerically the scalar coefficients of the tensor integral for every phase space point in n dimensions. PJFry reduction library uses QCDLoop [36] and OneLOop [37] to evaluate the scalar integrals in 4 dimensions. In order to validate our FORM codes, namely those ones that perform conversion of output of QGRAF to FORM readable symbolic expressions, reduction of tensor integrals to scalar coefficients and also to validate FORTRAN routines, which evaluate the virtual contributions numerically using PJFry, we recalculated the virtual corrections of the diphoton production process in both SM and BSM to order α s . We compared our results thoroughly against the results presented in [38][39][40] and found an excellent agreement between the two. Using our FORM codes and FORTRAN routines along with the publicly available packages, viz. QGRAF, PJFry, QCDLoop and OneLOop, we have evaluated the virtual contributions to the three photon production process at O(α s ) level. We find that after UV renormalisation, the IR poles, namely double and single poles in ε, are in accordance with the expectation. We express the virtual contribution of the three photon production in a form suitable for further analysis as follows: where α s is the strong coupling evaluated at the renormalisation scale μ R , s is the partonic center-of-mass energy and the colour factor is C F = 4/3 for SU (3). dσ denotes the finite virtual contribution that has been computed numerically. Note that the IR poles in ε are in agreement with the universal behaviour of soft and collinear partons.

Real emission contribution
Real emission contributions come from gluon emission from the Born processes as well as from the scattering of a quark/anti-quark and a gluon producing a quark/anti-quark and three photons. We use aMC@NLO [7,8] framework not only to compute these contributions along with the mass factorisation terms required to remove the initial state collinear singularities, but also to obtain the NLO results matched with PS. Within aMC@NLO, the stand-alone package MadGraph [41,42] generates all the required matrix elements both at LO as well as at NLO level. As already discussed in Sect. 2.1, we have prepared a set of external codes to deal with the virtual correction part and made an interface to implement it within MadFKS [13], which separates out the soft and collinear configurations in the real emission processes using the FKS subtraction scheme [43,44] and provides IR-divergent and IR-safe contributions sepa-rately along with the mass factorisation terms, which take part in removing the initial state collinear singularities coming from the virtual and the real emission processes. In the FKS subtraction scheme, the phase space is partitioned in such a way that each partition contains at most one soft and one collinear divergences. This is done by introducing a set of positive-definite S i j functions, where the S i j s' are chosen in such a way that they vanish in all singular limits not related to: (i) a particle i becoming soft, ( ii) particles i and j becoming collinear, obeying the restriction that the sum over all such pairs must be equal to identity. This ensures that each term of the sum is finite throughout the phase space, except when the energy of particle i goes to zero or particles i and j become collinear. Now, after finding out the exact position of the divergences for a given partition, the generalized plus distribution is used to regulate them. All these steps are systematically automated in MadFKS within the MadGraph5 environment. We have explicitly checked the cancellation of the soft and collinear divergences among the virtual, real and mass factorisation terms at different regions of the phase space thereby confirming the perfect implementation of all the above mentioned external inputs within the aMC@NLO framework. The events that are generated using aMC@NLO also include the Monte Carlo counter terms to take care of the MC@NLO matching and thereby preventing the occurrence of any double counting at the time of matching to PS. These events are then showered by HERWIG [45][46][47], PYTHIA [48] parton shower to get the realistic events.
Photons are produced not only at the partonic level, but also through the fragmentation of partons into photons and a jet of hadrons can often be collinear to them. This necessitates the inclusion of non-perturbative fragmentation functions. At NLO level, the QED collinear divergence can arise when one of the final state parton becomes collinear to a photon. This can be factorised in a universal manner and then removed by adding counter terms, which renormalise the fragmentation functions, thereby bringing in a scale dependence at the partonic cross sections through the fragmentation functions, which is known as fragmentation scale. An alternate isolation criterion has been proposed in [49], using which one can obtain an observable in which fragmentation contribution is minimised and at the same time, the IR-safety of that observable is guaranteed. We call it Frixione isolation here after and use this isolation for our analysis. It works in the following way: define a cone centered around each photon with a radius R in the rapidity-azimuthal angle (η-φ) plane, 2 . Now, it is demanded that the sum of hadronic transverse energy H (R) inside any concentric circle of radius R < R γ would be less than an amount given by the function H (R) max . This function can be chosen in such a way that less and less hadronic energy is allowed as we move closer to a given photon. Because of the fact that H (R) goes to zero as R → 0, the partons that are collinear to photon are removed while the soft partons are kept intact thereby guaranteeing the QCD IR-safety. For our analysis, we have taken the following canonical choice for H (R) max , i.e.: where E γ T is the transverse energy of the photon and R γ , γ , n are three parameters that are to be set while applying this isolation criterion. 1

Numerical results
In this section, we present the results for various kinematic distributions relevant to the production of three photon in SM at the LHC with the center-of-mass energy √ S = 14 TeV. Here we list the input parameters used for the whole computation: For the fixed order (N)LO calculation, we have taken the following choices of cuts: rapidity of each photon |η γ | < 2.5, separation between any two photons in the (η-φ) plane In addition, we have studied a variety of differential distributions applying two types of cuts on the transverse momentum of each photon i.e., P γ T > 20 GeV and P γ T > 30 GeV in the fixed order analysis. Unless stated otherwise, we consider P γ T > 30 GeV as the generic choice of the cut on the photons' transverse momenta. Parameters involved in the Frixione isolation are set as R γ = 0.7, γ = 1 and n = 2.

Fixed order analysis
In Table 1, we have shown the results of total cross sections for fixed order LO and NLO using the central choice of μ F and μ R for two different P γ T cuts. To begin with, we present some distributions of few selective kinematical variables at fixed order LO and NLO. Photons are ordered according to their transverse momentum. The hardest photon with maximum transverse momentum is denoted by γ 1 . Like wise, γ 2 represents the second hardest photon and the softest photon is labelled as γ 3 . In Fig. 1, we have shown transverse momentum distribution of γ 1 at LO and NLO in the left panel and in the right panel, distribution of invariant mass of the three photon system has been plotted. The lower insets show the bin-by-bin distribution of the K-factor for the corresponding observable. We find that, for low transverse momentum, the K-factor is large, as it is due to the fact that the recoil against the extra parton helps to fulfill the transverse momentum cut, which was not possible at LO. The left panel of Fig. 2 shows the distribution of the rapidity of the hardest photon, whereas the rapidity of the three photon system is shown in the right panel. The distribution of the K-factor is shown for the corresponding variables in the lower insets. Unlike Fig. 1, the K-factors in Fig. 2 appear to be mostly steady, indicating the affinity of these observables towards the photons having fairly high transverse momenta. All the above distributions show a substantial effect of radiative corrections on this process. This is mainly because of the inclusion of new subprocesses at the NLO, as quarkgluon subprocesses begin to contribute at this order and due to the enhancement in the phase space. In Fig. 3, we have plotted the separation between the ordered photons in the (η-φ) plane obeying the selection cut: R γ i γ j > 0.4, where i, j = 1, 2, 3. We have checked the rapidity differences between these photons are quite small. Therefore, the peaks arising in these distributions near the angle π (180 • ), suggest that the emitted photons are mostly back-to-back. The hardest photon γ 1 is separated form the softest one i.e., γ 3 , by at least R γ 1 γ 3 = 1.6 at LO, whereas at NLO they can be very close as permitted by the selection cut due to the emission of an extra radiation at this level.
Besides, we have checked the effect of variation of Frixione isolation parameters i.e., R γ , γ and n. Though Frixione isolation has no effect on the LO cross section, the dependency of the NLO cross section on these isolation parameters is shown in Table 2. From Eq. (8), it is evident that the NLO cross section increases when R γ decreases and it also increases with increasing γ . In Fig. 4, we have shown the transverse momentum distribution of the hardest photon by  Table 2, as well as from Fig. 4, that, for a fixed choice of n, the NLO cross section is large for R γ = 0.4 and γ = 1, whereas it Table 2 Total cross sections for the three photon production at the LHC for various Frixione isolation parameters. We have taken p becomes much smaller for R γ = 1 and γ = 0.1, indicating the fact that a smaller R γ increases the cross section when γ is larger [51]. It is also clear from Table 2 that the effect of varying γ , keeping n and R γ fixed, is quite minimal. Similar studies with changing the value n = 2, provide same kind of distributions analogous to Fig. 4.

Discussion on NLO+PS
In this section, we compare the fixed order NLO result with the NLO results matched with PS (NLO+PS) with two different showering algorithm, namely HW6 and PY6. For the showering purpose, parton level events are generated using very loose cuts: P γ T > 15 GeV, |η γ | < 2.7, R γ γ > 0.3 with the following Frixione isolation parameters: R γ = 0.4, γ = 1 and n = 2. We have explicitly checked that the events thus produced, remain unbiased in total rates and differential distributions after showering and hadronisation for this The scale dependencies of the results are calculated by varying μ F and μ R independently around the central value μ F = μ R = M γ γ γ via the following assignment: μ F = ξ F M γ γ γ and μ R = ξ R M γ γ γ , where ξ F and ξ R are varied in the range [1/2,2] independently. Various ratios of μ F , μ R and M γ γ γ that appear as arguments of logarithms in the perturbative expansion to NLO are within the range [1/2,2]. The scale uncertainty band is the envelope of the results obtained by varying this ξ F and ξ R within this range [40]. The PDF uncertainties are estimated with the Hessian method, as given by the MSTW [52] collaboration. We have plotted fractional uncertainty, which is defined as the ratio of the variation about the central value divided by the central value, being a good indicator of the uncertainties. These uncertainty bands can be generated automatically at the time of parton level event generation by storing additional information, sufficient to determine via a reweighting technique, at no extra CPU cost within the aMC@NLO framework as described in [53].
We have shown log 10 P γ γ γ T distribution for HW6 and PY6 together the fixed order NLO result, in Fig. 5. It is clear that at low P γ γ γ T values, NLO+PS (for both HW6 and PY6) result shows the effect of all order resummation of the large logarithms, hereby suppressing the cross section leading to a meaningful value, while the fixed order NLO result diverges for P γ γ γ T → 0. At low P γ γ γ T , PY6 result is different from the HW6 result as the soft and collinear emissions constituting the parton shower are treated differently. PYTHIA generates softer spectra than HERWIG in this region and as a result of this, these two showers show different behaviour as expected [54]. At high P γ γ γ T , the NLO fixed order and NLO+PS (for both HW6 and PY6) results are in agreement as in this region, the hard emissions are dominant and they are correctly described by the NLO hard cross section. In the middle and lower insets of Fig. 5, we have presented the fractional scale and PDF uncertainties of the NLO+PS result for HW6 and PY6, respectively, which increase with increasing P γ γ γ T [54]. We do not find any significant differences in the case of studying fractional uncertainties using these two different showers. Therefore, in the rest of the figures, we present the fractional uncertainty plots only for HW6.
We now present the results for various kinematical distributions to NLO accuracy, matched with PS (labelled as NLO+PS), for both HW6 and PY6 with the specified analysis cuts. We have adopted a consistent pattern for all the rest of the distributions. In each case, within the main frame, three curves corresponding to the distributions in fixed order NLO (solid red) and NLO+PS using HW6 (dashed blue) and NLO+PS using PY6 (dotted black) are shown. The middle inset shows fractional scale uncertainty (dashed cyan) and fractional PDF uncertainty (solid violet), while the lower inset shows the ratio of NLO+PS and NLO for HW6 (dashed blue) and for PY6 (solid black). In the left panel of Fig. 6, we have shown the plots for transverse momentum distribution of the hardest photon and the right panel shows the distribution of the invariant mass distribution of the three photons. We do not find much difference in the results of two showers HW6 and PY6. In both distributions, NLO results are very little larger than the NLO+PS results and this is due to the fact that the QCD radiation becomes softer when we demand all the photons to satisfy a high P γ T cut (i.e. P γ T > 30 GeV) and damping of PDFs at large Bjorken x values further subdue its effect at the parton level, whereas aMC@NLO produces more events with hard and central jets resulting in the suppression of these distributions after showering. In Fig. 7, we have depicted the plot showing rapidity distribution of the three photon system and as expected, we observe that the NLO result is slightly harder than the NLO+PS result. However, the ratio in the lower inset shows that PY6 generated events give a larger contribution than HW6 in the large rapidity region indicating that PY6 produces a significantly larger number of radiations than HW6 in the full kinematically available phase space.

Conclusion
Precise and realistic predictions of both signal and background processes at hadron colliders are now possible due to tremendous developments in the computational methods and the availability of the state of the art computational tools. We have used packages, namely QGRAF, PJFry and aMC@NLO, to study the three photon production process at the NLO level in QCD for the LHC taking into account the parton shower effects and realistic experimental cuts. In addition, we have developed some codes that build the interfaces among these different analytical and numerical tools. We have plotted different kinematic observables and discussed the consequences of showering the fixed order NLO results with the two different showering algorithms HERWIG and PYTHIA. We have also discussed the effects of scanning over the Frixione isolation parameters on the NLO cross section. We find that our predictions are less sensitive to scale uncertainties and choice of PDFs and hence more suited for direct comparison with the data from the experiments.