Light-by-Light Scattering at Next-to-Leading Order in QCD and QED

The recent experimental observation of Light-by-Light (LbL) scattering at the Large Hadron Collider has revived interest in this fundamental process, and especially of the accurate prediction of its cross-section, which we present here for the first time at Next-to-Leading Order (NLO) in both QCD and QED. We compare two radically different computational approaches, both exact in the fermion mass dependence, thus offering a strong cross-check of our results. The first approach is a fully analytic method to calculate compact and well-organized two-loop helicity amplitudes. The second one is entirely numerical and leverages the Local Unitarity construction. Our two calculations agree with each other and conclude that including the exact fermion mass contribution typically increases the size of the NLO corrections. Moreover, we find that the exact result converges slowly to the massless limit of the high-energy regime, thus emphasizing the importance of including the full mass dependence at NLO. We also compare our results with the ATLAS measurement of LbL in ultra-peripheral lead-lead collisions, and find that the inclusion of exact NLO corrections reduces, but does not eliminate, the existing tension with theoretical predictions.


Introduction
Light, a form of ElectroMagnetic (EM) radiation, has revolutionized our understanding of Nature, e.g., from Maxwell's equations to the establishments of the two pillars of modern physics -Einstein's theory of relativity and quantum mechanics.A natural question about light that needs answering is how it interacts with itself.The abelian nature of Quantum ElectroDynamics (QED) prevents any direct interaction, and due to Furry's theorem [1], the three-photon interaction vanishes in QED and is strongly suppressed in ElectroWeak (EW) theory [2,3].Therefore, the four-photon vertex is the most viable light self-interaction, and it can be directly probed by studying the Light-by-Light (LbL) scattering process γγ → γγ.LbL scattering predictions date back to the 1930s [4][5][6][7].However, its experimental evidence remained elusive and only recently emerged at the Large Hadron Collider (LHC) through Ultra-Peripheral heavy-ion Collisions (UPCs) [8][9][10][11], where the impact parameter is larger than the sum of the two ion radii and the ions scatter quasi-elastically.The experimental feasibility of observing LbL has revived broad theoretical interest in studying the process in the context of both the Standard Model (SM) and beyond the SM (BSM) of particle physics.In the literature, LbL has been suggested to be an ideal probe of the bound states of leptons [12,13] and gluons [14], the quartic anomalous gauge couplings [15], the axion-like particles [16], the graviton-like particles [17,18], the nonlinear Born-Infeld extensions of QED and SM [19], the photon self-interaction from the noncommutative QED [20], large extra dimensions [21,22], and supersymmetry [14].
The direct experimental observation of γγ → γγ was only achieved recently in heavy-ion UPCs at the LHC [8][9][10][11], thanks to the large coherent photon flux carried by the ultra-relativistic nucleus with a large charge number Z (e.g., the lead Pb ion has Z = 82).The cross-section in lead-lead (Pb-Pb) collisions benefits from a huge enhancement factor of Z 4 ≈ 4.5 • 10 7 com- 1 Strictly speaking, the ions do not need to be intact since they can be EM excited and emit neutrons.These forward neutrons can be tagged.However, in this paper, we only consider the inclusive case with respect to the forward neutrons.Therefore, the word "intact" here should be understood in a loose sense.duction for obtaining compact two-loop helicity amplitudes.The second approach considers the Local Unitarity (LU) construction [29,30] to directly compute the fully differential cross-section entirely numerically.We have cross-checked the results from these two methods against each other and found perfect agreement.We stress that this is the first time that LU, or any fully numerical method in momentum space, yields results for a yet unknown cross-section.Our results bring the accuracy of the theoretical prediction for the LbL crosssection down to the percent-level in the regime where the di-photon invariant mass ranges from a few GeV up to the W ± mass threshold 2m W .

Methodology
As already mentioned, our first approach is based on the derivation of the analytic one-and two-loop helicity amplitudes for the process γ(p 1 , λ 1 ) + γ(p 2 , λ 2 ) + γ(p 3 , λ 3 ) + γ(p 4 , λ 4 ) → 0, where all external four momenta p i are incoming and λ i 's are the helicities of the photons.Up to any loop order in QCD and QED, we can decompose the helicity amplitude onto the following tensor basis: where ⃗ λ = (λ 1 , λ 2 , λ 3 , λ 4 ), ε λ i ,µ i are the photon polarization vectors, and the coefficients A i , B i jk and C i jkl are functions of the Mandelstam variables s = (p 1 + p 2 ) 2 , t = (p 2 + p 3 ) 2 , and u = (p 1 + p 3 ) 2 , as well as the masses of the internal fermions.Taking advantage of the transversality of external photons, Bose symmetry and gauge invariance makes it possible to rewrite the original 138 form factors into only 7 independent ones A 1 (s, t, u), A 1 (t, s, u), A 1 (u, s, t), ∆B 1  11 (s, t, u), ∆B  1), all with implicit arguments (s, t, u), have also been rewritten into form factors with permutations of the arguments.The 5 independent helicity amplitudes only depend on the following 5 form factors A S (s, t, u), ∆ B1 11 (s, t, u), ∆ B1 11 (t, s, u), ∆ B1 11 (u, s, t), and ∆ Ĉ2111 (s, t, u) as follows: 11 (s, t, u) = ∆B 1 11 (s, t, u) + A 1 (s, t, u)/u, and ∆ Ĉ2111 (s, t, u) = ∆C 2111 (s, t, u) − A S (s, t, u)/(su).The form factors A S (s, t, u), su∆ Ĉ2111 (s, t, u), and the allplus and one-minus amplitudes M ++++ and M −+++ are fully symmetric in s, t, u.The last three amplitudes can be related to each other using crossing symmetry: All other helicity amplitudes can be obtained from these 5 independent ones using charge conjugation C, parity P, and time reversal T symmetries.
Our work only considers the two-loop QCD and QED amplitudes mediated through a fermion loop.The twoloop contribution mediated by the W ± boson would require the inclusion of the complete set of EW corrections, which is beyond the scope of this letter.

Analytical approach
For a given fermion species, the two-loop amplitudes, generated with QGRAF [32] and FEYNARTS [33], are reduced to a linear combination of master integrals using integration-by-parts identities as implemented in FINITEFLOW [34] and KIRA [35,36], in d = 4 − 2ϵ dimensions.In order to solve the master integrals analytically, we transform them into a canonical basis [37] that satisfies the ϵ-form of differential equations [38] with uniform transcendental weight.The solution can then be generally expressed through iterated integrals, that can be expressed either as multiple polylogarithms or retained as iterated integrals with dlog one-forms.Both of these representations allow fast numerical evaluations.Some of these integrals have also been studied in the context of other processes [37,[39][40][41][42].In addition, to simplify the final expression of the helicity amplitudes, we use the shuffle algebra of iterated integrals and the symbol techniques [43,44] to find relations among master integrals at a given transcendental weight.The rational coefficients are further simplified using FINITEFLOW and MULTIVARIATEAPART [45] as well as taking advantage of the symmetry properties of the amplitudes.The final two-loop amplitudes renormalized in the On-Shell (OS) scheme are short enough to be written in a few pages [31], which allows us to analytically verify pole cancellations, and also compute interesting kinematic limits.

Numerical approach
Our second computational approach relies on direct MC integration in momentum space using the LU construction.The starting point of this method is to consider all Forward-Scattering Graphs (FSG) relevant to a particular cross-section, and for each of them collect all Cutkosky cuts contributing to the scattering process definition of interest.For the LbL process, each FSG contains a single Cutkosky cut traversing two photons and separating the FSG into a two-loop (at NLO) diagram on the left of the cut, and a one-loop diagram on the right of the cut, thus forming one particular interference term contributing to 2ℜ M 2-loop M 1-loop ⋆ .To simplify our computation, we chose to write the 1-loop amplitude as an effective four-photon vertex with the exact 1-loop amplitude M µ 1 µ 2 µ 3 µ 4 as its Feynman rule.At LO (NLO), this leads to 2 ( 16) distinct non-isomorphic two-(three-)loop FSGs, see fig. 1 for an example FSG contributing at NLO. Loop integrals are computed numerically, together with the phase-space integral stemming from the integration of the momenta entering the Cutkosky cut and accompanied by an observable function.LU achieves this by projecting momenta onto the Cutkosky cut using the causal flow [29] and by writing loop integrals using the Loop-Tree Duality (LTD) [46][47][48][49][50][51][52] expression, obtained by analytic integration over the energy components of the loop momenta.In this work, the manifestly causal LTD expression [53] was used.Ultra-Violet (UV) singularities are locally subtracted in a way that immediately yields results renormalized in the OS scheme [30].Threshold singularities appear when √ s > 2m f and they are regularized using a contour deformation of the spatial loop momenta [54].To mitigate large gauge cancellations between FSGs, we wrote the photon spin-sums in the temporal axial gauge with n = (1, 0, 0, 0).The resulting 9-dimensional momentumspace integral is then evaluated using MC methods with importance sampling.

Results
We first set the notation used in the presentation of our results.For a given fermion f with mass m f we denote the one-and two-loop QCD and QED amplitudes as respectively.For the W ± boson, since we only consider the one-loop result, we denote its amplitude as M (0,0,W) ⃗ λ . The amplitude combinations that will be used in the cross-sections are then M (0,0) . Then, the LO partonic crosssections can be computed with where dΦ 2 is the 2-body phase-space measure, and the overlined sum denotes the sum over the helicity configurations of the external photons and the average over the initial photon polarizations.We can also define the NLO QCD and/or QED corrections to the partonic cross-sections as follows ) with (i, j) ∈ {(1, 0), (0, 1), (1, 1)}.Thus, the NLO QCD, NLO QED, and NLO QCD+QED cross-sections are σNLO QCD = σ(0,0) + σ(1,0) , σNLO QED = σ(0,0) + σ(0,1) , and σNLO QCD+QED = σ(0,0) + σ(1,1) respectively.We have analogous definitions for a specific fermion contribution by adding f in the superscripts.Thanks to the exclusiveness of the LbL process2 , we can even include partial (5) which is the prediction we provide here.
The physical LbL cross-sections of AB γγ → AγγB, with beam particles A and B, should convolve the partonic cross-section with the corresponding photon-photon flux L (AB) so that σ AB = dx 1 dx 2 L (AB) (x 1 , x 2 ) σ, where x 1 and x 2 are the longitudinal momentum fractions carried by the initial-state photons.Note that the photon-photon flux cannot simply be factorized into a product of two photon density functions of A and B in the case of UPCs [27] because of the presence of a non-trivial survival probability.In the special case of A = B = γ, we take with δ being the Dirac delta function.
In order to assess the mass effect, without losing generality, let us first consider the NLO QCD corrections to σ γγ with a single massive top quark loop of m f = 173 GeV.We fix the α value in the Thomson limit α(0) = 1/137.036and the strong coupling constant α s = 0.118.We also impose a transverse momentum cut on the final photons p γ T > √ s/100.In the center-of-mass frame of the two initial-state photons, this cut corresponds to limiting the scattering angle to be within 1.15 • ≲ θ γ ≲ 178.85 • , or equivalently the pseudo-rapidity to be |η γ | ≲ 4.6.The NLO QCD cross-section in terms of √ s/m f from the analytic amplitudes is given as the black curve in the upper panel of fig. 2. We also overlay with blue points the results from the numerical LU approach, with its corresponding MC uncertainty.Each of these results with LU has been obtained using 100M sample points across all FSGs, computed in 50 CPU-hours to reach below 1% in the NLO correction.The two approaches agree well within the errors except in the asymptotic limits √ s ≪ m f and √ s ≫ m f , where numerical instabilities prohibit a fair comparison.We also compare our full mass-dependent results with the two known approximations: the highenergy (HE) [24] (red) and the LE [23] (green) limits.Our exact results match well with these two approximations in their applicable regimes.In particular, the relative difference between the exact computation and the HE (LE) approximation reaches below 2% when √ s/m f > 37 ( √ s/m f < 0.52).A kink appears at the threshold √ s → 2m f , where the two-loop amplitudes suffer from the Coulomb singularity.Such a sin-gularity is however integrable and thus harmless when convolved with the realistic photon-photon flux that exhibits a wide spectrum.A proper treatment of this region requires a Coulomb resummation, which is beyond the scope of this letter.The K-factors of σ NLO QCD γγ /σ LO γγ are displayed in the lower panel of fig. 2. Our exact computation exhibits a non-trivial shape of the K-factor as a function of the collision energy, whereas the HE and LE approximations are close to the constants 1.124 and 1.512 respectively.In contrast to the LE limit, the exact mass-dependent K-factor approaches the HE limit slowly.
The fermion f loop: We now collect all SM contributions from the charged fermions and the W ± boson, and perform a data-theory comparison for the ATLAS Pb-Pb UPC measurement at √ s NN = 5.02 TeV [11] with our stateof-the-art computation.Since we have renormalized the internal fermion masses in the OS scheme, we set pole masses to the following values: m µ = 0.1134 GeV, m τ = 1.777GeV, m c = 1.5 GeV, m b = 4.75 GeV, m t = 172.69GeV and m W = 80.377 GeV, whereas the remaining SM charged fermions-electron, up, down and strange quarks -are kept massless.The EM coupling α associated to the external on-shell photons is set to α(0) = 1/137.036,but we consider the renormalization group running of the couplings α s and α in the loops.For the central value of the renormalization scale, we choose µ R,0 = m γγ /2.We then estimate missing higher order contributions by varying the renormalization scale by the factors 1/2 and 2. The scale evolution of the strong coupling α s (µ R ) follows the standard 2-loop modified minimal subtraction scheme with variable flavor number and α s (m Z ) = 0.118.The running of the coupling α(µ R ) is obtained from the photon propagator at momentum transfer µ R , with α(m Z ) = 1/127.955.We also considered the alternative EW G µ renormalization scheme for α and verified that it yields insignificant differences.For the Pb-Pb UPC photon-photon flux L (PbPb) , we consider the charge form factor (ChFF) available in GAMMA-UPC [27], and we tested that the resulting K-factors are mostly independent of the choice of flux but of course not the resulting cross sections.
The NLO ′ QCD+QED cross-section within the AT-LAS fiducial volume increases by 6.5% +2.1% −1.2% with respect to (wrt) the LO cross-section of 76 nb3 .Our final prediction of 81.2 +1.6 −0.9 nb4 stands 1.8σ below the ATLAS measurement of 120 ± 22 nb, which is a reduction of the tension observed when comparing to LO predictions.The quoted uncertainty is obtained from a (1/2, 2)-variation of the renormalization scale.For reference, we stress that the LE and HE approximations for the K-factors of the tau, charm quark and bottom quark contributions increase the LO cross-section by 13% and 0.7% respectively (see bottom layout of fig.3).This significantly differs from our prediction and thus highlights the importance of retaining exact fermion mass dependence in the computation of NLO corrections to the LbL cross-section.
Fig. 3 reveals that the tension between data and theory is largest in the first di-photon invariant mass bin m γγ ∈ [5,10] GeV.This observation motivates the studies of the impact from the C-even bottomonia [55] and the first experimentally-observed fully-charmed tetraquark state X(6900) [56].We have also simulated 6 Ceven bottomonium states as well as X(6900) with the HELAC-ONIA event generator [57,58].The inclusive cross-sections of these resonances are proportional to the square of their di-photon decay widths (see eq. ( 7) in [27]), which are mostly unconstrained by experiments so that we take their values from theoretical cal-culations [56,59].In particular, the X(6900) Breit-Wigner peak shown in fig. 3 is generated using the decay width Γ X(6900)→γγ computed assuming the vector meson dominance hypothesis.It turns out that the contributions of these resonances to the LbL cross-section are negligible.The study in ref. [56] reveals bridging the gap between the LO prediction for the LbL cross-section and ATLAS data using solely the X(6900) resonance, one would need to increase Γ X(6900)→γγ by two orders of magnitude wrt their calculation assuming vector meson dominance.Our state-of-the-art NLO ′ QCD+QED result (shown as the grey hatched band in fig. 3) reduces but does not eliminate this tension.It is interesting to note that the size of NLO QCD and QED corrections is largest in the lowest m γγ bin and decreases from 10% down to 2% in the highest mass bin.This can be explained by noting that the tau, charm, and bottom mass effects are expected to be larger at lower di-photon invariant masses, as it can be seen from the HE and LE approximations of these fermion contributions shown in the lower panel of fig. 3. The exact mass-dependent Kfactor is substantially different than these approximations.Indeed, the HE approximation strongly underestimates quantum corrections for small m γγ , whereas the LE approximation significantly overestimates them, especially for larger values of m γγ .The comparison of our NLO cross-section with ATLAS data for other observables can be found in Appendix B, as well as a more detailed comparison between the results from our numerical and analytical approaches in Appendix A. The additional fiducial cut on the acoplanarity distribution due to the non-zero but small virtuality of the initial photons further reduces our theoretical predictions by around 0.5%, which we have ignored here.

Conclusion
In this letter, we have carried out for the first time the computation of NLO QCD and QED contributions to the LbL cross-section, retaining exact dependence in the fermion masses.We considered two completely different and independent computational approaches that serve as a strong cross-check of our results.The first method is more traditional and results in compact analytical expressions for the two-loop helicity amplitudes, whereas the LU method is fully numerical and was introduced only recently.This stands as the first yetunknown cross-section obtained with LU and it highlights its potential as an alternative for tackling cuttingedge multi-loop computations.We have also compared our exact result to known approximations in the high and low energy regimes, and we find that they are only satisfactory in limited regions of the phase-space that do not cover the entire range of phenomenological interest.In particular, we find that our exact result converges slowly to the high-energy limit.Therefore, we conclude that retaining exact fermion mass dependence in the computation of higher-order corrections to the LbL cross-section is mandatory for reliable predictions.Finally, we compare our final prediction to the ATLAS measurement of LbL in Pb-Pb UPCs and find that the inclusion of exact NLO QCD and QED corrections reduces, but does not eliminate, the moderate tension observed for this process.
Acknowledgements.VH thanks Charalampos Anastasiou and Zeno Capatti for their insights regarding the transient singularity appearing in the LU integrand for massless internal fermions, and Ben Ruijl for his help with the symbolic treatment of the LU numerator.HSS is grateful to David d'Enterria for pointing out the impact of the acoplanarity cut on our results.This work is supported by the grants from the ERC (grant 101041109 'BOSON', grant 101043686 'LoCoMotive'), the SNSF (grant PCEFP2 203335), the French ANR (grant ANR-20-CE31-0015 'PrecisOnium'), the French LIA FCPPL, the European Union's Horizon 2020 research and innovation program (grant 824093 STRONG-2020, EU Virtual Access 'NLOAccess').

Appendix A. Comparisons between analytical and numerical results obtained with LU
In this section, we provide a more detailed comparison of the differential results obtained using both our analytical and numerical approaches.We turn off the convolution with the incoming photon flux and consider LbL cross-sections for three choices of fixed collision energy corresponding to distinct kinematic regimes: a) below the 2m f threshold at m f / √ s = 0.5767, b) when mass effects are important at m f / √ s = 0.173, and c) when the HE approximation starts becoming viable at m f / √ s = 0.0865.We consider two fully correlated observables: the transverse momentum p T of the final-state photons, and the cosine of the angle from their momenta wrt the beam cos θ.We report our findings in fig.A.4 (LO, σ (0,0, f ) γγ ) and fig.A.5 (NLO correction only, excl.LO contribution, σ (1,0, f ) γγ ), and we find perfect agreement across the range displayed for both observables, with relative errors on the fiducial cross-sections at, or below, 0.1% and pulls rarely exceeding 1σ at LO and 2σ for the NLO correction.
We note that numerical instabilities in our implementation of the analytic one-loop effective four-photon form factors in LU prevented us from reaching values much below p T / √ s = 0.1 (both at LO and NLO).For this reason, the cos θ range was also adjusted to be within [−0.98, 0.98].This could easily be addressed by promoting their implementation to higher arithmetic precision.
For massless internal fermions, the numerical approach of LU suffers from a transient infrared singularity in FSG featuring internal self-energy diagrams.A transient integrand singularity is non-integrable but yields no dimensional regularization pole at the integrated level.The origin of this transient singularity is well-understood and could be cured by either tensorreducing the internal self-energy or including additional Cutkosky cuts traversing the massless internal fermion loop, leading to the final-state of γ f f which can in principle not be distinguished from γγ in the collinear limit for massless fermions.In the absence of such a treatment, the presence of this transient singularity implies poor numerical convergence of the LU construction in the HE limit, and in practice, this renders the numerical computation impractical for √ s/m f > 15.We however stress that the inclusion of the additional γ f f cuts is natural in light of the final experimental resolution, and that tensor-reduction of internal self-energy can be achieved with minimal effort.Analogously, the implementation of the analytic two-loop amplitudes with exact fermion mass dependence suffers from numerical inaccuracies in the LE and HE limits that limit its range of applicabilities in particular phase space corners, but which could be alleviated using higher-precision arithmetics.
Finally, we point out the stark change in shape of the photon p T distribution for m f / √ s = 0.0865 wrt the LO distributions and also NLO correction distributions at lower energies.

Appendix B. Details on the comparison of our prediction to data
In this section, we present additional comparisons between our prediction and the ATLAS measurement [11] for other observables in fig.B.6.From left to right, the figure shows the single photon (averaged) transverse momentum p γ T , the absolute value of the rapidity of the photon pair y γγ , and the absolute value of the cosine of the pair polar angle cos θ γγ .θ γγ is the scattering angle in the rest frame of the initial two photons.In each plot, the layout is similar to the one of fig. 3. The K-factor decreases from 10% at the lowest p γ T value to around 1% at its highest.In contrast, the y γγ and cos θ γγ distributions exhibit rather flat quantum corrections, however with noticeable trends.Analogous to what we observed in fig.3, the HE (LE) approximation underestimates (overestimates) the size of the exact quantum corrections in all bins.
We now briefly explain how we constructed the LE and HE approximations in figs.3 and B.6 presented in this letter.In the HE limit, we take all fermions massless except for the top quark with m t = 172.69GeV and the W ± boson with m W = 80.377 GeV in both LO and NLO ′ QCD+QED, in order to obtain the Kfactors.Conversely, the LE limit is obtained by setting m τ = 1.777GeV, m c = 1.5 GeV, m b = 4.75 GeV, m t = 172.69GeV, m W = 80.377 GeV and m i = 0 for i ∈ {e, u, d, s, µ}.The NLO ′ LE partonic cross section is then defined to be: , (B.1) where and M (0,0, f ) ⃗ λ respectively.It is thus clear that the LO LE partonic cross-section is identical for the transverse momentum p T and the cosine of the angle-with-beam cos θ of the final-state photons.The result using the square of the analytical one-loop amplitudes is shown in green, and the fully numerical LU result is shown in blue.The latter was obtained using 1B sample points across both contributing 2-loop FSGs, with a run time of approximately 250 CPU hours per choice of √ s.
to the exact cross-section σ(0,0) .Note that these definitions of the LE and HE limits imply that the dependence on the masses of the top quark and the W ± boson are always accounted for exactly, although we stress that their contribution to the cross-section presented is small.(excl.LO contribution) to the LbL process at fixed incoming energies for the transverse momentum p T and the cosine of the angle-with-beam cos θ of the final-state photons.The result using the analytical two-loop amplitudes is shown in green, and the fully numerical LU result is shown in blue.The latter was obtained using 1B sample points across all 16 contributing 3-loop FSGs, with a run time of approximately 500 CPU hours for each choice of √ s.

Figure 1 :
Figure 1: Example of one of the 16 distinct 3-loop FSG contributing to the NLO correction of the LbL cross-section.The single Cutkosky cut contributing is shown in red.The effective four-photon vertex is denoted with a cross and is implemented with the exact 1-loop amplitude.The double line corresponds to a massive fermion.

Figure 2 :
Figure2: The NLO QCD LbL cross-section (upper) as well as the K factor (lower) for a given quark contribution.The black curve and the blue points denote the full mass calculations with the analytic and numeric approaches.The red and green curves represent the HE and LE approximations.

Figure 3 :
Figure 3:  The comparison between the LO (blue dashed), NLO ′ QCD+QED (grey hatched) with the ATLAS measurement[11] for the di-photon invariant mass distribution.The C-even resonances are also displayed in the upper panel, while the K factors from the two approximations are reported in the lower panel.

>2. 5 Figure B. 6 :
Figure B.6: The data-theory comparisons for the transverse momentum of the photon p γ T (left), the rapidity of the photon pair y γγ (middle) and the cosine of the scattering angle in the di-photon rest frame cos θ γγ (right).