Evidence for light-by-light scattering and searches for axion-like particles in ultraperipheral PbPb collisions at $\sqrt{s_\mathrm{NN}} =$ 5.02 TeV

Evidence for the light-by-light scattering process, $\gamma\gamma\to\gamma\gamma$, in ultraperipheral PbPb collisions at a centre-of-mass energy per nucleon pair of 5.02 TeV is reported. The analysis is conducted using a data sample corresponding to an integrated luminosity of 390 $\mu$b$^{-1}$ recorded by the CMS experiment at the LHC. Light-by-light scattering processes are selected in events with two photons exclusively produced, each with transverse energy E$_\mathrm{T}^{\gamma}>$ 2 GeV, pseudorapidity $ |\eta^{\gamma}|\lt $ 2.4, diphoton invariant mass $m^{\gamma\gamma}\gt$ 5 GeV, diphoton transverse momentum $p_\mathrm{T}^{\gamma\gamma}\lt$ 1 GeV, and diphoton acoplanarity below 0.01. After all selection criteria are applied, 14 events are observed, compared to expectations of 11.1 $\pm$ 1.1 (theo) events for the signal and 4.0 $\pm$ 1.2 (stat) for the background processes. The excess observed in data relative to the background-only expectation corresponds to a significance of 4.1 standard deviations, and has properties consistent with those expected for the light-by-light scattering signal. The measured fiducial light-by-light scattering cross section, $\sigma_\mathrm{fid} (\gamma\gamma \to \gamma\gamma)=$ 120 $\pm$ 46 (stat) $\pm$ 28 (syst) $\pm$ 4 (theo) nb, is consistent with the standard model prediction. The $m_{\gamma\gamma}$ distribution is used to set new exclusion limits on the production of pseudoscalar axion-like particles, via the $\gamma\gamma \to$ a $\to \gamma\gamma$ process, in the mass range $m_{\mathrm{a}} =$ 5-90 GeV.


Introduction
Elastic light-by-light (LbL) scattering, γ γ → γ γ , is a pure quantum mechanical process that proceeds, at leading order in the quantum electrodynamics (QED) coupling α, via virtual box diagrams containing charged particles (Fig. 1, left). In the standard model (SM), the box diagram involves contributions from charged fermions (leptons and quarks) and the W ± boson. Although LbL scattering via an electron loop has been indirectly tested through the high-precision measurements of the anomalous magnetic moment of the electron [1] and muon [2], its direct observation in the laboratory remains elusive because of a very suppressed production cross section proportional to α 4 ≈ 3 × 10 −9 . Out of the two closely-related processes-photon scattering in the Coulomb field of a nucleus (Delbrück scattering) [3] and photon splitting in a strong magnetic field ("vacuum birefringence") [4,5]-only the former has been clearly observed [6]. However, as demonstrated in Ref. [7], the LbL process can be experimentally observed in ultraperipheral interactions of ions, with impact parameters larger than twice the radius of the nuclei, exploiting the very large E-mail address: cms -publication -committee -chair @cern .ch. fluxes of quasireal photons emitted by the nuclei accelerated at TeV energies [8]. Ions accelerated at high energies generate strong electromagnetic fields, which, in the equivalent photon approximation [9][10][11], can be considered as γ beams of virtuality Q 2 < 1/R 2 , where R is the effective radius of the charge distribution. For lead (Pb) nuclei with radius R ≈ 7 fm, the quasireal photon beams have virtualities Q 2 < 10 −3 GeV 2 , but very large longitudinal energy (up to E γ = γ /R ≈ 80 GeV, where γ is the Lorentz relativistic factor), enabling the production of massive central systems with very soft transverse momenta (p T 0.1 GeV). Since each photon flux scales as the square of the ion charge Z 2 , γ γ scattering cross sections in PbPb collisions are enhanced by a factor of Z 4 5 × 10 7 compared to similar proton-proton or electron-positron interactions.
Many final states have been measured in photon-photon interactions in ultraperipheral collisions of proton and/or lead beams at the CERN LHC, including γ γ → e + e − [12 -21], γ γ → W + W − [22 -24], and first evidence of γ γ → γ γ reported by the ATLAS experiment [25] with a signal significance of 4.  citation denoted by the ( * ) superscript) survive the interaction and escape undetected at very low θ angles with respect to the beam direction ( Fig. 1, left). The dominant backgrounds are the QED production of an exclusive electron-positron pair (γ γ → e + e − ) where the e ± are misidentified as photons (Fig. 1, centre), and gluon-induced central exclusive production (CEP) [26] of a pair of photons ( Fig. 1, right).
The γ γ → γ γ process at the LHC has been proposed as a particularly sensitive channel to study physics beyond the SM. Modifications of the LbL scattering rates can occur if, e.g. new heavy particles, such as magnetic monopoles [27], vector-like fermions [28], or dark sector particles [29], contribute to the virtual corrections of the box depicted in Fig. 1. Other new spin-even particles, such as axion-like particles (ALPs) [30] or gravitons [31], can also contribute to the LbL scattering continuum or to new diphoton resonances. In addition, light-by-light cross sections are sensitive to Born-Infeld extensions of QED [32], and anomalous quartic gauge couplings [33].
We report a study of the γ γ → γ γ process, using PbPb collision data recorded by the CMS experiment in 2015 at √ s NN = 5.02 TeV and corresponding to an integrated luminosity of 390 μb −1 . A comparison of exclusive diphoton and dielectron yields, with almost identical event selection and reconstruction efficiencies, is provided as a function of key kinematic variables to check of the robustness of the analysis. The ratio of the LbL scattering to QED e + e − production cross sections is reported, so as to reduce the dependence on various experimental corrections and uncertainties. Using the measured m γ γ distribution, new exclusion limits are set on ALP production, in the mass range m a = 5-90 GeV.

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel (EB and HB) and two endcap (EE and HE) sections. Forward calorimetry (HF), based on a steel absorber and quartz fibres that run longitudinally through the absorber and collect Cherenkov light, primarily from the electromagnetic particles, complements the coverage provided by the barrel and endcap detectors up to pseudorapidity |η| = 5.2. Muons are measured in gas-ionisation detectors embedded in the steel flux-return yoke outside the solenoid. The silicon tracker measures charged particles within the |η| < 2.5 range. It consists of 1440 silicon pixel and 15 148 silicon strip detector modules. For nonisolated charged particles in the transverse momentum range 1 < p T < 10 GeV and |η| < 1.4, the track resolutions are typically 1.5% in p T and 25-90 (45-150) μm in the transverse (longitudinal) impact parameter [34]. The first level of the CMS trigger system [35], Level-1 (L1), composed of custom hardware processors, uses information from the calorimeters and muon detectors to select the most interesting events in a fixed time interval of less than 4 μs. The high-level trigger (HLT) processor farm further decreases the event rate before data storage. A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [36].

Simulation and reconstruction
The light-by-light signal is generated with the MadGraph v5 [37] Monte Carlo (MC) event generator, with the modifications discussed in Refs. [7,38] to include the nuclear photon fluxes and the elementary LbL scattering process. The latter includes all quark and lepton loops at leading order, but omits the W boson contributions, which are only important for diphoton masses m γ γ > 2m W .
Next-to-leading order (NLO) quantum chromodynamics and QED corrections increase σ γ γ →γ γ by just a few percent [39] and are also neglected here. Exclusive γ γ → e + e − events can be misidentified as LbL scattering if neither electron track is reconstructed or if both electrons undergo hard bremsstrahlung. This QED process is generated using the starlight v2.76 [40,41] event generator, also based on the equivalent-photon fluxes. Since the cross section for the QED e + e − background is four to five orders of magnitude larger than that for LbL scattering, and it relies on physics objects (electrons) that are very similar to those of the signal (photons), the exclusive dielectron background is analysed in depth in order to estimate many of the (di)photon efficiencies directly from the data, as well as to determine an LbL/(QED e + e − ) production cross sections ratio with reduced common uncertainties. The central exclusive production process, gg → γ γ , is simulated using superchic 2.0 [42], where the computed proton-proton cross section [26] is conservatively scaled to the PbPb case by multiplying it by where A = 208 is the mass number of lead and R g ≈ 0.7 is a gluon shadowing correction in the relevant kinematic range [43], and where the rapidity gap survival probability, encoding the probability to produce the diphoton system exclusively without any other hadronic activity, is assumed to be 100%. Given the large theoretical uncertainty of the CEP process for PbPb collisions, the absolute normalisation of this MC contribution is directly determined from a control region in the data, as explained later. All generated events are passed through the Geant4 [44] detector simulation, and the events are reconstructed with the same software as for collision data. The simulation describes the tracker material budget with an accuracy better than 10%, as established by measuring the distribution of reconstructed nuclear interactions and photon conversions in the tracker [34,45]. Photons and electrons are reconstructed using an algorithm based on the particle flow global event description (GED) [46]. The GED algorithm uses information from each subdetector system to provide charged-particle tracks, calorimeter clusters, and muon tracks. Electromagnetic showers from photons and electrons deposit 97% of their incident energy into an array of 5×5 ECAL crystals. The tracker material can induce photon conversion and electron bremsstrahlung and, because of the presence of the strong CMS solenoidal magnetic field, the energy reaching the calorimeter is thereby spread in φ. The spread energy is captured by building a cluster of clusters, or "supercluster" [47]. The GED algorithm allows for an almost complete recovery of the energy of the photons and electrons, even if they initiate an electromagnetic shower in the material in front of the ECAL. Nonetheless, in the case of photons, in order to keep to a minimum the e ± contamination, we require them to be unconverted in the tracker. The reconstructed energy of this supercluster is used to define the energy of the photon. Since the default CMS photon reconstruction algorithm is optimised for γ and e ± with transverse energies E T = E sin θ > 10 GeV, whereas the cross section for photons and electrons from exclusive production peaks in the lower E T ≈ 2-10 GeV range, a version of the GED algorithm optimised for this transverse energy range is employed. The threshold for the energy of photons, electrons, and superclusters, is lowered to 1 GeV, instead of the 10-15 GeV threshold used in the standard CMS analyses [47]. The full analysis is independently repeated using a different "hybrid" photon/electron reconstruction algorithm [47], obtaining reconstruction efficiencies and final results fully consistent with those derived with the default GED approach. Additional particle identification (ID) criteria are applied, in order to remove photons (mostly) from high-p T neutral pion decays, based on a shower shape analysis that requires the width of the electromagnetic shower along the η direction to be below 0.02 (0.06) units in the ECAL barrel (endcap).
Electron candidates are identified by the association of a charged-particle track from the collision vertex with superclusters in the ECAL. The association takes into account energy deposits both from the electron and from bremsstrahlung photons produced during its passage through the inner detector. Additional electron identification criteria (isolation, number of tracker hits, HCAL/ECAL energy deposit) are applied, as discussed in Ref. [13]. The electron energy scale is verified using a sample of γ γ → e + e − events, comparing the energy of the supercluster E to the momentum of the track p. The electron E/p ratio is within 5% of unity in the barrel and 15% in the endcaps. A good agreement is found between data and simulation, both in the energy scale and resolution. In addition, the LbL simulation is also used to validate the photon energy scale. The reconstructed supercluster energy and generated photon energy agree within a few percent, confirming that the reconstructed supercluster energy is well calibrated.

Event selection and background estimation
The exclusive diphoton candidates are selected at the trigger level with a dedicated L1 algorithm that requires at least two electromagnetic clusters (L1 EG) with E T above 2 GeV and at least one of the HF detectors with total energy below the noise threshold. No additional selection is applied in the HLT. Data are also recorded with single-photon triggers with E T thresholds above 5 and 10 GeV, and used in this analysis to estimate the efficiency of the first trigger via a tag-and-probe procedure [48], as described below. In the offline analysis, events are selected with exactly two photons, each with E T > 2 GeV and |η| < 2.4, that satisfy further selection requirements described below. Neutral and charged exclusivity selection criteria are applied to reject events having any additional activity over the range |η| < 5.2. First, all events with reconstructed charged-particle tracks with p T > 0.1 GeV are removed from further analysis. Second, events are required to have no activity in the calorimeters, above energy noise thresholds (ranging between about 0.6 GeV in the barrel, to 4.9 GeV in the HF), outside a region η < 0.15 and φ < 0.7 in the barrel ( η < Fig. 2. Acoplanarity distribution of exclusive e + e − events measured in data (circles), compared to the expected QED e + e − spectrum in the starlight MC simulation (histogram), scaled as described in the text. The curve shows a χ 2 fit to the sum of two exponential distributions corresponding to exclusive e + e − plus any residual (nonacoplanar) background pairs. Error bars around the data points indicate statistical uncertainties, and hashed bands around the histogram include systematic and MC statistical uncertainties added in quadrature. The horizontal bars around the data symbols indicate the bin size. 0.15 and φ < 0.4 in the endcap) around the two photons. The noise thresholds are determined from no-or single-bunch crossing events. To eliminate nonexclusive backgrounds, characterised by a final state with larger p T and larger diphoton acoplanarities, A φ = (1 − φ γ γ /π ), than the back-to-back exclusive γ γ events, the transverse momentum of the diphoton system is required to be p γ γ T < 1 GeV, and the acoplanarity of the pair to be A φ < 0.01. The chosen values of the pair p T and acoplanarity selections, similar to those originally suggested in Ref. [7], are motivated by previous CMS studies of exclusive dilepton production [12][13][14]. The two dominant exclusive background sources potentially remaining in the LbL scattering signal region, γ γ → e + e − and CEP gg → γ γ , are studied in detail next.

QED e + e − background
In order to have a full control of the QED e + e − background in the LbL scattering signal region, the same analysis carried out for the LbL events is done on exclusive dielectron candidates, applying the same criteria as described above for diphoton events, with the exception that exactly two opposite-sign electrons are reconstructed, instead of exactly two photons, and no additional track with p T > 0.1 GeV should be present in addition to the two tracks corresponding to the electrons. Fig. 2 shows the acoplanarity distribution measured in QED e + e − events passing all selection criteria compared to the MC expectation.
The curve is a binned χ 2 fit of the data to the sum of two exponential functions representing the exclusive QED e + e − production plus any residual background in the high-acoplanarity tail. In the region of acoplanarity below 0.01, 9570 dielectron events are reconstructed with a purity of P = 0.960 ± 0.002 (stat), obtained from the ratio of amplitudes of the two exponential functions fitted to the data. The yellow histogram shows the same distribution obtained directly from the starlight MC simulation, scaled to the total number of events in data, multiplied by the purity. The corresponding kinematic distributions of the selected γ γ → e + e − events in the A φ < 0.01 region are shown in Fig. 3, together with the corresponding MC predictions normalised in the same manner. The hashed band around the MC histograms include the systematic uncertainties (trigger, electron reconstruction and identification, and MC statistical uncertainties added in quadrature) discussed in Section 6, estimated as a function of electron E T and η. A good data-to-simulation agreement is found, thereby confirming the quality of the electromagnetic particle reconstruction, and of the exclusive event selection criteria, as well as of the MC predictions [7,41] for exclusive particle production in ultraperipheral PbPb collisions at the LHC. Small systematic differences between the central values of the exclusive dielectron data and the MC prediction are seen in tails of some of the distributions (at increasing acoplanarity and p T ) due to the presence of slightly acoplanar events in data, likely from γ γ → e + e − events where one (or both) electrons radiate an extra soft photon, that are not explicitly simulated by the MC event generator. These small discrepancies have no impact on the final extracted cross sections integrated over the whole range of distribution(s).
The QED dielectron background is then directly estimated from the starlight MC simulation by counting the number of such e + e − events that pass all LbL scattering selection criteria. The charged exclusivity condition, requiring no track in the event above the p T = 0.1 GeV threshold, is successful in removing this background almost entirely. This tracking efficiency is controlled in events with exactly two reconstructed photons and exactly one track, finding good data-MC agreement. The QED background in the signal region is estimated to be N ee,data = 1.0 ± 0.3 (stat) events, where the assigned uncertainty corresponds to the event count in the simulated samples.

Central exclusive diphoton background
Although the LbL and CEP processes share an identical final state, their kinematic distributions are different. Diphotons from quasireal γ γ fusion processes are produced almost at rest in the transverse plane and, thus, the final-state photons are emitted back-to-back with balanced pair transverse momentum p On the other hand, typical CEP photon pairs are produced in diffractive-like gluon-mediated processes [26,42] with larger momentum exchanges leading to a diphoton transverse momentum distribution peaking at p T ≈ 0.5 GeV, after selection criteria, and moderately large tails in the azimuthal acoplanarity distribution.
Thus, the requirement on diphoton acoplanarity (A φ < 0.01) also significantly reduces the gg → γ γ background. Since the MC prediction for CEP gg → γ γ has large theoretical uncertainties, and in order to account for any other remaining backgrounds resulting in photons that are not back-to-back (such as γ γ → e + e − γ (γ ) events passing the analysis selection criteria), for which no event generator is currently available, the CEP background is normalised to match the data in the region A φ > 0.02, where the contribution from γ γ → γ γ is negligible (Fig. 4). The background normalisation factor is then obtained from f norm (1) and found to be f norm nonacoplanar = 1.06 ± 0.35 (stat). The number of events due to CEP plus any residual backgrounds is thus estimated to be 3.0 ± 1.1 (stat). The statistical uncertainties quoted in both values are driven by the size of the data sample remaining at high acoplanarities, after all selection criteria have been applied. Table 1 Number of diphoton candidates measured in data and expected from MC simulation for LbL scattering, QED e + e − production, and from the CEP+other contributions, after each event selection step (cumulative) described in the text. The yields of the simulated processes are scaled according to their theoretical cross sections and the integrated luminosity of the analysed data set. The CEP+other values are normalised from the high-acoplanarity tail with a scale factor estimated from the data as described in the text. The LbL scattering simulation uncertainty quoted is that of the theoretical uncertainty of the prediction, whereas the uncertainties in the QED e + e − and CEP+others yields are statistical.

Selection criteria Data
LbL MC QED e + e − MC CEP MC+other (normalised to data) Charged exclusivity 648 Neutral exclusivity 108

Light-by-light signal distributions
The exclusive diphoton signal is extracted after applying all selection criteria described above and estimating the amount of residual QED e + e − and CEP+other backgrounds. Table 1 shows the number of events remaining after each selection criterion. The main selection requirement corresponds to two photons each with E T > 2 GeV, |η| < 2.4 (excluding photons falling in the η ≈ 0.1 gap region between the EB and EE, 1.444 < |η| < 1.566), and diphoton invariant mass greater than 5 GeV. The numbers of events measured in data and expected from the sum of all MC contributions in the first two rows do not match because these selection requirements accept a fraction of nonexclusive backgrounds that are not included in the simulation. Once the full exclusivity selection criteria are applied, the data-to-simulation agreement is very good. We observe 14 LbL scattering candidates, to be compared with 9.0 ±0.9 (theo) expected from the LbL scattering signal, 3.0 ± 1.1 (stat) from central exclusive plus any residual diphoton backgrounds, and 1.0 ± 0.3 (stat) from misidentified QED e + e − events.
An extra selection criterion has been also studied by further requiring that the candidate LbL scattering events have no signal above the noise threshold in the pixel tracker layers. This more stringent selection is sensitive to charged particles down to ∼40 MeV, and results in a number of reconstructed LbL scattering signal counts (and even more reduced QED backgrounds) consistent with the MC predictions. However, since the efficiency of such a tight selection is difficult to assess from a control region in data, the default analysis is kept with the charged-particle track p T > 0.1 GeV exclusivity requirement.

Cross section extraction
Given the low signal yield available for an extraction of differential cross section distributions, an integrated fiducial cross section for LbL scattering above a diphoton mass m γ γ = 5 GeV is calculated instead. The ratio R of cross sections of the light-bylight scattering over the QED e + e − processes is measured, thereby reducing the uncertainties related to trigger and reconstruction efficiencies, and integrated luminosity. Efficiency uncertainties partially cancel in the ratio, as described later, thanks to a similar selection applied to photons and to electrons; and the integrated luminosity dependence fully cancels out. The ratio R is defined as Here σ fid (γ γ → γ γ ) is the LbL scattering fiducial cross section (i.e. passing all the aforementioned p T , η, m γ γ kinematic selection criteria for the single photons and for the photon GeV) is the total cross section for the QED e + e − process for masses above 5 GeV; Acc ee = N gen (p gen T is the dielectron acceptance for the fiducial single-electron kinematic selections determined from the starlight MC generator; N γ γ ,data is the number of diphoton events passing the selection in data; N γ γ ,bkg is the estimated number of background events passing all selection criteria; N ee,data is the number of dielectron events passing our selection in data; P is the purity of the estimated fraction of QED e + e − signal among these dielectron events; and C γ γ and C ee are the overall efficiency correction factors, for the γ γ and e + e − selections, respectively, that are determined as discussed in the next section.

Diphoton analysis efficiencies
The C γ γ correction factor in Eq. (2) is obtained through the factorised expression where the diphoton efficiency ε γ γ is determined using the LbL scattering MC simulation. This efficiency receives contributions from triggering, photon reconstruction and identification, and neutral and charged exclusivity criteria that are directly determined from the data via independent data-to-simulations scale factors, SF = ε data /ε MC , as explained below.
The diphoton efficiency is first derived from the LbL scattering simulation via: where the selection in the numerator and denominator applies to exactly two photons required in each event, which are also within the fiducial kinematic region in diphoton p T , mass, and acoplanarity. It is found to be ε γ γ = (20.7 ± 0.4)%, mostly driven by the inefficiencies of the single photon reconstruction and identi-fication, and of the trigger (ε γ ,reco+ID , ε γ γ ,trig. ≈ 70%). The quoted uncertainty here is statistical only, reflecting the finite size of the LbL scattering MC sample. The second term of Eq. (3), the photon reconstruction and identification efficiency correction (SF γ ,reco+ID ), is extracted from data by selecting γ γ → e + e − (γ ) events, where one of the electrons emits a hard bremsstrahlung photon due to interaction with the material of the tracker. The p T of the two electrons in γ γ → e + e − events being approximately equal, if one of the electrons emits a hard bremsstrahlung photon, it may not reach the ECAL to be identified as an electron but it can still be reconstructed in the tracker as a charged particle. In a first step, hard-bremsstrahlung events are selected among events passing a trigger requiring one L1 EG cluster with E T > 5 GeV, that have exactly two oppositely charged particle tracks and exactly one electron reconstructed. Among those events, we then look for exactly one photon compatible with a hard bremsstrahlung, as described below. Such events are used to estimate the efficiency in a tag-and-probe procedure, via ε γ ,reco+ID, hard-brem data = N reco+ID,hard-brem probe N reco+ID, hard-brem passing , (5) where denominator and numerator are defined as follows: • N reco+ID, hard-brem passing : Electrons are selected if (i) their direction matches with one of the two reconstructed tracks within a radius R = √ ( η) 2 + ( φ) 2 < 1.0 (where η and φ are those of the electron track), (ii) they have E T above 5 GeV, and (iii) their associated ECAL supercluster is matched within R < 0.1 to an L1 EG cluster with E T > 5 GeV. The p T of the track that is not matched with the electron should be below 2 GeV, since we assume that track to be generated by the electron after bremsstrahlung emission. The p unmatched track T < 2 GeV requirement ensures that this low-p T charged particle is sufficiently bent by the magnetic field, and thus the expected photon (extrapolated to the ECAL) and the second electron are sufficiently separated. Events entering the denominator are not required to have a reconstructed photon.
• N reco+ID, hard-brem probe : Events from the denominator are also included in the numerator if a photon is found with E T > 2 GeV that passes the identification criteria.
The efficiency is extracted using a fit to the acoplanarity distribution between the electron and the charged-particle track, and amounts to ε γ ,reco+ID, hard-brem data = (86.5 ± 7.0)%, to be compared with ε γ ,reco+ID, hard-brem MC = (82.5 ± 2.0)% in the MC simulation, where uncertainties are statistical (as well as all other uncertainties quoted in this section). The ratio of these efficiencies is used to define the corresponding SF γ ,reco+ID, hard-brem = 1.05 ± 0.09 scale factor. We note that this procedure checks not only the reconstruction and identification efficiency in data, but also effectively includes bin migrations outside the fiducial p T range due to the effects of photon energy scale and resolution. The impact of bin migrations in the final diphoton cross section is found to be below the 1% level.
Events in the study above comprise exactly two chargedparticle tracks, corresponding to the two electrons. They do not probe the possibility that the photon, reconstructed in the ECAL, has previously also interacted in the tracker generating an e + e − pair that has been also reconstructed as one or two additional displaced low-p T charged-particle tracks. In the case of an LbL scattering event, such a genuine signal event would be discarded by the strict charged exclusivity criterion, which is applied independently of the proximity of the tracks to the photon, in order to keep the QED e + e − background to a minimum. The modeling of this efficiency loss in simulation is checked using hard-bremsstrahlung events with a similar selection as above, except that up to two additional charged-particle tracks are now allowed in the event. We check the fraction of events where no additional track, more displaced than the one with p T < 2 GeV required in the selection, is found in a window | η| < 0.15, | φ| < 0.7 around the photon. This efficiency amounts to ε γ ,tk veto data = (89.9 ± 1.7)% in data, and ε γ ,tk veto MC = (91.1 ± 1.2)% in the MC QED e + e − simulation. The ratio of these efficiencies gives SF γ ,tk veto = 0.99 ± 0.02. The final overall scale factor for reconstruction and identification, accounting for the modeling of photon conversions in the MC simulation and the efficiency to reconstruct the associated displaced tracks, is then SF γ ,reco+ID = SF γ ,reco+ID, hard-brem SF γ ,tk veto = 1.04 ± 0.09.
The third term of Eq. (3) accounts for the trigger selection efficiency. Exclusive diphoton events are selected using an L1 trigger requiring two electromagnetic clusters with E T > 2 GeV, and no activity (above noise thresholds) in at least one of the HF calorimeters. These two components of the trigger, the electromagnetic cluster selection and the HF energy veto, are verified independently in data. The efficiency for reconstructing an L1 EG cluster with E T > 2 GeV is verified using a tag-and-probe technique on QED e + e − events, where the dielectron acoplanarity is fitted to extract the signal and measure the efficiency. The same selection criteria used in the main analysis are applied, including the exclusivity requirements. Events are further selected using a supporting trigger requiring one L1 EG cluster with E T > 5 GeV with the same HF energy veto as the analysis trigger. The L1 EG cluster used in the trigger is matched (using the same criterion mentioned above) to one of the two electrons reconstructed offline, called the tag. The other electron in the event is the probe, and it qualifies as a passing probe if it is matched to an L1 EG cluster with E T > 2 GeV.
The efficiency is then the fraction of probes that are also passing probes, and it is in the 45-100% range, with the lowest efficiency found close to the E T = 2 GeV threshold and at high |η|. Scale factors are determined from the data-MC differences, as a function of E T , in two |η| bins. Applying them to the LbL simulation, we find an integrated scale factor of 1.12 ± 0.31 (stat). The same QED e + e − sample is used to test the HF veto component of the analysis trigger. This time, we apply the nominal dielectron selection, including exclusivity requirements, but for a data sample collected with a trigger requiring a single-EG object in the HLT with E T > 10 GeV and |η| < 1.5 plus a small amount of energy in the HFs corresponding to about 50% of the most peripheral PbPb events. Both electrons in the event are then matched to an L1 EG cluster with E T > 2 GeV. We find that (100 +0 −3 )% of the selected events also pass the analysis trigger, i.e. satisfy the HF veto in the trigger, in perfect agreement with the result predicted from the MC simulation. Combining the results of the studies above, the scale factor SF γ γ ,trig. = 1.12 ± 0.31 is obtained for the ratio of the product of trigger efficiencies in data to that obtained from the MC simulation.
The last two terms of Eq. (3) account for the efficiency of the exclusivity selections. The fraction of events passing the QED dielectron selection, with the exception of the charged and neutral exclusivity criteria, are analysed. Using the acoplanarity distribution to extract the signal, we find that (92.5 ± 0.3)% of the events feature no additional track in the event, to be compared to (99.3 ± 0.1)% in simulation. We deduce that the corresponding scale factor is SF ch.excl. = 0.93 ± 0.01. A similar strategy is used for the neutral exclusivity selection, this time in events passing the corresponding requirements. This efficiency is found to be (89.9 ± 1.4)% in data, and (96.9 ± 1.3)% in simulation. This scale factor is then SF neut.excl. = 0.93 ± 0.02. Differences between the exclusivity efficiencies in data and MC simulation are likely due to the presence of nonexclusive events, such as γ γ → e + e − processes with a small hadronic overlap of the lead ions, whose modeling is currently not available in the Monte Carlo generators. The incorporation of such nonexclusive events in the definition of the signal is irrelevant, because both SF ch.excl. and SF neut.excl. cancel out in the R ratio, as explained in Section 6.

Dielectron analysis efficiencies
For the exclusive dielectron analysis, the efficiency is estimated using the starlight MC simulation via

Table 2
Summary of the overall cross section measurement efficiencies C γ γ ,ee , efficiencies from simulation ε γ γ ,ee , and individual data-to-simulation scale factors SF γ γ ,ee , obtained for the diphoton and dielectron analyses. "Reco. and ID" stands for reconstruction and identification. All quoted uncertainties are systematic. where the kinematic criteria in the numerator and denominator are applied to exactly the two electrons required in the event. The different components of the electron efficiency are again checked using data, via a factorised expression for the corresponding correction factors: Most of the scale factors are common with those used in the diphoton analysis (since they are computed using the larger statistical sample of electrons in data), except for the reconstruction and identification efficiency, which we check again for electrons. For the latter, a tag-and-probe technique using a fit to the acoplanarity distribution in QED e + e − events is used, as done for the diphoton case, except that now the probe is a charged-particle track that is a passing probe if it is matched to an electron passing the reconstruction and identification criteria. We find an efficiency of (89.4 ± 1.2)% in data, consistent with (90.4 ± 1.3)% in the MC simulation, corresponding to a scale factor of SF e, reco+ID = 0.99 ± 0.02.
The scale factor for the trigger efficiency is also recomputed using the p T spectrum in the QED e + e − MC simulation, using the same p T -and |η|-dependent scale factors as for SF γ γ ,trig. , leading to SF ee,trig. = 1.09 ± 0.16.

Summary of the efficiencies
The overall cross section measurement efficiencies, efficiencies in simulation, as well as the individual data-to-simulation scale factors, obtained for the diphoton and dielectron analyses are summarised in Table 2. Since the data-to-simulation scale factors are consistent with unity, they are not included in the numbers listed in Table 1 nor in the results plotted in Figs. 2-5, but they are used to obtain the results in Section 7. The overall diphoton cross section efficiency, Eq. (4), is about 20% compared with about 10% for dielectrons, Eq. (6). The dielectron analysis is a factor of two less efficient than the diphoton one, because each single electron has a relatively larger probability of losing energy by bremsstrahlung before reaching the ECAL, and therefore their probability to pass the trigger selection threshold and/or their energy be properly reconstructed is smaller. Such efficiency losses are further enhanced as they enter squared for two electrons to pass the trigger or be concurrently reconstructed above the p T and mass thresholds.

Systematic uncertainties
The main sources of uncertainty in the LbL scattering and QED e + e − production measurements are related to the trigger and single γ , e ± reconstruction efficiencies ( Table 2). The uncertainty in Table 3 Summary of the systematic uncertainties in the ratio of the fiducial LbL scattering to total QED e + e − cross sections. the latter is doubled in the total cross section, since we consider diphoton and dielectron final states. No additional uncertainty in the photon energy scale and resolution is considered, since possible data-simulation differences are already included in the derivation of the reconstruction and identification scale factors. Systematic uncertainties have been estimated for the different terms defining the ratio R of the LbL scattering over QED e + e − production cross sections given by Eq. (2), and are summarised in Table 3. Because the scale factors used for the trigger efficiency are common to the diphoton and dielectron analyses, their associated uncertainty cancels partially in the ratio. However, because of different reconstruction and identification efficiencies, the E T spectrum of photons is different from that of the electrons, leading to only incomplete cancellation of the uncertainty. Assuming the uncertainty in each individual p T -binned scale factor is fully correlated between the photon and electron cases, but with no correlation between the scale factors for different p T , we propagate these uncertainties simultaneously to the numerator and denominator of the ratio SF γ γ ,trigger /SF ee,trigger , resulting in a 12% uncertainty in that ratio. The charged and neutral exclusivity scale factors, common to the diphoton and dielectron measurements, are assumed to cancel in the ratio R. The rest of the SF terms listed in Table 2 have small statistical uncertainties from the finite size of the MC samples used to derive them, which propagate into percent uncertainties in the final cross section, and are neglected here. Among the other parameters in Eq. (2), the normalisation of the CEP and QED e + e − backgrounds in the signal region propagates into a 16% uncertainty in the background yield (resulting in an 6% uncertainty in the cross section measurement), accounting for the finite size of the MC samples. An additional uncertainty of 25% (10% in the final cross section), reflecting the finite size of the data sample at high acoplanarities used for the absolute normalisation of the CEP plus residual nonexclusive backgrounds, is considered as a statistical uncertainty rather than a systematic one.
The final systematic uncertainty is obtained from adding in quadrature the individual uncertainties and is listed in the last row of Table 3.

Light-by-light cross section
The compatibility of the data with the background-only hypothesis has been evaluated from the measured acoplanarity distribution (Fig. 4), using a profile-likelihood ratio as a test statistic, including all systematic uncertainties as nuisance parameters with log-normal priors [49,50]. The uncertainty due to the finite size of the MC samples is also included as an additional nuisance parameter for each bin of the histogram. The significance of the excess at low diphoton acoplanarity in data, estimated from the expected distribution of the test statistic for the background-only hypothesis obtained with pseudo-experiments, is 3.7 standard deviations (3.5 standard deviations expected). If using only the total number of events observed and expected in the region A φ < 0.01, we obtain a significance of 3.4 standard deviations (3.2 expected).

Exclusion limits on axion-like particle production
The measured invariant mass distribution (Fig. 5, center right) is used to search for possible narrow diphoton resonances, such as pseudoscalar axion-like particles produced in the process γ γ → a → γ γ [30]. The LbL, QED, and CEP+other continuum processes are considered as backgrounds in this search. Fully simulated starlight samples for various ALP masses, m a , ranging from 5 to 90 GeV are reconstructed with the same code used for the LbL analysis in order to estimate the ALP acceptance and efficiency, as well as the expected reconstructed diphoton mass template distributions. Corrections to the efficiency estimated in the MC simulation are derived based on data, and applied in the same way as for the LbL analysis. A binned maximum likelihood fit of the signal and background contributions is performed on the data, where systematic uncertainties are included as nuisance parameters with a log-normal prior. The CL s criterion [53,54], with a profile likelihood ratio as test statistic [55], is used to extract exclusion limits in the σ (γ γ → a → γ γ ) cross section at 95% confidence level (CL). Limits on σ (γ γ → a → γ γ ) cross section for axion-like particles with masses 5-90 GeV are set in the 1500-20 nb range (Fig. 6). The 68 and 95% CL bands around the expected limits are obtained using pseudo-experiments. The cross section limits shown in Fig. 6 are used to set exclusion limits in the g aγ vs, m a plane, where g aγ ≡ 1/ is the ALP coupling to photons (with being the energy scale associated with the underlying U(1) symmetry whose spontaneous breaking generates the ALP mass). Two scenarios are considered where the ALP couples to photons F μν alone, or also to hypercharge B μν with operators: a F F /4 and aB B/(4 cos 2 θ W ) (where θ W is the Weinberg angle), respectively [30]. The derived constraints on the ALP mass and its coupling to photons are compared in Fig. 7 to those obtained [30,56]

Summary
Evidence for light-by-light (LbL) scattering, γ γ → γ γ , in ultraperipheral PbPb collisions at a centre-of-mass energy per nucleon pair of 5.02 TeV has been reported, based on a data sample corresponding to an integrated luminosity of 390 μb −1 recorded by the CMS experiment at the LHC in 2015. Fourteen LbL-scattering candidate events passing all selection requirements have been observed, with photon transverse energy above 2 GeV and pseudorapidity |η| < 2.4, diphoton invariant mass greater than 5 GeV, diphoton transverse momentum lower than 1 GeV, and diphoton acoplanarity below 0.01. Both the measured total yields and kinematic distributions are in accord with the expectations for the LbL scattering signal plus small residual backgrounds that are mostly from misidentified exclusive dielectron (γ γ → e + e − ) and gluon-induced central exclusive (gg → γ γ ) processes. The observed (expected) significance of the LbL scattering signal over the background-only expectation is 3.7 (3.5) standard deviations. The ratio of the fiducial LbL scattering to the total QED dielectron cross sections is R = (25.0 ± 9.6 (stat) ± 5.8 (syst)) × 10 −6 .
From the theoretical γ γ → e + e − cross section prediction, we derive a fiducial light-by-light scattering cross section, σ fid (γ γ → γ γ ) = 120 ± 46 (stat) ± 28 (syst) ± 12 (theo) nb, consistent with the standard model expectation. The measured exclusive diphoton invariant mass distribution is used to set new exclusion limits on the production of pseudoscalar axion-like particles (ALPs), via the process γ γ → a → γ γ , over the m a = 5-90 GeV mass range. Fig. 7. Exclusion limits at 95% CL in the ALP-photon coupling g aγ versus ALP mass m a plane, for the operators a F F /4 (left, assuming ALP coupling to photons only) and a B B/(4 cos 2 θ W ) (right, including also the hypercharge coupling, thus processes involving the Z boson) derived in Refs. [30,56] from measurements at beam dumps [60], in e + e − collisions at LEP-I [56] and LEP-II [57], and in pp collisions at the LHC [13,58,59], and compared to the present PbPb limits.
For ALPs coupling to the electromagnetic (and electroweak) current, the derived exclusion limits are currently the best over the m a = 5-50 GeV (5-10 GeV) mass range.

Acknowledgements
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centres and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses.