Exploring QED Effects to Diphoton Production at Hadron Colliders

In this article, we report phenomenological studies about the impact of O(α) corrections to diphoton production at hadron colliders. We explore the application of the Abelianized version of the qT-subtraction method to efficiently compute NLO QED contributions, taking advantage of the symmetries relating QCD and QED corrections. We analyze the experimental consequences due to the selection criteria and we find percent-level deviations for Mγγ>1TeV. An accurate description of the tail of the invariant mass distribution is very important for new physics searches which have the diphoton process as one of their main backgrounds. Moreover, we emphasize the importance of properly dealing with the observable photons by reproducing the experimental conditions applied to the event reconstruction.


Introduction
Direct photon production in hadronic collisions provides an excellent opportunity to test Standard Model (SM) predictions and explore potential new physics effects, due to its very clean experimental signature. Photons are also very relevant to study the internal structure of colliding particles [1,2] because they can shed light on the parton level kinematics, leading to new interesting paths to unveil the hidden symmetries of nature.
In this phenomenological context, one process of special relevance is the production of a pair of prompt photons (i.e., diphoton production), which was one of the key channels for the discovery of a new particle compatible with the SM-Higgs boson [3,4]. This indicates that direct photon production might also play an important role in the precision physics program for future colliders [5][6][7][8][9][10][11][12][13].
Due to its relevance, the direct diphoton production has been extensively studied and its perturbative corrections were computed in several models. At the lowest order (LO), only the qq partonic channel is available. The first quantitative correction appears at the next-to-leading order (NLO) in the strong coupling constant α S [14][15][16][17], which includes the qg channel and provides a noticeable enhancement due to the gluonic luminosity in the low-x region. At the next-to-next-to-leading order (NNLO) in α S (i.e., O(α 2 S )), the gg channel becomes available. Its contribution is comparable to the Born cross-section, but the total NNLO correction is still dominated by the qg channel. The NNLO QCD corrections were first calculated in Ref. [18] and later in Ref. [19]. The QCD resummed corrections have been recently obtained up to NNLL accuracy [20], thus improving the phenomenological description by including the effects due to low-energy gluon radiation. Recently, articles on diphoton production at NNLO were published, studying in detail the isolation criteria [21][22][23][24][25][26][27][28], the weight of the fragmentation component [21,22,24,26], the reliability of the prediction under symmetric cuts on single final photons [21], the scale sensitivity of the process and the hybrid isolation criteria [29,30], among other features of this important process. In addition, the diphoton background could be potentially interesting as an indirect measurement of the top-quark mass [31,32].
Besides the direct photon production from the hard subprocess, photons can be generated through the fragmentation mechanism, which involves non-perturbative transitions from QCD partons. The complete NLO single-and double-fragmentation contributions are implemented in DIPHOX [14]. The measured photon cross-sections use isolation prescriptions in order to reduce the large reducible background (associated to photons that are faked by jets or produced by hadron decays). Two of these prescriptions are the standard cone isolation and the smooth cone isolation [33] prescriptions. While measured cross-sections rely on the standard cone criterion, theoretical calculations can be greatly simplified with the smooth cone prescription. It is worth noticing that both algorithms produce similar results [22,34] for those isolation parameters implemented by the experiment.
Other important effect to be considered in the context of photon production is the presence of electroweak (EW) corrections. Leaving aside the large gluon luminosity at the LHC, a simple power counting shows that the NLO QED corrections could compete with the NNLO QCD contributions because α 2 S ≈ α. For this reason, the higher-order QED and EW effects deserve a serious treatment. Several studies about the EW corrections to gauge boson production [35][36][37][38][39][40] and diphoton production [41][42][43][44] showed small but still non-negligible contributions, which might play a crucial role in the context of the precision physics program. Moreover, we have explored the impact of pure QED and mixed QCD-QED corrections in the DGLAP equations [45][46][47] and we found that they could provide non-negligible effects to the PDF evolution. In particular, the NLO QED corrections to the diphoton production process are sensitive to the photon PDF through the qγ channel, which might propagate the information related to this distribution to the final cross-section result.
From the point of view of new physics searches, the high-energy region of the diphoton spectrum could be used to impose constraints on many models. Since no new particles are predicted by the SM beyond the TeV scale, the inclusion of higher-order EW and QCD effects could only lead to an enhancement/decrease of the distributions without introducing any peak. Thus, a proper understanding of all the SM effects for this region might allow to impose tighter constraints to BSM models involving a continuum spectrum at high-energies (for instance, large extra-dimensions or Randall-Sundrum with composite fermions [48]).
The purpose of this article is to describe interesting phenomenological features of the NLO QED corrections to diphoton production, which might be useful for deeper studies of this process including higher-order mixed QCD-EW corrections. We start in Section 2 describing the theoretical framework used to perform the computation, centering the discussion in the different partonic channels contributing to the process. Special emphasis is put on the treatment of the electromagnetic coupling, because it leads to noticeable corrections to the theoretical predictions. Then, in Section 2.1, we explain the experimental cuts applied, focusing on the photon isolation and reconstruction algorithms. In Section 3, we present our predictions for LHC, highlighting the effects due to the ordering algorithms used to deal with the radiated photons. Besides that, we briefly comment about the implementation of a clustering (i.e., merging) algorithm for reconstructing the photons and the dependence of the results on the use of different PDF sets, and the photon PDF in particular. Finally, in Section 4 we present our conclusions and a brief discussion about the importance of higher-order QED corrections in future experiments.

Computational Details
In our calculations we focus on the differential cross-section for the production of a prompt-diphoton system in hadron-hadron collisions. By virtue of the factorization theorem, the differential NLO QED cross-section for the process pp → γγ + X is given by where f a/h (x i , µ 2 F ) denotes the PDF for partons of flavor a = {q, g, γ} to be found inside the hadron h, x i are the longitudinal momentum fractions and µ F is the factorization scale. Additionally, we introduced the supraindices (i, j) to denote the perturbative expansion in QCD-QED; explicitly, (i, j) corresponds to the O(α i S α j ) correction to the Born process. dσ ab→γγ+c are the real corrections to the Born subprocess due to QED extraradiation, where the symbol dσ denotes the regularized partonic cross-section, i.e., without the infrared (IR) and ultraviolet (UV) singularities. In Equation (1), S i is the measurement function that defines the IR-safe observable at parton level for the Born (i = 2) and realradiation (i = 3) kinematics. At Born-level, the only available subprocess is qq → γγ, whilst the real-radiation contributions at NLO are associated to the subprocesses therefore the sum performed over the c parameter in Equation (1) refers to the precedent two channels. Notice that, at O(α 3 ), the qγ channel starts to contribute and the cross-section becomes sensitive to the photon PDF. In order to explicitly implement this computation numerically, we must properly define the regularized cross-sections, which implies canceling the UV divergences through renormalization and the IR singularities with a proper real-virtual combination. To tackle the IR regularization, we will rely on the q T -subtraction formalism [49,50], accordingly modified to deal with QED corrections in a fully consistent way. Within this formalism, the cross-section associated to the production of any neutral final state F with transverse momentum q T is separated into a regular and a singular contribution in the limit q T → 0. Whilst the regular part is finite and process-dependent, the singular contribution possesses a universal structure, given by [50] dσ (sing.) F where S c is the Sudakov form factor, [H F C 1 C 2 ] cc,a 1 a 2 is the hard-collinear factor and [dσ cc→F ] is related to the Born level cross-section. The Sudakov factor embodies all the information related with the emission of low-energy photons from the parton c entering in the hard-process cc → F. It is obtained from the exponentiation of photon emissions, and is given by which depends on the flavor of the parent parton c and the resummation coefficients A c (α) and B c (α). In Equation (4), Q 2 is evaluated at scales typically present in the hard process, b is the impact parameter and b 0 = 2e γ E (where γ E = 0.5772 . . . is the Euler-Mascheroni constant). The resummation coefficients A c (α) and B c (α) can be computed in perturbation theory, It is worth noticing that these ideas can be applied for dealing with a mixed QCD-QED expansion, which justifies the inclusion of a double-superscript to keep track of the perturbative order in each theory [40]. In the case c = q, the first terms of the expansions are with e q the electromagnetic (EM) charge of the quark q; A For the purpose of computing NLO QED corrections, we only need A (0,1) q : the remaining higher-order coefficients are available in Refs. [40,51].
On the other hand, the symbolic factor H F C 1 C 2 in Equation (3) for the qq annihilation channel is given by where the functions C q a are universal and the renormalised hard function H F q contains process-dependent contributions related to the virtual amplitudes. In both cases, they can be expanded in perturbative series, i.e., where we are using the double-superscript notation mentioned before. The only nonvanishing C (0,1) coefficients are with N C = 3 the number of colors. It is worth emphasizing that the structure of Equation (7) is different when dealing with vector particles (i.e., c = g or c = γ) due to the presence of spin correlations. More details about the mixed QCD-QED q T -subtraction/resummation formalism can be found in Refs. [40,51]. For the particular case of the diphoton production, we can show that the differential cross-section can be written as where the finite process-dependent pieces coming from the virtual amplitudes (and the Born cross-section) are contained in the hard coefficient function H (0,1) γγ . In this equation, the symbol ⊗ is understood as a convolution of momentum fractions and sum over flavour indices of the partons [49][50][51]. dσ While in the QCD formalism the undetected final state X can contain any number of QCD partons (i.e., quarks and/or gluons), in the QED case X contains any number of photons. The two terms in the r.h.s. of Equation (11) are finite, which implies that they can be numerically computed. However, both dσ (0,0) γγ+X and dσ (0,1) CT diverge in the limit where the transverse-momentum of the system γγ vanishes (i.e., p γγ γγ is extracted in a process-independent way from the O(α) process-dependent corrections to the scattering amplitudes. In this case, the only contribution to this hard coefficient comes from the qq channel. Thus, the regularized hard coefficient is given by with v = −û/ŝ the ratio of the partonic Mandelstam variables. This expression is obtained from the hard factor provided in Ref. [50] after applying the Abelianization algorithm described in Refs. [45,46]. Finally, it is important to mention that we are using n F = 5 massless quarks as active flavors and we include all the SM charged-fermions inside any closed fermionic loops. To be more precise, the QED beta function reads where q (l) denotes all the possible flavors of quarks (leptons) in the SM with their corresponding EM charges, e q (e l ).

Event Selection, Isolation Prescription and Setup of the Calculation
The present NLO QED computation is implemented by modifying the numerical program 2γNNLO [18], which originally includes up to NNLO QCD corrections to promptdiphoton production. At this perturbative order in QED, the Born + virtual component only contains two photons, whilst the real-emission part introduces contributions with up to three observable photons in the final state (i.e., through the subprocess qq → γγγ). Since we are interested in computing pp → γγ + X, we order the transverse-momenta of the triphoton system, thus reproducing the selection criteria applied by the experiment in direct diphoton measurements. It is interesting to emphasize that the momentum ordering of the photons introduces a dynamical cut in the available real-emission phase space, as we will explain in Section 3.
As it was stated in the Introduction, we apply the smooth cone isolation prescription [33]. This criterion differs from the one applied by the experiment (standard cone) but, for commonly used isolation parameters by the LHC, both criteria lead to similar results [22,34]. Explicitly, given a final state photon, we build a cone around it, whose radius in the (η − φ) plane is where we are using the well-known notation η and φ for the pseudo-rapidity and azimuthal angle, respectively (see Ref. [21] for more details). Then, we require that the partonic energy deposited inside the cone fulfils [33] with ∆R ≥ r the radius of the fixed outer cone used to initialize the isolation algorithm, E T,max the maximum allowed deposited transverse energy and n an arbitrary isolation parameter. In this work, we choose n = 1. We show results for two different center-of-mass energies using the guidelines applied by the ATLAS and CMS experiments. Explicitly, for √ s = 7 TeV, we choose those events where p harder T ≥ p H = 25 GeV and p softer T ≥ p S = 22 GeV, restricting both photon rapidities to satisfy |y γ | ≤ 1.37 and 1.52 ≤ |y γ | ≤ 2.37 [52]. The isolation parameter was fixed to E T,max = 4 GeV. On the other hand, for √ s ≥ 13 TeV we required p harder T ≥ p H = 40 GeV and p softer T ≥ p S = 30 GeV whilst the rapidity range was slightly modified to |y γ | ≤ 1.37 and 1.56 ≤ |y γ | ≤ 2.37, and we imposed E T,max = 11 GeV in the isolation criterion [53]. For all these setups, the separation between the hardest photons in the (η − φ), ∆R γγ plane must to be greater than ∆R = 0.4.
In order to reproduce the experimental measurement procedure, we consider the implementation of a photon-clustering algorithm. From the experimental point of view, it is not possible to distinguish events with quasi-collinear photons due to the finite granularity of the detectors. Thus, if they are produced within a cone of radius r < ∆R min γγ , they are identified as an unique particle. For instance, for the ATLAS detector [54], this value ranges from 0.05 to 0.075, although a conservative estimate of the minimal detectable separation could be set to ∆R min γγ = 0.1. From the theoretical side, it is possible to implement such a procedure by working with the parton-level kinematics. In particular, notice that this algorithm will be activated only for the qq → γγγ channel, and that only one pair of quasi-collinear photons is allowed due to the other kinematical cuts (i.e., angular separation between resolved particles). Schematically, we consider two variants of the clustering/merging procedure, which are defined in the following way: 1.
Compute the distance r γ i γ j among the final-state photons, where i, j = 1, 2, 3 at NLO.

2.
If min r γ i γ j ≥ ∆R min γγ for all {i, j} pairs, all the photons can be resolved and the process is described with a full 2 → 3 kinematics as normal.

3.
If the photon k is isolated but r γ i γ j ≤ ∆R min γγ , then the photons i and j are detected inside the same bin. Hence, we define the merged momenta as which corresponds to applying the E 0 and E-schemes for QCD jets, respectively, [55,56].
In both cases, the definition of the original algorithm was slightly modified to deal with photons instead of jets.

4.
Order the final state particles according to their new transverse momentum, and impose the corresponding cuts to select the events.
Since the matrix element is finite when two photons become collinear, there is some freedom in the implementation of the algorithm. In the E 0 -scheme, we have p 2 ij = 0 and the energy measured in the detector agrees with the total sum of the diphoton system. Moreover, the direction of motion of the merged particle corresponds to the sum of momenta, with a proper re-scaling. In this way, the reconstructed photon is identified as an on-shell particle with a well defined three-momentum. However, this approach does neither fulfill momentum conservation nor Lorentz invariance. On the other hand, the Escheme is compatible with the last two properties, but the merged momenta is not massless (i.e., the on-shell condition is not fulfilled).
Finally, we comment on the treatment of the electromagnetic coupling α and its running. In the first place, it is a well-known fact that the running effects arise from the renormalization of the interaction vertex. If the Born contribution only contains on-shell final-state photons, when collecting all the possible unrenormalized one-loop contributions we realize that there is a partial cancellation of UV singularities between self-energies and vertex corrections. This cancellation avoids the presence of UV-logarithmic terms within the renormalized one-loop amplitudes; including a running coupling would add O(α 2 ) corrections proportional to log(µ R ) (with µ R the renormalization scale) that might alter the high-energy behavior of the physical predictions.
On the other hand, when the LO process involves initial-state photons, the NLO QED corrections should include the one-loop running for α evaluated at the factorization scale µ F [57]. Otherwise, some terms proportional to log(µ F ) would remain un-canceled and enhance artificially the scale-dependence of the hadronic cross-section. As pointed out in Ref. [57], this behaviour is due the presence of the photon PDF, which introduces a dependence on µ F , even if photons entering the partonic processes are kinematically on-shell. This dependence is controlled by the splitting kernel P γγ , that at the lowest order is given by [46] Notice that this expression involves the QED beta function at one-loop, since this splitting function is directly extracted from the one-loop photon self-energy. In consequence, there is an interplay between the renormalization and factorization corrections that require the proper inclusion of the running effects to achieve a consistent result.
As a consequence of the precedent observations applied to the particular case of diphoton production, we claim that we should fix α = 1/137 in order to provide a conservative estimate of the NLO QED corrections. However, it is worth emphasizing that we can introduce a partial running for α to estimate effects due to missing higher-order terms. Explicitly, since the LO is O(α 2 ), we propose to modify the NLO behavior of the coupling according to thus including the running effects without overestimating the correction by introducing additional logarithms. So, we consider the frozen coupling scenario as the default one, and we also take into account the full hadronic evolution given by Refs. [58][59][60] with the purpose of studying the propagation of uncertainties due to the choice of α.

Results and Discussion
In this section, we present the numerical results for diphoton production comparing the strength of the NLO QED corrections with the well-known NLO QCD ones. In the following, we denote diphoton system as the system composed of the two hardest photons. We show results for the LHC at two different energies, √ s = 7 TeV and √ s = 13 TeV, making use of two different sets of cuts for the transverse momentum of the diphoton pair. For √ s = 7 TeV we use nearly symmetrical cuts, whilst for √ s = 13 TeV we apply asymmetrical ones as detailed in the previous section.
Concerning the PDF sets, we perform the QED computations with LUXqed [61,62] unless otherwise specified. Since LUXqed uses PDF4LHC [63] to describe the QCD parton distributions, we consistently use it when computing the pure QCD contributions. Another set containing the photon PDF is NNPDF3.0QED [64,65]. We provide also results computed with the later set in order to compare the effects of different implementation of photon PDFs.
The default renormalization (µ R ) and factorization (µ F ) scales are set to the value of the transverse invariant mass of the diphoton system, i.e., As explained in Section 2, the default value for the QED coupling is fixed to α = 1/137 and we consider a partial running to estimate the effects due to missing O(α 4 ) corrections. On the other hand, we also study the uncertainties originated by missing QCD-QED higherorder contributions. Starting from the formal additive result for the NLO QED corrections, we define an estimate of the mixed QCD-QED corrections with the multiplicative approach, i.e., dσ QCD×QED The multiplicative ansatz is described in Reference [66] and is motivated by the fact that the IR-divergent structure of the mixed QCD-QED terms factorizes as in Equation (23). However, this approach fails to describe the finite pieces of the total contribution because it does not take into account non-factorizable QCD-QED terms present in the full computation.
After describing the details of the implementation, we start the analysis of the results. In first place, we consider the invariant mass distribution of the two hardest photons for √ s = 7 TeV and √ s = 13 TeV in Figures 1 and 2, respectively. We compare the impact of the NLO QED corrections (δσ QED NLO = dσ LO δ QED NLO ) with the NLO QCD (δσ QCD NLO = dσ LO δ QCD NLO ) and the corresponding Born-level contribution. The two different channels that are part of the total NLO QED correction, i.e., qq and qγ, are separately shown in the plots. The NLO QED corrections are moderate to tiny in the whole invariant-mass range. To quantitatively describe them, we define the ratio that reaches up to 4 % for M γγ 2 TeV. σ QED NLO is dominated by the qq channel in almost the entire M γγ range. The only exception is the low-mass region, which is forbidden at the LO due to kinematical constraints. The kinematics of the LO subprocess is determined by the two final-state photons. Since the only kinematical configuration allowed by the Born contribution is such that the two photons are back-to-back, p γγ T = 0 (i.e., p harder T = p softer T ) and therefore, only the value of p harder T is effective as cut at LO, Specifically, M γγ ≥ 50 GeV for photon-pair production at √ s = 7 TeV (and its associated LHC cuts) and M γγ ≥ 80 GeV for diphoton production at √ s = 13 TeV. Perturbative calculations in regions around unphysical fixed-order thresholds are known [67] to be generally affected by perturbative instabilities at higher orders. The invariant mass distribution is very steep in the kinematical region around the LO threshold and even the effect of little instabilities is amplified by the large slope of dσ/dM γγ . Therefore, in order to obtain reliable predictions around the LO threshold in dσ/dM γγ it is enough to consider a bin size of roughly 2 GeV since we are dealing with integrable singularities [21,67].
The suppression of the qq channel in the low-mass region is an artifact of the ordering in the transverse momentum of the final state photons. This ordering procedure reproduces the selection algorithm implemented by the LHC. It is worth noticing that we can describe quantitatively the interplay among the kinematical cuts, the ordering and the position of the threshold in the qq channel. Let us consider a system with three massless on-shell particles in the center-of-mass frame. After ordering, we have with p i = p T,i + p L,i and p T,i · p L,i = 0. In fact, the transverse momentum conservation implies p 2 T,3 = | p T,1 + p T,2 | 2 , and the ordering condition leads to where θ 12 is the angular separation in the transverse plane. Since p T,2 ≤ p T,1 by definition, we conclude that 2π/3 ≤ θ 12 ≤ π. So, we explicitly find a constraint in the angular separation of the transverse momenta of the hardest particles by imposing an ordering condition. Of course, such a restriction is not present when selecting randomly two photons and using their momenta to classify the event, as we did for the QCD corrections to pp → γγ + X. Moreover, from Equation (27) we can appraise that p T,2 ≤ p T,1 ≤ 2p T,2 , which constitutes another constraint. On the other hand, it is crucial to notice that imposing both the transverse momenta ordering and the p T cuts detailed in Section 2, we end up with a cut in M γγ . To find the minimum allowed value for M γγ distributions measured at ATLAS with √ s = 7 TeV, we choose p T,1 = 25 GeV and p T,2 = 22 GeV. Then, inserting these values in Equation (27), we find θ min 12 ≈ 0.6924 π , that increases the minimal angular separation between the transverse momenta of the hardest photons. Besides this, the invariant mass of the hardest subsystem is given by where φ 12 is the angle between p 1 and p 2 . This angle is measured in the plane containing both three-vectors, which is not equivalent to the transverse plane to the colliding axis: thus, in general, φ 12 = θ 12 . However, we impose φ 12 = θ 12 in order to find the minimal invariant mass configuration. This choice is equivalent to settle the production of the triphoton system in the transverse plane, which leads to E 1 = p T,1 , E 2 = p T,2 and M min γγ = 2 p T,1 p T,2 (1 − cos θ min 12 ) ≈ 41.53 GeV .
Below these limit, the qq → γγγ channel is not compatible with the experimental cuts. Analogously, when dealing with ATLAS measurements at √ s ≥ 13 TeV, we find These results explain the behavior of the distributions shown in Figures 1 and 2, in particular, the steeply suppression of the cross-section in the region M γγ ≈ 41 GeV for √ s ≥ 7 TeV and M γγ ≈ 63 GeV for √ s ≥ 13 TeV. In Figure 3, we present the results for the transverse momentum distribution of the two hardest photons for √ s = 13 TeV. We compare the NLO QCD and NLO QED corrections, as well as the individual channels contributing to the last one. Notice that in the whole p γγ T range (with the exception of p γγ T = 0 GeV), these are effective LO contributions. As we appraised from Figures 1 and 2, the qγ channel for the invariant mass distribution was sub-dominant (and almost negligible) in all kinematical regions except the low-mass region.

10^-5
10^-10 10^-15 Figure 3. Transverse momentum distribution of the two hardest photons for √ s = 13 TeV. We present results for the NLO QCD prediction (black dots) and the NLO QED distribution (blue dashed triangles). In the case of the NLO QED prediction, we also present its two different channels: the qq channel (red squared dots) and the qγ channel (green diamonds).
This can be regarded as a hint of what we are obtaining after the shoulder (p We present also two different ratios between the NLO QCD corrections and possible estimates of the higher-order QED corrections in Figure 4. Besides K QED/QCD , we define where δ NLO QED NLO (µ R ) denotes the NLO QED correction including running effects (as described in Section 2.1). Both K QED (µ R ) and K QED are similar for 0 ≤ p γγ T ≤ 1 TeV. It is a well-known fact that QCD transverse momentum resummation is requested to recover the reliability of the calculation and improve the description of data in the low p γγ T region [20]. The same consideration applies here for the NLO QED corrections and the resummation QED program, which will be presented in a forthcoming paper.  In order to explore additional measurements that could be more sensitive to QED corrections, we define a jet-veto algorithm to partially suppress the QCD components. In the first place, we consider only jets produced by QCD partons since we assume that photons can be unambiguously identified. Of course, from the experimental point of view, this is not completely true: highly energetic photons might decay into hadrons and they could be reconstructed as a jet. However, we claim that these effects are related to higher-order corrections beyond NLO QED, and that the jet identification is well-defined within the theoretical computation (at least at this perturbative order). So, after generating the event, we request the jet to be generated by a quark or gluon (but not by a photon) and to be observable within the detector. Then, we reject all the events with p jet T ≥ 50 GeV and | η jet |≤ 2.37. It is worth emphasizing that this definition is IR safe at this perturbative order, since all the IR-singular configurations automatically pass the cuts. (This jet-veto algorithm properly works for computing NLO corrections, since only one extra-parton is present in the final state. At NNLO and beyond, some subtleties could arise and the algorithm should be carefully redefined).
After this discussion, we consider the diphoton production for √ s = 13 TeV with an additional jet-veto in Figure 5. In the left panel, we show the effects on the different contributions of the total cross section. The NLO QCD correction is noticeably reduced beyond the veto threshold (i.e., p jet T = p γγ T = 50 GeV), as well as the qγ channel contribution to the NLO QED correction. We can appraise a discontinuity in the distribution at p γγ T = 50 GeV because of the modification of the event selection criteria imposed by the jetveto. Of course, this is something artificial and not related to any underlying physical phenomena. Furthermore, as expected, since the O(α) corrections to the qq channel are associated to the process qq → γγγ, this process is not affected by the jet-veto and is not suppressed. In order to quantify these effects, we extend the definition of the ratios given in Equations (24) and (32) introducing where the additional superscript denotes that we are taking into account the jet-veto.
The results are shown in the right panel of Figure 5. While the low p γγ T region is similar to the one without imposing the jet-veto (see Figure 3), at large values of transverse momentum (in particular for p γγ T 600 GeV) the qq channel dominates the total cross section, and the ratios reach O(10 %). On the contrary, we notice strong effects due to the jet-veto in the p γ,S T distribution: an important suppression of the total NLO QCD corrections and the qγ channel are shown in Figure 6 (right panel). Since the NLO QED cross section is dominated by the qq channel, which is not affected by the jet-veto, an enhancement of the QED-to-QCD ratios takes place. The estimate of the mixed QCD-QED corrections is similar to the one obtained when considering the running of α. The exception to this behavior is found in the low transverse momentum region (i.e., p In Figure 7, we present our results for the transverse momentum distribution of the two hardest photons but using the distributions of the set NNPDF3.0QED. We compare the NLO QCD cross section with the NLO QED contribution, considering also in particular its two channels: qγ and qq. The main differences with the results using the LUXqed PDFs, detailed in Figure 3, are due to the qγ channel. Even if there is a disagreement among them, due to an enhancement of the qγ channel in the middle and high-energy region, we should consider the predictions obtained with NNPDF3.0QED as a qualitative estimate in that region. . Transverse momentum distribution of the two hardest photons for √ s = 13 TeV using NNPDF3.0QED PDFs. In the right panel we show the ratios between the NLO QCD corrections and the different estimates of the QED corrections, as a function of the diphoton transverse momentum p γγ T . We consider the default choice with a frozen EM coupling (black line) and an alternative one including running effects (red dots). We appraise that the NLO QED cross-sections obtained with NNPDF3.0QED are much larger than those for LUXqed in the high-energy region due to the qγ channel.
We would like to emphasize here that the purpose of the previous comparison is to highlight that different methodologies on the extraction of PDFs could have a noticeable impact of the phenomenological analysis. In particular, a purely data-driven extraction of PDFs suffers from the lack of good-quality data-points in some regions of the parameter space, which directly impact in the error propagation. The methodology applied for extracting NNPDF3.0QED is more sensible to the lack of experimental points for very high energies, which leads to huge errors in the x ≈ 1 region. The approach followed by LUXqed uses all-order theoretical predictions to constrain the photon PDF with datapoints in the low-energy region, thus reducing the fitting errors. It is worth noticing that the updated versions of NNPDF combine the methodology developed by LUXqed in order to improve the predictions for the photon PDF [68][69][70], leading to an impressive increase on the precision achieved.
Finally, we comment about two other aspects of this computation. In the first place, we explore the effects due to the merging algorithm for photons, as described in Section 2.1. For 0.05 ≤ ∆R min γγ ≤ 0.1, both the E 0 and E-schemes present significant deviations neither between these nor compared with the distributions obtained without any clustering algorithm. The precedent consideration is valid for the total cross section, as well as for all the differential distributions presented in this work. Thus, we conclude that the implementation of a clustering algorithms with the current experimental conditions has a negligible impact in the measurements, which is compatible with the experimental and theoretical uncertainties.

Conclusions
A discussion about the computation of NLO QED corrections to diphoton production was presented, putting special emphasis in the phenomenological aspects for hadron colliders. Besides the size of these contributions and the (partial) estimate of the mixed QCD-QED terms, we emphasize that keeping under control all the possible uncertainties in this process is crucial since it is the main background for many new physics searches. Taking under control the QED corrections at the LHC allows to potentially decide about the origin of future discrepancies between data and theory. Moreover, the importance of high-precision predictions becomes completely indispensable when considering the next generation of high-energy and high-luminosity experiments [6][7][8][9][10][11][12][13].
In this article, NLO QED corrections were obtained through the application of the QED version of the q T -subtraction method, after a proper Abelianization of the NLO QCD implementation available in the program 2γNNLO. It is important to notice that we carefully studied the cancellation of IR singularities in the q T → 0 limit, in order to guarantee the theoretical consistency of the approach.
In general, we found small corrections for the M γγ distribution, whilst the ratio K QED reaches O(10 %) for the p γγ T spectrum. This is because only the real-emission terms contribute to the cross-section for p γγ T > 0: a simple power counting shows that the QEDto-QCD ratio behaves like α/α S ≈ 10 %. Besides that, we studied the effects due to the treatment of the EM running coupling, which could also introduce percent-level effects in all the distributions considered.
An interesting observation is that the qγ channel is not always negligible, specially in some kinematical regions for distributions involving the transverse momentum. For instance, it dominates the NLO QED corrections to the p γγ T distribution for moderate values of transverse momentum, i.e., p γγ T 80 GeV, as we showed in Figure 3. From this results, we conclude that it is always recommended to include the photon-initiated processes when computing differential distributions.
Besides that, using the approach described in Ref. [66], we provided an estimate of the mixed QCD-QED corrections. Other interesting effects were found when considering the dependence on p γγ,S T (the transverse momentum of the second hardest photon), exhibiting O(10 %) corrections relative to the NLO QCD contributions. These effects are enhanced when including a jet-veto to suppress the QCD component, although the M γγ distribution was slightly modified.
Finally, we would like to highlight that QED corrections cannot be neglected in the context of precision high-energy physics. Moreover, such computations must be done within a fully consistent framework, which carefully includes all the ingredients without introducing additional sources of uncertainties.
Author Contributions: Implementation of the code, L.C.; formulation of the problem, L.C. and G.S.; writing, L.C. and G.S. All authors have read and agreed to the published version of the manuscript.